指數族分布與廣義線性模型－劉應興的部落格

指數族分布 (exponential family (of distributions)) 和指數分布族 (exponential distributions) 是兩回事，不過也許中文名稱會搞混了。我們說指數分布（族），就像說常態分布（族），二項分布（族）之類的，指的是一類分布稱指數分布，只是其中有參數（位置參數及尺度參數）未定。「分布族」這樣的稱呼可大可小，如指數分布族、均勻分布族等，一族之不同分布只差在位置參數與尺度參數（指這兩分布族），甚至有時候只允許其中一個參數（位置或尺度）可改變；分布族一詞也可代表一個龐大的機率分布族群，如這裡談的「指數族」，其 p.d.f. 或 p.m.f. 形式如下：

f(x; θ) = C(θ) Q(x) e^{b(θ)．t(x)}, x in A, A 與 θ 無關

或完全寫成指數形式

f(x; θ) = e^{b(θ)．t(x) + c(θ) + q(x)}, x in A, A 與 θ 無關

其中 x, θ 可以是向量，t(x), b(θ) 也可以是向量，故 b(θ)．t(x) 其實是向量的點積。如前述指數分布可以是也可以不是指數族的一員：

f(x; a, b) = (1/b) e^{-(x-a)/b}, x ≧ a

如果位置參數 a 是可變的，機率密度正值的範圍依參數 a 而變，它就不符前面指數族要求 A 與 θ 無關的定義了；但若位置 a 是固定的，只有尺度參數 b 允許改變，這樣的指數分布就是指數族的一員。很多統計學上常見的分布，如離散型的二項、Poisson，連續型的常態、beta 都是指數族成員；但也有很多不是，如前述帶未定「位置」的指數、Cauchy、超幾何等，都不是指數族。

廣義線性模型 (generalized linear model) 和一般線性模型 (general linear model) 又是不同，後者就是線性模型，稱「一般」只是有別於古典的線性模型假設了

Y = X β + ε, E[ε] = 0, Cov(ε) = σ^2 I

甚至有時還假設了 ε 服從常態分布。而「一般」線性模型是假設誤差項可以有相關，可以變異數不等，因此 ε 的共變異矩陣是 σ^2 V。以單一觀測值來看（假設各觀測值獨立）：

E[Yi] = μ_i, μ_i = x_i' β, Yi ～常態

而廣義線模是假設：

E[Yi] = μ_i = g^{-1}(η_i), η_i = g(μ_i) = x_i' β, Yi ～指數族

和線模不同的是 μ_i 不是直接等於 x_i' β, 而是透過 g 的反函數與「線性預測子 (linear predictor)」連上關係。另外是反應變數 (response variable) 一般假設是指數族，而因為包含了離散型，廣義線模直接考慮反應變數的分布，不像線性模型中以誤差項的分布來描述。總而言之，廣義線模把模型分解成三個成分：

(1) 隨機成分：反應變數 Y 的分布；

(2) 線性預測子：η = x'β，解釋變數的組合；

(3) 連結函數 (link function)：線性預測子 η 和（反應）平均數 μ 的關係式，η = g(μ)。

有時還會有一個「過度分散參數」(overdispersion)。「分散」這個中文詞，有的寫「離散」，必須小心分辨它的意思：筆者儘可能用「分散」、「散布」表示隨機變數或資料分布、散布的範圍的一個指標，如變異數、標準‵差、四分位差之類的。而「離散」指的是隨機變數或資料的一種型態，是分布在一些互有間隔相互分離的點上各點有其機率，而不是像連續型變數或資料，我們通常需要以一個區間來考慮分布的機率。當然數學上離散型分布可能不像這裡說的分布在外觀上「相互分離」「互有間隔」的點，因為那些「離散」的點可能是稠密的 (dense)，不過實務上我們可能很少遇到，就不談它了。過度分散本來是

Y 似應服從某一分布，故 Var[Y] = V(μ) 其中 μ = E[Y]，
但實際上 Y 的分散度 Var[Y] > V(μ)。

與過度分散的情況相反的是「偏低分散 (underdispersion)」現象。過度分散比較常見，由條件變異數恆等式

Var[Y] = E[Var[X|Z]] + Var[E[X|Z]]

如果 Y 的分布其實依賴另一隱藏變數，如果

E[Var[X|Z]] ＝ V(E[X]) = V(E[E[X|Z]])，

則 Var[E[Z|X]] 成為額外的部分；另一方面，如果 E[Var[X|Z]] 比 V(E[X]) 小，則加上 E[X|Z] 的變異仍然可能偏低，這就產生偏低分散現象。不過，在此種混合模型，顯然較容易產生過度分散而不容易產生偏低分散的現象，特別是在 Z 條件下 Y 服從假設的模型情況，其無條件分布（邊際分布）可能不是原先所以為的分布，

回到我們的指數族: f(x; θ) = e^{b(θ)．t(x) + c(θ) + q(x)}，其中 c(θ) 對 p.d.f./p.m.f. 而言是常數；但對概似函數而言，q(x) 才是常數。由

∫_A f(x; θ) dx = 1 (或 Σ f(x; θ) = 1)

如果是單參數，對 θ 微分，則（以積分式表示）

∫_A (b'(θ) t(x) + c'(θ)) f(x; θ) dx = 0

即 b'(θ) E[t(X)] = - c'(θ)，故 E[t(X)] = - c'(θ)/b'(θ)。如果有 k 個參數，則 b(θ) 也要 k 項，而且同樣 k 項 t(x), 或說 t_j(x) 線性獨立，

Σ_j Ｄi(b_j(θ)) E[t_j(X)] = - Di(c(θ))，i = 1, ..., k.

式中 Di 代表對 θ_i 偏微。由 k 個方程式解 k 個 E[t_j(X)]。有時為了方便，會重新參數化，使 b_j(θ) = θ_j,，也會為了方便，取 t_j(x) = x_j，如在 GLM 中，常是如此。簡化的指數族成為：

f(x; θ) = e^{θ．x + c(θ) + q(x)}, x in A

也可以加上分散參數：

f(x; θ, φ) = e^{(θ．x + c(θ))/a(φ) + q(x, φ)}, x in A

對 θ_i 偏微分及二階偏微分後積分，也就是對恆等式

∫_A e^{(θ．x + c(θ))/a(φ) + q(x, φ)} dx = 1

做積分號下微分，可得

E[X_i] = - Di(c(θ)), Cov(X_i X_j) = - a(φ) Dij(c(θ))

所以 a(φ) 就是過度或偏低分散的乘數。而

μ = E[X] = - ▽ c(θ)，當 a(φ) = 1 時 Cov(X) = - ▽▽' c(θ) = ▽ μ'

運算符 ▽ 是對諸 θ_i 偏微的行向量，▽▽' 就是二階偏微矩陣，是對諸 θ_i 偏微後轉置成列向量，再對其中每一元素求對諸 θ_j 的偏導數，構成二階偏微矩陣。由於 μ 是 c(θ) 的一次偏微，Cov(X) 是 μ 對 θ 的偏微，如果 θ 可以表示成 μ 的函數，則 Cov(X) 也是 μ 的函數，寫成 Cov(X) = V(μ)。

假設 X1, ..., Xn i.i.d. 上述帶分散參數的指數族，資料是 x1, ..., xn，這裡假設是多參數多變量模型，則對數概似函數為

l(θ, φ; x) = Σ [(θ．xi)/a(φ) + c(θ) + q(xi, φ)]

對 θ 做偏微分並令其為 0，得方程式

( Σ xi + n ▽ c(θ) )/a(φ) = 0

對 φ 偏微得

-(a'(φ)/a(φ)^2) [θ．(Σ xi) + c(θ)] + Σ D_φ q(xi, φ) = 0

可以令 a(φ) = φ，這可以說不失其一般性，因為 q(x, φ) 對於 φ 的依賴也是透過 a(φ)，所以重參數化也就是令 a(φ) = φ 並沒失去什麼。

但在上列兩方程式組成立的臨界點是否得（對數）概似函數的極大，還要看 Hessian 矩陣，由二階偶導數構成的矩陣是否為負確定 (negative definite) 或負半確定 (negative semidefinite)。對數概似函數對 θ 的二階偏微是 n(▽▽'c(θ)), 一個非正確定矩陣 (nonpositive definite)，因 Cov(X) = -a(φ) ▽▽'c(θ) 是非負確定 (nonnegative definite)。也就是說：對於任意 φ，概似方程式給了 θ 的 MLE。至於 φ，在 a(φ) = φ 之下，對數概似函數對 φ 的二階偏微是

2[θ．(Σ xi) + c(θ)]/φ^3 + Σ D_φ^2 q(xi, φ)

就這一般式難以確定其正負；不過回到原（對數）概似函數，（在 θ 給定的情況下）可以用其他方式尋找其最大值及對應的 φ 值。

考慮 Yi 相互獨立服從指數族 ( φ = 1 )，E[Yi] = μ_i, η_i = g(μ_i) = β'x_i, i = 1, ..., n。則對數概似函數為

l(θ; Y) = Σ_i (θ_i Yi + c(θ_i) + q(Yi))

其中 c'(θ_i) = μ_i = g^{-1}(η_i) = g^{-1}(Σ_j x_{ij}β_j)，x_{ij} 是第 x_i 第 j 解釋變數的值，則 l(θ; Y) 對 β_j 的偏導數為

Ｄ_j l(β; Y) = Σ_i (Yi - μ_i)/(V(μ_i) g'(μ_i)) x_{ij}

首先是對 θ_i 微分，然後 θ_i 對 μ_i 微分，μ_i 對 η_i 微分，最後 η_i 對 β_j 做偏微。其中 θ_i 對 μ_i 微分是 μ_i 對 θ_i 微分的導數，後者是 -c"(θ_i) = Var[Yi] = V(μ_i)；μ_i 對 η_i 微分是 η_i 對 μ_i 微分的倒數，後者是 g'(μ_i)。故概似方程式以行向量表示：

▽ l(β) = Σ_i (Yi - μ_i)/(V(μ_i) g'(μ_i)) x_i = 0

有時對各觀測值有額外的加權 Wi，則用下列估計方程式：

Σ_i Wi (Yi - μ_i)/(V(μ_i) g'(μ_i)) x_i = 0

令 Wi(μ_i) = Wi/(V(μ_i) g'(μ_i))，並以之構成對角線矩陣 W，又令 X = [x_1 ... x_n]'，如同線性模型的 X 矩陣，Ｙ, μ 分別是諸 Yi 和 μ_i 構成的行向量，則概似方程式可寫成

X' W (Y - μ) = 0

類似線性模型用加權最小平方法得到的估計方程式，不過 μ 不是 X β, 而是 X β = g(μ)，其意是對諸 μ_i 取 g 函數轉換；另外是 W 和 μ 有關，也就和 β 有關。把上列方程式左邊寫成 β 的函數 d(β)，這是 β 的非線性函數，考慮 Newton-Raphson 迭代法

β* = β - (▽ d(β)')^{-1} d(β)

給予 β 適當初值，而後依上列公式修正成 β*，然後以之代替 β 再修正，如此持續迭代至收斂到一個難以再改進的結果。所謂「難以再改進」有兩個方式來看：一是 d(β) 接近 0，二是 β 的改進夠小。這其中又有不同的評量方式，例如 d(β) 看其成分絕對值之和或平方和；β 的改進可以看絕對值或相對值。在上列迭代式中，d(β)' 的微分結果是矩陣 X'W*X，其中 W* 是對角線 n 階方陣，其第 i 元素

Wi* = - Wi/[V(μ_i) g'(μ_i)^2]
+ Wi (Yi - μ_i) (V'(μ_i) g'(μ_i) + V(μ_i) g"(μ_i))/[V(μ_i)^2 g'(μ_i)^3]

有一個問題值得注意：上述推導是假設 a(φ) = φ = 1，如果 φ 不是 1 呢？對數概似函數對 β 偏微只是在 (Yi - μ_i) 底下多個 φ 或 a(φ) 當分母，但這是常數，在概似方程式中完全可以拿掉而不影響結果。

回顧對數概似函數對 β 的微分得

Σ_i (Yi - μ_i)/(V(μ_i) g'(μ_i)) x_i

對單一觀測值的話是

(Y - μ)/(V(μ) g'(μ)) x

如果是對 μ 微分的話得到的是 (Y - μ)/V(μ)。如果我們只知 Var(Y) = V(μ), 變異數和平均數 μ 的關係，並把 (Y-μ)/V(μ) 當做對數概似函數對 μ 的導數，則對數概似函數應為

l*(μ) = ∫ (Y-μ)/V(μ) dμ 或取 ∫_Y^μ (Y-t)/V(t) dt

這是否可當做是 μ 的對數概似函數？也就是說：可否假設

f*(y) = e^{q(y) + ∫_y^μ (y-t)/V(t) dt}

是 Y 的 p.d.f.？因為 q(y) 可以選擇使 f*(y) 在某個集合 A 之中積分為 1, 故 f* 視為一個 p.d.f. 沒有問題：

1 = ∫_A e^{q(y) + ∫_y^μ (y-t)/V(t) dt} dy

積分號下對 μ 微分得

0 = ∫_A [(y-μ)/V(μ)] f*(y) dy

0 = ∫_A {[(y-μ)/V(μ)]^2 - 1/V(μ) - (y-μ)V'(μ)/V(μ)^2} f*(y) dy

由第一式得 E*[Y] = μ, 將其代入第二式得 E*[(Y-μ)^2] = V(μ)，此處 E*[．] 表示在假設的 p.d.f. f* 之下計算期望值。這表示 f*(y; μ) 可當作 μ 的概似函數，l*(μ) 則是對數概似函數。由於它不是真正 Y 的機率分布導得的概似函數，稱之為「準概似函數 (quasi-likelihood function)」, 由 Var[Y] = V(E[Y]) 的形式建構準概似函數進行推論的方法，稱為準概似法。

今有資料 Yi, i = 1, ..., n，各對應解釋變數資料 x_i, 與 μ_i = E[Yi] 有如下關係：

η_i = g(μ_i) = x_i'β, i = 1, ..., n

若 Var[Yi] = V(μ_i) 函數形式已知，則得對數準概似函數

l*(β) = ∫_{Yi}^{μ_i} (Yi - t)/V(t) dt

故得求 β 之 qMLE 的概似方程式

Σ_i [(Yi - μ_i)/(V(μ_i) g'(μ_i))] x_i = 0

或每項前面加權量 Wi。這與前面正確指數族的 MLE 法相同，前面關於 β 之 MLE 算法也同樣適用。

若 Yi, i = 1, ..., n 為多變量，例如多項群體抽樣、重複量測資料，

E[Yi] = μ_i, η = g(μ) = B'x_i, Cov[Yi] = V(μ_i)

第 i 個體資料是一 k×1 行向量，其期望值也是，經一個連結函數轉成線性預測子 η，如果要求它和 μ 可互轉，則 g 是 R^k 到 R^k 的可逆變換，它與預測變數資料 x_i 之間的連繫不只是一個與 x_i 同階的向量，而是一個 p×k 矩陣 B。所有 Yi 來自典型指數族，則 B 的對數概似函數為

l(B; Y) = Σ _i [(θ_i．Yi) + c(θ_i) + q(Yi)]

先考慮 n = 1, 把註標 i 拿掉，

▽_θ l(B; Y) = Y + ▽_θ c(θ)
▽_μ l(B; Y) = ▽_μ(θ') (Y - μ)

第一式是概似函數（純量）對參數 θ 偏微，結果第二項即 - μ；第二式改為對 μ 偏微，則應用連鎖律，先逐一考慮 l 函數對 θ_j 的偏導，乘以 θ_j 對 μ 的偏導，用矩陣表示就是第一式的右邊再左乘 θ' 對 μ 的偏導構成的矩陣，由反函數定理， ▽_μ(θ') = [▽_θ(μ')]^{-1}，即 ▽_θ ▽_θ' C(θ) = V(μ)。類似，