劉應興的部落格

本文要談的問題是：假設 k 個物件，如商品品牌、旅遊地之類的，有一單變量指標影響個體（消費者、旅遊者、評價者）的偏好選擇，結果會呈現怎樣的偏好？這裡假設選擇是隨機的，因此「偏好」最後是以機率來表現。

以 k = 2 為例，假設個體在選擇時是對兩群體獨立評量，Xi ～ N(μ_i, σ^2_i), i = 1, 2。評量的過程是群體 i 取得觀測值 Xi，而後取觀測值較大的為其最後選擇：不是群體（物件） 1 就是群體 2。

p = P{選擇物件 1} = P[X1 > X2] = Φ([(μ1-μ2)/√(σ^2_1+σ^2_2)])

因此，由 p 可以反推 μ1 和 μ2 的差距：

μ1 - μ2 = z(p) √(σ^2_1 + σ^2_2)

此處 z(p) 正式的符號是 Φ^{-1}(p), 而 Φ 是標準常態分布的（累積）分布函數 (cumulative distribution function, CDF)。如果對兩物件的評估不是獨立的，也就是說 X1 與 X2 不獨立，假設它們服從雙變量常態分布，則

p = Φ([(μ1 - μ2)/√(σ^2_1 - 2ρ σ_1 σ_2 + σ^2_2)])

或其反向：μ1 - μ2 = z(p) √(σ^2_1 - 2ρ σ_1 σ_2 + σ^2_2), 式中 ρ 是 X1 與 X2 的相關係數。由 p 反推 μ1 - μ2 適用於我們有 p 的估計值，而希望換成「計量的」或數值的差異量。但如果沒有辦法估計 σ_1, σ_2 與 ρ 的數值，只能以標準值 z(p) 呈現。

如果  k = 3，如果兩兩比較，有三對比較，上面的方式可以算出

p(ij) = 1 - p(iji) =Φ((μ_i - μ_j)/√(σ^2_i - 2ρ(i,j) σ_i σ_j + σ^2_j))

如果三物件並呈，令 Y_1 = X_i - X_t, Y_2 = X_j - X_t，則

P[X_i > X_j > X_t] = P[Y_1 > Y_2 > 0]

其中 Y_1, Y_2 聯合服從雙變量常態分布，

E[Y_1] = μ_i - μ_t = δ_1,    E[Y_2] = μ_j - μ_t = δ_2
Var[Y_1] = σ^2_i + σ^2_t - 2 σ(it),    Var[Y_2] = σ^2_j + σ^2_t - 2 σ(jt)
Cov[Y_1, Y_2] = σ(ij) - σ(it) - σ(jt) + σ^2_t

假設 X_1, X_2, X_3 相互獨立，則

P[Xi > Xj > Xt]
   = ∫_(-∞, ∞) ∫_(-∞, x_i) ∫_(-∞, x_j) f_i(x_i)f_j(x_j)f_t(x_t) dx_t dx_j dx_i

在常態性的假設下，不論上面的直接積分，或是化成 Y_1, Y_2 如前列 P[Y_1 > Y_2 > 0]，都不是一個簡單的計算式。由此可知，若是直接 k > 2 個物件的排序，其機率計算是很複雜的事，那麼如果偏好資料的蒐集是採用 k 物件直接排序，或以三個一組的方式，如果想以常態模式轉換成數值指標，是很困難的。

假設有 k 物件，給 n 個案以成對比較法做 C(k,2) = k(k-1)/2 對的排序，又如何定出 k 物偏好順序，如何轉成數值？由 k = 2 的情形可知實際上由 p(ij) = P[Xi>Xj] 只能得 μ_i - μ_j 而不能得到 μ_i 與 μ_j，即使現在有 C(k,2) 對的比較，也不可能得到諸 μ_i。就單一個案評價而言，物件 i 與 j 比較，不是選 i 就是選 j，以 s(p,i,j) = 1 表示個案 p 選 i 而 0 表示選 j，而 s(p,i,j) + s(p,j,i) = 1，故

s(p,i,.) = Σ_j s(p,i,j) = 個案 p 評判物件  i 優於其他物件數

類似地，s(p,.,j) = Σ_i s(pi,j) 是個案 p 評削優於物件 j  的物件數。加總式 Σ_j 是對所有 j≠i 加總；類似地，Σ_i 是對‵所有 i≠j 加總，共 n-1 項。所以 s(p,i,.) 愈大表示物件 i 愈被偏愛，s(p,.,j) 愈小則是物件 j 愈受喜愛。

s(p,.,i) = Σ_j (1 - s(p,i,j)) = (k-1) - s(p,i,.)

所以 s(p,i,.) 可視為個案 p 對物件 i 的喜好評分，而 s(p,.,i)+1 則是倨案 p 對物件 i 的喜好排序，

Σ_i s(p,i,.) = Σ _i s(p,.,i) = ΣΣ_{i≠j} s(p,i,j) = ΣΣ_{i<j} (s(p,i,j)+s(p,j,i)) = k(k-1)/2

另一方面，把 s(p,i,j) 對 p 加總得

n(i,j) = Σ_p s(p,i,j) = 物件 i 與 j 的比較中偏好 i 的樣本數

而 n(i,j) + n(j,i) = n 為個案數，故 p(i,j) = n(i,j)/n 為物件 i 與 j 的比較中偏好 i 的樣本比例。樣本比例 p(i,j) 可轉換成 z(p(i,j))，當 σ^2_i ≡ 1  且假設個案評價物件 i 與 j 時所以據的背後變數 X_i, X_j 相互獨立，則 μ_i - μ_j 的估計是 Z(ij) = √2 z(p(i,j))，如果 p(i,j) 是區間 (0,1) 之間均勻分布的隨機變數，則 Z(ij) 將服從常態分布；不過此處 p(i,j) 並非 U(0,1) 隨機變數，而是二項分布隨機變數之緊縮版：n(i,j)/n，一個二項分布隨機變數除以一個常數，只在 n 足夠大時其分布接近常態但變異數又小，使得 Z(ij) 也接近常態分布，並且

Var[Z(ij)] ≒ (2/φ^2(z(p(i,j))))p(i,j)(1-p(i,j))/n

注意其變異數和 E[Z(ij)] ≒ √2 z(p(i,j)) 有關；對於樣本比例，消弭變異數隨期望值而變的「變異數穩定化變換 (variance stabilizing transformation)」應是「（平方根）反正弦變換 (arcsine transform)」

z = sin^{-1}(√p*),   p* = (n(i,j)+3/8)/(n+3/4)

上式是一個修正版本，原始版本是直接對 √p(i,j) 做 arcsine 變換。

如前述，將 p(i,j) 轉換成 Z(i,j) = √2 z(p(i,j)) 是估計 μ_i - μ_j。但總共有 k(k-1)/2 對「獨立的」觀測，分別估計的 μ_i - μ_j 將會有不一致的情形，也就是說估計結果將不會滿足

μ_i - μ_k = ( μ_i - μ_j) + ( μ_j - μ_k)

因此，不妨把它看成線性模型

Z(ij) =  μ_i - μ_j + ε(ij),  i < j

我們不知個案在評估諸物件 i 對比物件 j 的偏好時是否完全相互獨立，但既然考慮兩兩比對的評量方式，似乎可以假設在做諸對間評比時互不影響，故上列線性模型可以假設諸 ε(ij) 間是相互獨立的，所以問題只是諸 ε(ij) 之變異數如前述是和 p(i,j) 相關連的。不過如果我們只想得到諸 μ_i 的一個可用估計，似乎（誤差項）變異數不等並不是一個問題：我們無需分析誤差，也不求估計之最優性，只想得到一組 μ_i 的值，足以代表諸個案偏好之彙總指標。而對於上列線性模型 μ_i 不可估，只能估得諸  μ_i - μ_j 及其組合式函數式，可任意定個原點，例如 min{μ_i, i = 1,..., k} = 0 或 Σ_i μ_i = 0 即可解決。由這些 μ_i 的一致性結果，可以再反算回

p*(i,j) = Φ(( μ_i - μ_j)/√2)

成為以偏好機率表示的成對偏好評估。

用常態分布描述群體分布只是關連偏好比例和數值性差異的一種假設模型，我們沒有理由不允許群體是另一種分布，例如 Gumbel 分布。若 X_i, i = 1, ..., k, 相互獨立服從

f_i(x) = (1/β) e^{-[(x-μ_i)/β + e^{-(x-μ_i)/β}]},   -∞ < x < ∞

其平均數為 μ_i + γβ，其中 γ 是 Eular 常數。則 X_i - X_j 服從 logistic 分布，其平均數 δ = μ_i - μ_j，而其 p.d.f. 為

g(x) = e^{-(x - δ)/β}/{β [1 + e^{-(x - δ)/β}]^2,   -∞ < x < ∞

分布函數 F(x) = 1/[1 + e^{-(x - δ)/β}], 故

P[X_i > X_j] = P[X_i - X_j > 0] = e^{(μ_i - μ_j)/β}/[1 + e^{(μ_i - μ_j)/β}]

令 p(i,j) = P[X_i > X_j]，則 logit(p(i,j)) = ㏑(p(i,j)/p(j,i)) = (μ_i - μ_j)/β，反過來

μ_i - μ_j = β logit(p(i,j)) = ㏑(p(i,j)/p(j,i))

或 p(i,j)/p(j,i) = e^{μ_i}/e^{μ_j}。令 π_i = e^{μ_i}, 則

p(i,j) = π_i/(π_i + π_j)

這稱為 Bradley–Terry 模型，是一個成對偏好的 logit 模型，可以用 ML（最大概似度法）由成對比較偏好資料估計諸 π_i 或設定固定 β 值估計諸 μ_i。不論 π_i 或是 μ_i 都可以當作偏好度的指標，前者不是機率，但都是正值，可把它化成似乎在 k 物件擇其一的機率；後者是數值化的偏好度指標，但沒有固定的 0 點，就像先前機率值變換 Z(ij) = z(p(i,j)) 一般只能確定 μ_i - μ_j 而 μ_i 等卻不可估。

前面兩種變換：Z(ij) = z(p(i,j)) 與 logit(p(i,j)) 明確顯示：假設評價者依據的指標 X_i 有不同分布，在成對比較中，其偏好機率 p(i,j) 可能與兩群體（兩物件評價指標）平均數之差 μ_i - μ_j 有關， 但關聯形式不同。再看一個例子，假設

X_i ～ (1/μ_i) e^{-x/μ_i},   x > 0

則 P[Xi > X_j] = μ_i/(μ_i + μ_j) = (μ_i/μ_j)/[1+(μ_i/μ_j)]。此例 p(i,j) 不再與 μ_i - μ_j 相關連而是 μ_i/μ_j = e^{㏑(μ_i)-㏑(μ_j)}。其所以如此，是因在此指數分布，期望值 μ 是尺度參數 (scale parameter)；在常態分布，期望值  μ 是位置參數 (location parameter)；在 Gumbel 分布，期望值是 μ 加上一個和 β 有關的常數，前面我們假設 β 相同，而 μ 是位置參數。就一般來說，兩相互獨立隨機變數差別在於位置參數，則其一大於另者的機率是位置參數差的函數；若其差別是在於尺度參數，則機率是尺度參數比的函數。

如同在 Bradley–Terry 模型用一組 π_i 來表示物件 i 的偏好度，在指數群體成對比較模型中，也可以固定一組數值 μ_i 代表物件 i 的偏好指標，由 p(i,j) 與 μ_i/μ_j 的關係，得

μ_i = c/[Σ_j p(j,i)/p(i,j)]

若用最大概似法，可選擇一組 μ_i 極大化

ΣΣ_{i<j} [n(ij) ㏑(μ_i) + n(ji) ㏑(μ_j) - n ㏑(μ_i + μ_j))]

式中 n(ij) 是在 n 次比較物件 i 與 j 中，偏愛 i 的次數。對 μ_i 微分，得概似方程式

(1/μ_i) Σ_{j≠i} p(i,j) = Σ_{j≠i} 1/(μ_i + μ_j)

式中 p(i,j) 是指其粗估計，n(ij)/n。故

μ_i = [Σ_{j≠i} p(i,j)]/[Σ_{j≠i}  1/(μ_i + μ_j)]

解 μ_i 需要迭代計算。

如果底層分布是如上的指數分布，我們不只可得兩兩比較的

p(i,j) = P[X_i > X_j] = μ_i/(μ_i + μ_j) = (1/μ_j)/(1/μ_i + 1/μ_j)

還可計算三三比較中特定偏號順序之機率得：

p(i,j,t) = P[X_i > X_j > X_t]
      = [(1/μ_j)/(1/μ_j + 1/μ_i)] [(1/μ_t)/(1/μ_t + 1/μ_j + 1/μ_i)]
      = (μ_i^2 μ_j)/[(μ_i + μ_j)(μ_i + μ_j + μ_t)]

甚且，此結果可推至 k 物件一次排序：

p(i_1, ..., i_k) = P[X(i_1) > ... > X(i_k)]
       = [(1/μ(i_2))/(1/μ(i_1) + 1/μ(i_2))] . . .
                   + [(1/μ(i_k))/(1/μ(i_1) + ... + 1/μ(i_k))]

實務上如採兩兩比較法蒐集樣本個案對 k 物件的偏好評價，可得 C(k,2) 個獨立的偏好比例；採三三比較，可得 6 C(k,3) 個偏好順序比例，其中 5 C(k,3) 是獨立的；要求個案做 k 物件完全排序，則有 k! - 1 獨立的樣本比例。不管哪種方式，要做的無非是統合各特定順序之樣本比例，套入以上 p(...) 等與諸 μ_i 之關係式做 μ_i 之估計，結果可用以代表樣本個案對 k 物件之總合偏好。這要如何進行？我們有 p(i,j,...) 與諸 μ_i 的關係式：

p(i,j,...) = P[X_i > X_j > ...; μ_1, ..., μ_k]

有諸 p(i,j,...) 的自然不偏估計，即樣本比例 p*(i,j,...)，選擇 μ_i 等以極小化目標

Q(p(i,j,...), p*(i,j,...))

目標函數 Q 中，p(i,j,...) 是諸 μ_i 的函數式，即使如常態等不能如指數分布那樣直接寫出其計算式的也無妨。實際計算大概需要迭代：先設定諸 μ_i 之初值，代入 Q 計算其函數值，再根據極小化 Q 所要湍足的條件修改 μ_i 的設定值，再次評估、改進，直到結果滿意或可接受。