1. 廣義線性模型 1.1 從簡單線性回歸開始機器學習的任務是從已知的數據中提取知識, 從而完成對新數據的識別與預測, 即應用知識. 但是, 數據本身不會給予研究者想要的答案, 大多數時候人們很難通過觀察或者簡單的計算提取有用的結論. 為了從歷史數據中整合知識, 人們提出了"模型"這一概念, 認為 ...

1. 廣義線性模型

1.1 從簡單線性回歸開始

機器學習的任務是從已知的數據中提取知識, 從而完成對新數據的識別與預測, 即應用知識. 但是, 數據本身不會給予研究者想要的答案, 大多數時候人們很難通過觀察或者簡單的計算提取有用的結論. 為了從歷史數據中整合知識, 人們提出了"模型"這一概念, 認為數據是從模型中生成的, 遵循一定的規律, 如果規律是可知的, 那麼機器從數據中學習規律就是可能的.

由於數據不會說話, 最初的模型由人們假設(assumption)而來. 下麵用大家熟知的簡單線性回歸為例, 描述一個數據生成的過程. 假設數據為\(\{(y_i,x_i)\}_{i=1}^{n}\), 其中\(x_i\)是一個向量, 描述樣本點\(i\)的特征(features). 簡單線性回歸假設所有的樣本點由下麵的模型生成:

\[y=\theta^Tx+\epsilon \]

這裡\(\theta\)是一個與\(x\)同型的向量, 代表了模型的可知部分, \(\epsilon\)是一個不可知的正態零均值隨機誤差, 方差為\(\sigma^2\). 上述模型也可以寫成下麵的概率分佈形式:

\[y\sim N(\theta^Tx,\sigma^2). \]

我們獲得了從這個模型中生成的數據集\(\{(y_i,x_i)\}_{i=1}^{n}\), 希望從中學習到關於可知部分\(\theta\)的知識, 因為知道了\(\theta\), 就可以找到數據生成的過程, 從而對任何新觀測的特征\(x_{n+1}\), 即使不能知道對應的隨機誤差, 也能夠大概瞭解標簽\(y_{n+1}\)的期望\(\mathrm{E}[y_{n+1}]\), 這就是把從歷史數據中學到的知識應用於新數據的場景.

現在的目標是儘可能準確地估計一個\(\theta\)的值\(\hat{\theta}\)得到估計模型, 因為真實的\(\theta\)是不會直接呈現給我們的, 只能從歷史數據中估計. 對於每一個估計的\(\hat{\theta}\), 得到一個估計的概率分佈\(N(\hat{\theta}^Tx,\sigma^2)\), 即我們視為數據集是從這個分佈中產生的, 在期望意義下, \(h_{\hat\theta}(x)=\mathrm{E}[N(\hat{\theta}^Tx,\sigma^2)]=\hat{\theta}^Tx\)是從這個概率分佈中抽樣最應該出現的值.

如果模型準確, \(h_\hat{\theta}(x)\)和真實數據\(y\)的差異肯定儘可能小, 因此只需要想辦法刻畫模型結果和真實數據的差異並最小化之. 對兩個在\(\mathbb{R}^{n}\)上取值的向量, 衡量其差異的方法有很多, 選擇均方誤差作為其差異的度量, 要最小化兩個向量之間的均方誤差, 意味著我們的目標(objective)是

\[\hat\theta=\min_{\hat\theta}(h_\hat\theta(x)-y)^T(h_{\hat\theta}(x)-y). \]

優化不屬於我們討論的範疇, 我們只需要知道梯度下降法可以優化這個目標函數, 得到一個\(\hat\theta\)即可. 這樣, 我們就得到了估計模型.

但是, 簡單線性回歸模型顯得過於簡單了. 如果我們想獲得一個更加有應用價值的模型, 至少有下麵幾點是可以考慮的:

在數據生成機制上, 沒有任何的非線性運算元使得數據的生成機制被我們描述得過於理想化, 但過多的非線性也會增大估計難度. 因此可以考慮保留核心的線性生成部分\(\theta^Tx\), 但對結果施加一個非線性函數.
模型可能不由正態分佈所生成, 可以嘗試其他分佈.
我們也許感興趣的不是隨機數\(y\)本身, 而是其某種變換\(T(y)\).

考慮以上三點, 就從簡單線性回歸模型進入了廣義線性模型.

1.2 指數族

現在先考慮前兩點, 嘗試某種其他的分佈\(\mathrm{Distribution}(\eta)\), 其中核心未知參數是\(\theta^Tx\)的某個變換, 即

\[y\sim \mathrm{Distribution}(a(\theta^Tx)). \]

假定我們知道分佈的具體形式, 也獲得了數據集\((y,x)\), 那麼對於任何\(\hat{\theta}\), 如果數據來源於模型\(\mathrm{Distribution}(a(\hat\theta^Tx))\), 那麼此分佈的期望和真實數據\(y\)必定差異較小(這裡我們用"距離"來表示, 以區別於一般的作差), 也就是說我們的目標依然是

\[\min_{\hat{\theta}} \mathrm{Distance}(y,\mathrm{E}[\mathrm{Distribution}(a(\hat\theta^Tx))]). \]

在簡單線性回歸中, 距離函數被我們選擇為均方差, 現在我們且不論距離的選擇. 現在擺在面前的問題是如何計算分佈的期望, 如果有可能的話, 最簡單的期望形式必然是\(a(\theta^Tx)\)的某個映射, 這讓我們不能對所選分佈太過隨意. 有一類特殊的分佈滿足我們的需求: 指數族(exponential family). 由於指數族定義不依賴於模型, 先用\(\eta\)代替\(\theta^Tx\).

如果某個參數分佈族的密度函數滿足下麵的形式:

\[p(y;\eta)=b(y)\exp(\eta T(y)-a(\eta)), \]

就稱此分佈族是指數分佈族, 這裡\(a(\cdot)\), \(b(\cdot)\)和\(T(\cdot)\)都是可微函數. 我們關註它的期望:

\[\mathrm{E}[p(y;\eta)]=\int_{\mathbb{R}}yp(y;\eta)\mathrm{d}y=\int_{\mathbb{R}} yb(y)\exp(\eta T(y)-a(\eta))\mathrm{d}y, \]

直接求期望較困難, 註意到對分佈族中的任意概率分佈, 密度函數的積分為\(1\), 因此我們可以嘗試將指數族密度函數的積分對參數\(\eta\)求導, 其導數值應恆為\(0\). 指數族積分和求導可交換, 為計算提供了便利([2]中給出了不使用可交換性條件下的證明):

\[\begin{aligned} \frac{\partial }{\partial \eta}\int_{\mathbb{R}} p(y;\eta)\mathrm{d}y&=\int_{\mathbb{R}}\frac{\partial b(y)\exp(\eta T(y)-a(\eta))}{\partial \eta}\mathrm{d}y\\ &=\int_{\mathbb{R}}(T(y)-a'(\eta))b(y)\exp(\eta T(y)-a(\eta))\mathrm{d}y\\ &=T(\mathrm{E}[p(y;\eta)])-a'(\eta)\\ &=0, \end{aligned} \]

這樣我們求出了\(\mathrm{E}[p(y;\eta)]\)的一個變換的值, 如果\(T(\cdot)\)是可逆函數, 那麼\(\mathrm{E}[p(y;\eta)]=T^{-1}(a'(\eta))\); 特別地如果\(T(\cdot)\)是恆等變換, 那麼\(\mathrm{E}[p(y;\eta)]=a'(\eta)\). 這裡我們得到了指數族的期望, 證明指數族正是我們所關註的分佈族.

現在回到機器學習任務上, 假如我們選擇的分佈是指數族分佈\(\mathrm{ExponentialFamily}(\theta^Tx)\), 並取\(T(\cdot)\)為恆等映射, 那麼根據在簡單線性回歸一節中所討論的, 我們只需要:

對每一個估計的\({\theta}\), 代入指數族中, 得到期望觀測向量\(a'(\theta^Tx)\).
選擇一個合適的距離函數, 最小化\(\mathrm{Distance}(a'(\theta^Tx),y)\).
如果最小化由迭代完成, 就重覆上述兩個步驟; 否則用其他方式最小化.

1.3 廣義線性回歸實例

現在回顧簡單線性回歸模型, 由\(N(\theta^Tx,\sigma^2)\)也就是\(N(\eta,\sigma^2)\)生成, 我們考慮將正態分佈化為指數族的形式:

\[p(y;\eta)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{y^2-2\eta y+\eta^2}{2\sigma^2} \right)=\frac{\exp(-\frac{y^2}{2\sigma^2})}{\sqrt{2\pi\sigma^2}}\cdot\exp\left(\frac{\eta\cdot y-\eta^2}{\sigma^2} \right), \]

這裡存在另一個參數\(\sigma^2\), 這關係到指數族一個更廣泛的形式:

\[p(y;\eta,\tau)=b(y,\tau)\exp\left(\frac{\eta^TT(y)-a(\eta)}{c(\tau)} \right). \]

這裡\(\eta=\theta^Tx\)是向量, \(\theta\)是參數矩陣, \(T(y)\)是一個向量函數, 這種表述下文還會提到. 上式中\(\tau\)稱為擴散繫數, 也就是正態分佈中的\(\sigma^2\). 另外, 指數分佈族中包含眾多分佈: 正態分佈, 兩點分佈, 二項分佈, 泊松分佈, Gamma分佈等.

以兩點分佈\(B(1,p)\)為例, 它常常被視為二分類任務的生成模型. 要註意不能直接用\(\eta\)替換\(p\), 至少因為\(p\)的取值範圍只有\([0,1]\)而\(\eta=\theta^Tx\in\mathbb{R}\), 具體\(\eta\)應該如何表現還應當先將兩點分佈的概率函數(在這裡為分佈列)改寫為指數族的形式.

\[p(y;p)=p^{y}(1-p)^{1-y}=\exp\left(y\log p+(1-y)\log(1-p) \right)=\exp\left(y\cdot\log\frac{p}{1-p}+\log(1-p) \right). \]

這個密度形式與指數族的標準形式不同, 但只要施加變換

\[\eta=\log\frac{p}{1-p}, \]

同時有

\[p=\frac{e^{\eta}}{1+e^\eta}, \]

就可以把密度變為

\[p(y;p)=\tilde p(y;\eta)=\exp\left(y\cdot\eta-\log(1+e^{\eta}) \right), \]

變成指數族標準形式, 且\(a(\eta)=\log(1+e^{\eta})\). 對比前面給出的結論, 從\(B(1,p)\)中抽出的樣本\(y\)應當與\(a'(\eta)\)一致, 即最小化

\[\mathrm{Distance}(y,a'(\eta))=\mathrm{Distance}\left(y,\frac{e^{\theta^Tx}}{1+e^{\theta^Tx}} \right)=\mathrm{Distance}\left(y,\mathrm{sigmoid}(\theta^Tx) \right). \]

這就是二分類器中, 用線性函數sigmoid變換作概率估計的由來. 同時, 對分類問題, 儘管均方差用作距離函數依然是可行的, 但梯度下降中均方差表現不好, 因此我們往往使用交叉熵函數作距離度量(原因這裡不談):

\[\mathrm{CrossEntropy}(y,\hat y)=-\sum_{i=1}^{n}y_i\log_2 \hat{y}_i. \]

最後, 回到上文我們提到的\(\theta\)為參數矩陣的情況, 此時對應的分佈族是多項分佈\(B(1;p)\), 其中\(p=(p_1,\cdots,p_k)'\), 每次生成一個單位向量\(y=(y_1,\cdots,y_k)\)其中有且僅有一個分量\(y_i=1\)(以概率\(p_i\)), 其餘分量為\(0\), , 它常常被認為是多分類任務的生成模型.

嘗試寫出多項分佈的概率函數, 註意, 一個生成\(k\)維向量的多項分佈中, 未知參數只有\(k-1\)個, 這是因為\(\sum_ip_i=1\)且\(\sum_iy_i=1\).

\[\begin{aligned} p(y;p)&=p_1^{y_1}p_2^{y_2}\cdots p_k^{y_k} \\ &=\exp\left(\sum_{i=1}^{k-1}y_i\ln p_i+\left(1-\sum_{i=1}^{k-1}y_i \right)\ln\left(1-\sum_{i=1}^{k-1}p_i \right) \right)\\ &=\exp \left(\sum_{i=1}^{k-1}y_i\ln\frac{p_i}{1-\sum_{i=1}^{k-1}p_i}+\ln\left(1-\sum_{i=1}^{k-1}p_i \right) \right) \end{aligned} \]

至此, 令

\[\eta=\left(\ln\frac{p_1}{p_k},\cdots,\ln\frac{p_{k-1}}{p_k},0 \right), \]

可得\(p_i=p_ke^{\eta}(i=1,\cdots,k-1)\), 結合\(\sum_ip_i=1\)可得

\[p=\left(\frac{e^{\eta_1}}{1+\sum_{i=1}^{k-1}e^{\eta_i}},\cdots,\frac{e^{\eta_{k-1}}}{1+\sum_{i=1}^{k-1}e^{\eta_i}},\frac{1}{1+\sum_{i=1}^{k-1}e^{\eta_i}} \right), \]

上面兩個繁瑣的式子實際上就是\(\eta_i=\ln(p_{k-1}/p_k)\)以及\(p_i=e^{\eta_i}/\sum_je^{\eta_j}\). 此時密度函數為

\[p(y;p)=\tilde p(y;\eta)=\exp\left(y^T\eta-\ln \left(\frac{e^{\eta_k}}{\sum_{i=1}^{k}e^{\eta_i}} \right)\right):=\exp(y^T\eta-a(\eta)) \]

這裡\(a'(\eta)\)就是對其每個分量求導, 即

\[a'(\eta)=\begin{bmatrix} \frac{e^{\eta_1}}{\sum_i e^{\eta_i}} \\ \vdots \\ \frac{e^{\eta_k}}{\sum_{i}e^{\eta_k}} \end{bmatrix}=\begin{bmatrix} \frac{\exp(\theta_1^Tx)}{\sum_i\exp(\theta_i^Tx)}\\ \vdots \\ \frac{\exp(\theta_k^Tx)}{\sum_i\exp(\theta_i^Tx)} \end{bmatrix}. \]

此時\(\theta\)是參數矩陣, \(\theta_i\)代表矩陣的第\(i\)列.

1.4 廣義線性模型代碼

statsmodels庫中提供了單參數指數族對應的廣義線性模型. 以下假設我們有一個從泊松分佈中生成的模型, 這裡先驗證泊松分佈族是指數族.

\[p(y;\lambda)=\frac{\lambda^y}{y!}e^{-\lambda}=\frac{1}{y!}\exp(y\log \lambda-\lambda)\xlongequal{\eta=\log\lambda}\frac{1}{y!}\exp(y\eta-e^{\eta}), \]

也即數據從\(P(e^{\theta^Tx})\)中生成, 生成數據的期望為\(a'(\eta)=e^{\eta}=e^{\theta^Tx}\). 現在我們先生成數據. 假設\(x\)是三維數據, \(\theta=(1, 2, 5)^T\).

import numpy as np
import statsmodels.api as sm
np.random.seed(2000)

## Generate data from poisson distribution

theta = np.array([1, 2, 5])
x = np.random.normal(size=(100, 3))
y = np.random.poisson(lam=np.exp(np.dot(x, theta))) + np.random.normal(scale=1, size=(100,))

註意, GLM()函數是不會自己添加截距項的, 這裡簡化處理, 我們也認為模型不含截距項. 因為是泊松分佈族生成的數據, 因此要選擇分佈族為Poisson().

glm = sm.GLM(y, x, family = sm.families.Poisson())
res = glm.fit()
print(res.summary())

得到的結果為\(\hat{\theta}=(0.9964,1.9946,5.0010)'\), 非常接近真實值.

在GLM()函數中可以設定聯結函數(link function), 每一個分佈族會自帶預設的聯結函數, 這個函數就是連接分佈參數\(\lambda\)和自然參數\(\eta=\theta^Tx\)的函數, 本例中為\(\log\), 因為\(\eta=\log\lambda\). 有關聯結函數的更多信息, 可參考[3].