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].

Reference

[1] https://cs229.stanford.edu/notes2021fall/cs229-notes1.pdf

[2] https://xg1990.com/blog/archives/304

[3] https://www.statsmodels.org/stable/glm.html

轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/546431.html

標籤：其他

上一篇：Python gdal讀取MODIS遙感影像并結合質量控制QC波段掩膜資料

下一篇：快 40 歲，剛被裁。。