【机器学习】1. 广义线性模型

1. 广义线性模型

指数族

若概率分布\(p(y;\eta)\)可以写成下面的形式（\(\eta\)不一定是分布本身指定的参数，可以经过变换）：

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

则称\(p(y;\eta)\)为指数型分布族，且

• 称\(\eta\)为自然参数或规范参数(canonical parameter)，不同的\(\eta\)对应不同的概率分布.
• \(T(y)\)是分布族的充分统计量. 在许多情况下\(T(\cdot)\equiv I(\cdot)\).
• \(a(\eta)\)是配分函数，用以使分布族的积分为\(1\).

对指数族\(Y\sim p(y;\eta)\)，若\(T\)为可逆函数，则\(\mathrm{E}(Y)=T^{-1}(a'(\eta))\). 直接证明不易，但已有证明表明指数族的积分和求导可交换，直接使用此结论，求\(\int_{\mathbb{R}}p(y;\eta)\mathrm{d}y\equiv 1\)关于\(\eta\)的偏导数，可以得到

\[\begin{aligned} &\quad \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} \]

若模型本身参数\(\theta\)与自然参数\(\eta\)不等，则link function \(g\)将模型参数联系到自然参数，有

\[g(\theta)=\eta. \]

如，二项分布\(B(1,p)\)是指数族，但其参数\(p\)不是自然参数，因为二项分布的密度为

\[p(y;p)=p^y(1-p)^{1-y}=\exp\left(y\ln \left(\frac{p}{1-p} \right)+\ln (1-p) \right), \]

因此取\(g(p)=\ln(\frac{p}{1-p}):=\eta\)，将二项分布改写为指数族的自然形式，就有

\[p(y)=\exp\left(y\eta+\ln\left(1-\frac{1}{1+e^{-\eta}}\right)\right). \]

广义线性模型

广义线性模型依赖于一个指数族\(\mathrm{ExponentialFamily}(\eta)\)，这个分布族决定了数据如何生成. 在一个广义线性模型中有以下三个要素：

1. 数据集为\(\{(x^{(i)},y^{(i)})\}_{i=1}^{n}\)，其中\(x^{(i)}\in\mathbb{R}^{d}\)，\(y^{(i)}\in\mathbb{R}^{K}\)，大部分情况下\(K=1\).
2. \(y|x;w \sim \mathrm{ExponentialFamily}(\eta)\)，这里的\(w\in\mathbb{R}^{d}\)是线性模型参数. 若\(K>1\)则\(W\in\mathbb{R}^{K\times d}\)为线性模型参数. 此时用\(W_1,\cdots,W_K\)表征其每一行.
3. \(\eta=w'x\)，即广义线性模型中线性部分拟合的是自然参数\(\eta\)，如果\(\eta\in \mathbb{R}^{K}\)则\(\eta=Wx\).
4. 模型的目标是估计\(\mathrm{E}[T(y)|x]\). 尽管许多情况下\(T(y)=y\)，但有时不同（如多项分布）.
5. 最后得到的模型将使用\(\eta\)来构造\(\mathrm{E}[T(y)|x]\).

以多项分布为基础的Softmax回归

下面以多项分布为例，构造以多项分布为基础的广义线性模型. 首先明确多项分布\(\mathrm{Multinomial}(\theta)\)是指数族，模型参数\(\theta\in\mathbb{R}^{K}\)，当\(y\sim \mathrm{Multinomial}(\theta)\)时假设从中抽取出的单个样本\(y^{(i)}\sim \{e_1,\cdots,e_K\}\)，这里\(e_i\)是除了第\(i\)个分量为\(1\)，其他分量均为\(0\)的标准正交基.

\[\begin{aligned} p(y;\theta)&=\theta_1^{y_1}\theta_2^{y_2}\cdots\theta_K^{y_K}\\ &=\exp\left(y_1\ln\theta_1+\cdots+y_{K-1}\ln\theta_{K-1}+\left(1-\sum_{k=1}^{K-1}y_k\right)\ln\theta_K \right)\\ &=\exp\left(y_1\ln\frac{\theta_1}{\theta_K}+\cdots+y_{K-1}\ln\frac{\theta_{K-1}}{\theta_K}+\ln\theta_K \right), \end{aligned} \]

这就说明\(\mathrm{Multinomial}(\theta)\)是一个指数族，\(b(y)=1\)，\(a(\eta)=\ln(\theta_K)\)，

\[T(y)=\begin{pmatrix} y_1 \\ \vdots \\ y_{K-1} \\ y_K \end{pmatrix},\quad \eta=\begin{pmatrix} \ln\frac{\theta_1}{\theta_K} \\ \vdots \\ \ln\frac{\theta_{K-1}}{\theta_K} \\ 0 \end{pmatrix}. \]

这是一个\(T(y)\ne y\)的例子，且为了形式美观令\(T(y)\in\mathbb{R}^{K}\)，尽管\(\mathbb{R}^{K-1}\)已经足够（参数\(\theta\)实际上只有\(K-1\)个自由变量）. 这样就得到了数据的生成模型，接下来使用极大似然法估计线性模型参数\(W\in\mathbb{R}^{}\). 在此之前，还要先知道自然参数\(\eta\)和分布参数\(\theta\)之间如何转化，这里\(\theta_K=\frac{1}{\sum_{k=1}^{K}e^{\eta_k}}\)且\(\theta_k=\theta_Ke^{\eta_k}\)，所以

\[\theta=\begin{pmatrix} \frac{e^{\eta_1}}{ \sum_{k=1}^{K}{e^{\eta_k}}} \\ \vdots \\ \frac{e^{\eta_{K}}}{\sum_{k=1}^{K}e^{\eta_k}} \end{pmatrix}=\mathrm{softmax}(\eta). \]

接下来对任意样本点\((x,y)\)，有\(y^{(i)}\sim \mathrm{Multinomial}(W x^{(i)})\)，故似然函数为

\[L(W)=\prod_{k=1}^{K}\theta_k^{y_k}=\prod_{k=1}^{K}\left(\frac{e^{\eta_k}}{\sum_{j=1}^{K}e^{\eta_j}}\right)^{y_k}=\prod_{k=1}^{K}\left(\frac{\exp(W_kx)}{\sum_{j=1}^{K}\exp(W_jx)} \right)^{y_k}, \]

接下来最大化对数似然，即最小化\(\ell(W)=-\ln L(W)\)，有

\[\ell(W)=-\sum_{k=1}^{K}y_k\left(W_kx-\ln\left(\sum_{j=1}^{K}\exp(W_jx) \right) \right). \]

实际训练模型时常采用梯度下降法，这里\(\ell(W)\)对各行的更新与样本所属的类有关，假设\(y=e_C\)即样本属于第\(C\)类，则

\[\begin{gather} \frac{\partial \ell(W)}{\partial W_C}=\frac{\exp(W_Cx)}{\sum_{k=1}^{K}\exp(W_jx)}\cdot x-x=\left(\frac{\exp(W_Cx)}{\sum_{k=1}^{K}\exp(W_jx)}-1 \right)x, \\ \frac{\partial \ell(W)}{\partial W_k}=\frac{\exp(W_kx)}{\sum_{j=1}^{K}\exp(W_jx)}x,\quad k\ne C, \end{gather} \]

容易发现，Softmax回归的优化仅对样本所属的类特殊处理，而对其他非样本类一视同仁. 上式还可以合并成

\[\frac{\partial \ell(W)}{\partial W}=(\mathrm{softmax}(W_Cx)-e_C)\cdot x'. \]

Softmax回归的简易代码实现

下面给出了Softmax回归的一个简易实现，这里训练不是批次的而是逐样本的. 模型对一个人造数据集和真实数据集iris都有较为不错的表现.

### Softmax Regression

import numpy as np
from sklearn.datasets import load_iris

class SoftmaxRegression:
    def __init__(self, X, y, num_classes):
        # X (n, d); y (n,)
        self.X = np.array(X)
        self.label = np.array(y)
        self.num_classes = num_classes
        self.n_samples, self.dim = self.X.shape
        self.params = np.random.normal(loc=0, scale=np.sqrt(2/(num_classes + self.dim)), size=(self.num_classes, self.dim))

    def normalize_Y(self):
        self.Y = np.zeros(shape=(self.n_samples, self.num_classes))
        self.dict = []
        for i, c in enumerate(self.label):
            if c not in self.dict:
                self.dict.append(c)
            ind = self.dict.index(c)
            self.Y[i, ind] = 1

    def _softmax(self, X):
        if len(X.shape) == 1:
            X = X - X.max()
        elif len(X.shape) == 2:
            X = X - X.max(axis=0, keepdims=True)
        return np.exp(X) / sum(np.exp(X))
    
    def predict_proba(self, samples):
        samples = np.array(samples)
        if len(samples.shape) == 1:
            return self._softmax(np.dot(self.params, samples))
        elif len(samples.shape) == 2:
            return self._softmax(np.dot(self.params, samples.T)).T
    
    def predict(self, samples):
        proba = self.predict_proba(samples)
        if len(proba.shape) == 1:
            return np.argmax(proba)
        elif len(proba.shape) == 2:
            return np.argmax(proba, axis=1)

    def fit(self, lr=0.05, max_iter=1e5):
        ## train for single sample
        self.normalize_Y()
        itter_time = 0
        ind = 0
        while itter_time < max_iter:
            y_pred = self.predict_proba(self.X[ind])
            grad = np.outer(y_pred - self.Y[ind], self.X[ind])
            self.params -= lr * grad
            ind = (ind + 1) % self.n_samples
            itter_time += 1
            if True in np.isnan(self.params):
                break
            if np.linalg.norm(grad) < 1e-10:
                break
    
def acc(pred, label):
    return (pred == label).mean()

## For synthetic datasets

np.random.seed(0)

k = 3; d = 5; n = 300
theta = np.random.uniform(0, 1, size=(k, d))
X_syn = [np.random.multivariate_normal(mean=np.random.normal(size=d), cov=np.identity(d), size=(100,)) for _ in range(3)]
X_syn = np.vstack(X_syn)
y_syn = np.repeat([0, 1, 2], 100)
clf = SoftmaxRegression(X_syn, y_syn, k)
clf.fit()
print(f'acc = {acc(clf.predict(X_syn), y_syn)}')

## For true datasets (iris)

df = load_iris()
iris_X, iris_y = df['data'], df['target']
clf_iris = SoftmaxRegression(iris_X, iris_y, 3)
clf_iris.fit()
print(f'acc of iris = {acc(clf_iris.predict(iris_X), iris_y)}')