广义线性模型（Generalized_Linear_Model）

1. 线性回归

1.1 多元线性回归模型

给定训练数据集
D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯   , ( x i , y i ) , … , ( x N , y N ) } \begin{aligned} \\& D = \left\{ \left( \mathbf{x}_{1}, y_{1} \right), \left( \mathbf{x}_{2}, y_{2} \right), \cdots, \left(\mathbf{x}_i,y_i\right),\dots, \left( \mathbf{x}_{N}, y_{N} \right) \right\} \end{aligned} ​D={(x1​,y1​),(x2​,y2​),⋯,(xi​,yi​),…,(xN​,yN​)}​
其中，${\mathbf{x}}_i\in \mathcal{X}\subseteq {\mathbb{R}}^n,y_i\in \mathcal{Y}\subseteq \mathbb{R}$。

多元线性回归模型：
$$f\left(\mathbf{x}\right)=\mathbf{w}\cdot \mathbf{x}+b=\sum _{i=1}^n{w}^{\left(i\right)}\cdot x^{\left(i\right)}+b$$
其中，$\mathbf{x}\in \mathcal{X}\subseteq {\mathbb{R}}^n$是输入记录，$\mathbf{w}={\left(w^{\left(1\right)},w^{\left(2\right)},\dots ,w^{\left(n\right)}\right)}^{\top }\in {\mathbb{R}}^n$和$b\in \mathbb{R}$是模型参数，$\mathbf{w}$称为权值向量，$b$称为偏置，$\mathbf{w}\cdot \mathbf{x}$为$\mathbf{w}$和$\mathbf{x}$的内积。

当$n=1$时，模型为一元线性回归模型：
$$f\left(x\right)=w\cdot x+b$$
其中，$w\in \mathbb{R}$和$b\in \mathbb{R}$为模型参数。

令
$$\begin{gathered} \hat{\mathbf{w}}={\left(\mathbf{w},b\right)}^{\top } \\ \hat{\mathbf{x}}={\left(\mathbf{x},1\right)}^{\top } \end{gathered}$$
则多元线性回归模型可简化为
$$f\left(\hat{\mathbf{x}}\right)=\hat{\mathbf{w}}\cdot \hat{\mathbf{x}}$$
其中，$\hat{\mathbf{x}}$为增广特征向量，$\hat{\mathbf{w}}$为增广权重。

1.2 多元线性回归参数学习——经验风险最小化与结构风险最小化

损失函数：平方损失损失函数
$$L\left(y,f\left(\mathbf{x}\right)\right)={\left(y-f\left(\mathbf{x}\right)\right)}^2$$

经验风险
$$\begin{array}{r}{R}_{emp}\left(f\right)=\frac{1}{N}\sum _{i=1}^{N}L\left(y_i,f\left({\mathbf{x}}_i\right)\right)\end{array}$$

模型参数最优解：
$$\begin{array}{rl}{\hat{\mathbf{w}}}^{\ast }& ={\mathrm{argmin}}_{\hat{\mathbf{w}}}\sum _{i=1}^N{\left(y_i-f\left({\hat{\mathbf{x}}}_i\right)\right)}^2\\ & ={\mathrm{argmin}}_{\hat{\mathbf{w}}}\sum _{i=1}^N{\left(y_i-\hat{\mathbf{w}}\cdot {\hat{\mathbf{x}}}_i\right)}^2\end{array}$$

基于均方误差最小化来进行模型求解的方法称为“最小二乘法”（least square method）。

等价的，模型参数最优解：
$${\hat{\mathbf{w}}}^{\ast }={\mathrm{argmin}}_{\hat{\mathbf{w}}}{\left(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}\right)}^{\top }\left(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}\right)$$
其中，
$$\begin{gathered} \mathbf{X}=\left(\begin{array}{cc}{\mathbf{x}}_1^{\top }& 1\\ {\mathbf{x}}_2^{\top }& 1\\ ⋮& ⋮\\ {\mathbf{x}}_N^{\top }& 1\end{array}\right)=\left(\begin{array}{c}{\hat{\mathbf{x}}}_1^{\top }\\ {\hat{\mathbf{x}}}_2^{\top }\\ ⋮\\ {\hat{\mathbf{x}}}_N^{\top }\end{array}\right) \\ \mathbf{y}={\left(y_1,y_2,\dots ,y_N\right)}^{\top } \end{gathered}$$

令$E_{\hat{\mathbf{w}}}={\left(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}\right)}^{\top }\left(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}\right)$，对$\hat{\mathbf{w}}$求偏导，得
$$\frac{\partial E_{\hat{\mathbf{w}}}}{\partial \hat{\mathbf{w}}}=2{\mathbf{X}}^{\top }\left(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y}\right)$$

当${\mathbf{X}}^{\top }\mathbf{X}$为满秩矩阵或正定矩阵时，令上式为零可得最优闭式解
$${\hat{\mathbf{w}}}^{\ast }={\left({\mathbf{X}}^{\top }\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$$

当上述条件不满足时，可使用主成分分析（PCA）等方法消除特征间的线性相关性，再使用最小二乘法求解。
或者通过梯度下降法，初始化${\hat{\mathbf{w}}}_0=0$，进行迭代
$$\hat{\mathbf{w}}\leftarrow \hat{\mathbf{w}}-\alpha {\mathbf{X}}^{\top }\left(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y}\right)$$
其中，$\alpha$是学习率。

结构风险：
$$\begin{array}{r}{R}_{str}=\frac{1}{N}\sum _{i=1}^{N}L\left(y_i,f\left({\mathbf{x}}_i\right)\right)+\lambda J\left(f\right)\end{array}$$

岭回归（Ridge Regression）：
$$\begin{array}{r}{R}_{str}=\frac{1}{N}\sum _{i=1}^{N}L\left(y_i,f\left({\mathbf{x}}_i\right)\right)+\alpha \Vert \mathbf{w}{\Vert }^2,\alpha \ge 0\end{array}$$

套索回归（Lasso Regression）：
$$\begin{array}{r}{R}_{str}=\frac{1}{N}\sum _{i=1}^{N}L\left(y_i,f\left({\mathbf{x}}_i\right)\right)+\alpha \Vert \mathbf{w}{\Vert }_1,\alpha \ge 0\end{array}$$

弹性网络回归（Elastic Net）：
$$\begin{array}{r}{R}_{str}=\frac{1}{N}\sum _{i=1}^{N}L\left(y_i,f\left({\mathbf{x}}_i\right)\right)+\alpha \rho \Vert \mathbf{w}{\Vert }_1+\frac{\alpha \left(1-\rho \right)}{2}\Vert \mathbf{w}{\Vert }^2,\alpha \ge 0,1\ge \rho \ge 0\end{array}$$

1.3 多元线性回归模型应用

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model, model_selection

def load_data():
    diabetes = datasets.load_diabetes()
    return model_selection.train_test_split(diabetes.data,diabetes.target,test_size=0.25,random_state=0) 

def test_LinearRegression(*data):
    X_train,X_test,y_train,y_test=data
    regr = linear_model.LinearRegression()
    regr.fit(X_train, y_train)
    print('Coefficients:%s, intercept %.2f'%(regr.coef_,regr.intercept_))
    print("Residual sum of squares: %.2f"% np.mean((regr.predict(X_test) - y_test) ** 2))
    print('Score: %.2f' % regr.score(X_test, y_test))

if __name__=='__main__':
    X_train,X_test,y_train,y_test=load_data() 
    test_LinearRegression(X_train,X_test,y_train,y_test) 

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model,model_selection

def load_data():
    diabetes = datasets.load_diabetes()
    return model_selection.train_test_split(diabetes.data,diabetes.target,
        test_size=0.25,random_state=0) 
def test_Lasso(*data):
    X_train,X_test,y_train,y_test=data
    regr = linear_model.Lasso()
    regr.fit(X_train, y_train)
    print('Coefficients:%s, intercept %.2f'%(regr.coef_,regr.intercept_))
    print("Residual sum of squares: %.2f"% np.mean((regr.predict(X_test) - y_test) ** 2))
    print('Score: %.2f' % regr.score(X_test, y_test))
    
def test_Lasso_alpha(*data):
    X_train,X_test,y_train,y_test=data
    alphas=[0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50,100,200,500,1000]
    scores=[]
    for i,alpha in enumerate(alphas):
        regr = linear_model.Lasso(alpha=alpha)
        regr.fit(X_train, y_train)
        scores.append(regr.score(X_test, y_test))
    
    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.plot(alphas,scores)
    ax.set_xlabel(r"$\alpha$")
    ax.set_ylabel(r"score")
    ax.set_xscale('log')
    ax.set_title("Lasso")
    plt.show()
    
if __name__=='__main__':
    X_train,X_test,y_train,y_test=load_data() 
    test_Lasso(X_train,X_test,y_train,y_test) 
    test_Lasso_alpha(X_train,X_test,y_train,y_test) 

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model,model_selection

def load_data():
    diabetes = datasets.load_diabetes()
    return model_selection.train_test_split(diabetes.data,diabetes.target,
        test_size=0.25,random_state=0) 

def test_Ridge(*data):
    X_train,X_test,y_train,y_test=data
    regr = linear_model.Ridge()
    regr.fit(X_train, y_train)
    print('Coefficients:%s, intercept %.2f'%(regr.coef_,regr.intercept_))
    print("Residual sum of squares: %.2f"% np.mean((regr.predict(X_test) - y_test) ** 2))
    print('Score: %.2f' % regr.score(X_test, y_test))
    
def test_Ridge_alpha(*data):
    X_train,X_test,y_train,y_test=data
    alphas=[0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50,100,200,500,1000]
    scores=[]
    for i,alpha in enumerate(alphas):
        regr = linear_model.Ridge(alpha=alpha)
        regr.fit(X_train, y_train)
        scores.append(regr.score(X_test, y_test))
    
    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.plot(alphas,scores)
    ax.set_xlabel(r"$\alpha$")
    ax.set_ylabel(r"score")
    ax.set_xscale('log')
    ax.set_title("Ridge")
    plt.show()
    
if __name__=='__main__':
    X_train,X_test,y_train,y_test=load_data() 
    test_Ridge(X_train,X_test,y_train,y_test) 
    test_Ridge_alpha(X_train,X_test,y_train,y_test) 

2. 逻辑斯谛回归

2.1 sigmoid函数与二分类逻辑斯谛回归模型

sigmoid函数：
$$sigmoid\left(z\right)=\sigma \left(z\right)=\frac{1}{1+e^{-z}}$$
其中，$z\in \mathbb{R}$，$sigoid\left(z\right)\in \left(0,1\right)$。

sigmoid函数的导数：
$${\sigma }^{\prime }\left(z\right)=\sigma \left(z\right)\left(1-\sigma \left(z\right)\right)$$

二分类逻辑斯谛回归模型是如下的条件概率分布：
$$\begin{array}{rl}P\left(y=1|\mathbf{x}\right)& =\sigma \left(\mathbf{w}\cdot \mathbf{x}+b\right)\\ & =\frac{1}{1+\exp \left(-\left(\mathbf{w}\cdot \mathbf{x}+b\right)\right)}\\ & =\frac{\exp \left(\mathbf{w}\cdot \mathbf{x}+b\right)}{1+\exp \left(\mathbf{w}\cdot \mathbf{x}+b\right)}\\ P\left(y=0|\mathbf{x}\right)& =1-\sigma \left(\mathbf{w}\cdot \mathbf{x}+b\right)\\ & =\frac{1}{1+\exp \left(\mathbf{w}\cdot \mathbf{x}+b\right)}\end{array}$$
其中，$\mathbf{x}\in {\mathbb{R}}^n$， y ∈ { 0 , 1 } y \in \left\{ 0, 1 \right\} y∈{0,1}，$\mathbf{w}\in {\mathbb{R}}^n$是权值向量，$b\in \mathbb{R}$是偏置，$\mathbf{w}\cdot \mathbf{x}$为向量内积。

可将权值权值向量和特征向量加以扩充，即增广权值向量$\hat{\mathbf{w}}={\left(w^{\left(1\right)},w^{\left(2\right)},\cdots ,w^{\left(n\right)},b\right)}^{\top }$，增广特征向量$\hat{\mathbf{x}}={\left(x^{\left(1\right)},x^{\left(2\right)},\cdots ,x^{\left(n\right)},1\right)}^{\top }$，则逻辑斯谛回归模型：
$$\begin{array}{rl}& \\ & P\left(y=1|\hat{\mathbf{x}}\right)=\frac{\exp \left(\hat{\mathbf{w}}\cdot \hat{\mathbf{x}}\right)}{1+\exp \left(\hat{\mathbf{w}}\cdot \hat{\mathbf{x}}\right)}\\ & P\left(y=0|\hat{\mathbf{x}}\right)=\frac{1}{1+\exp \left(\hat{\mathbf{w}}\cdot \hat{\mathbf{x}}\right)}\end{array}$$

2.2 二分类逻辑斯谛回归参数学习——最大似然估计

给定训练数据集
D = { ( x ^ 1 , y 1 ) , ( x ^ 2 , y 2 ) , ⋯   , ( x ^ N , y N ) } \begin{aligned} \\& D = \left\{ \left( \hat{\mathbf{x}}_{1}, y_{1} \right), \left( \hat{\mathbf{x}}_{2}, y_{2} \right), \cdots, \left( \hat{\mathbf{x}}_{N}, y_{N} \right) \right\} \end{aligned} ​D={(x^1​,y1​),(x^2​,y2​),⋯,(x^N​,yN​)}​
其中， x ^ i ∈ R n + 1 , y i ∈ { 0 , 1 } , i = 1 , 2 , ⋯   , N \hat{\mathbf{x}}_{i} \in \mathbb{R}^{n+1}, y_{i} \in \left\{ 0, 1 \right\}, i = 1, 2, \cdots, N x^i​∈Rn+1,yi​∈{0,1},i=1,2,⋯,N。

设
$$\begin{array}{rl}& \\ & P\left(y=1|\hat{\mathbf{x}}\right)=\sigma \left(\hat{\mathbf{w}}\cdot \hat{\mathbf{x}}\right),P\left(y=0|\hat{\mathbf{x}}\right)=1-\sigma \left(\hat{\mathbf{w}}\cdot \hat{\mathbf{x}}\right)\end{array}$$
似然函数
$$\begin{array}{rl}L\left(\hat{\mathbf{w}}\right)& =\prod _{i=1}^{N}P\left(y_i|{\hat{\mathbf{x}}}_i\right)\\ & =\prod _{i=1}^N{\left[\sigma \left(\hat{\mathbf{w}}\cdot {\hat{\mathbf{x}}}_i\right)\right]}^{y_i}{\left[1-\sigma \left(\hat{\mathbf{w}}\cdot {\hat{\mathbf{x}}}_i\right)\right]}^{1-y_i}\end{array}$$

因为似然函数累乘会可能出现下溢的情况，可以转换为对数似然函数（累加）
$$\begin{array}{rl}& \\ l\left(\hat{\mathbf{w}}\right)& =\log L\left(\hat{\mathbf{w}}\right)\\ & =\sum _{i=1}^N\left[y_i\log \sigma \left(\hat{\mathbf{w}}\cdot {\hat{\mathbf{x}}}_i\right)+\left(1-y_i\right)\log \left(1-\sigma \left(\hat{\mathbf{w}}\cdot {\hat{\mathbf{x}}}_i\right)\right)\right]\end{array}$$

最大似然估计
$${\hat{\mathbf{w}}}^{\ast }={\mathrm{argmax}}_{\hat{\mathbf{w}}}l\left(\hat{\mathbf{w}}\right)$$

最小负对数损失
$${\hat{\mathbf{w}}}^{\ast }={\mathrm{argmin}}_{\hat{\mathbf{w}}}-l\left(\hat{\mathbf{w}}\right)$$

令${\hat{y}}_i=\sigma \left(\hat{\mathbf{w}}\cdot {\hat{\mathbf{x}}}_i\right)$，则对数似然函数$l\left(\hat{\mathbf{w}}\right)$关于$\hat{\mathbf{w}}$的偏导数
$$\begin{array}{rl}\frac{\partial l\left(\hat{\mathbf{w}}\right)}{\partial \hat{\mathbf{w}}}& =-\sum _{i=1}^N\left(y_i\frac{{\hat{y}}_i\left(1-{\hat{y}}_i\right)}{{\hat{y}}_i}{\hat{\mathbf{x}}}_i-\left(1-y_i\right)\frac{{\hat{y}}_i\left(1-{\hat{y}}_i\right)}{1-{\hat{y}}_i}{\hat{\mathbf{x}}}_i\right)\\ & =-\sum _{i=1}^N\left(y_i\left(1-{\hat{y}}_i\right){\hat{\mathbf{x}}}_i-\left(1-y_i\right){\hat{y}}_i{\hat{\mathbf{x}}}_i\right)\\ & =-\sum _{i=1}^N{\hat{\mathbf{x}}}_i\left(y_i-{\hat{y}}_i\right)\end{array}$$

采用梯度下降法，初始化$\hat{\mathbf{w}}=0$，进行迭代
$${\hat{\mathbf{w}}}_{t+1}\leftarrow {\hat{\mathbf{w}}}_t+\alpha \sum _{i=1}^N{\hat{\mathbf{x}}}_i\left(y_i-{\hat{y}}_i^{{\hat{\mathbf{w}}}_t}\right)$$
其中，$\alpha$是学习率，${\hat{y}}_i^{{\hat{\mathbf{w}}}_t}$是当参数${\hat{\mathbf{w}}}_t$时模型的预测输出。

2.3 逻辑斯谛回归模型应用

from math import exp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split

def create_data():
    iris = load_iris()
    df = pd.DataFrame(iris.data, columns=iris.feature_names)
    df['label'] = iris.target
    df.columns = ['sepal length', 'sepal width', 'petal length', 'petal width', 'label']
    data = np.array(df.iloc[:100, [0,1,-1]])
    return data[:, :2], data[:, -1]

X, y = create_data()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

class LogisticRegressionClassifier:
    
    def __init__(self, max_iter=200, learning_rate=0.01):
        self.max_iter = max_iter
        self.learning_rate = learning_rate
        
    def sigmoid(self, x):
        return 1 / (1 + exp(-x))
    
    def data_matrix(self, X):
        data_mat = []
        for d in X:
            data_mat.append([1.0, *d])
        return data_mat
    
    def fit(self, X, y):
        data_mat = self.data_matrix(X)
        self.weights = np.zeros((len(data_mat[0]), 1), dtype=np.float32)
        
        for iter_ in range(self.max_iter):
            for i in range(len(X)):
                result = self.sigmoid(np.dot(data_mat[i], self.weights))
                error = y[i] - result
                self.weights += self.learning_rate * error * np.transpose([data_mat[i]])
        print('LogisticRegression Model(learning_rate={}, max_iter={})'.
              format(self.learning_rate, self.max_iter))
    
    def score(self, X_test, y_test):
        right = 0
        X_test = self.data_matrix(X_test)
        for x, y in zip(X_test, y_test):
            result = np.dot(x, self.weights)
            if (result > 0 and y == 1) or (result < 0 and y == 0):
                right += 1
        return right / len(X_test)

lr_clf = LogisticRegressionClassifier()
lr_clf.fit(X_train, y_train)

lr_clf.score(X_test, y_test)

x_points = np.arange(4, 8)
y_ = -(lr_clf.weights[1] * x_points + lr_clf.weights[0]) / lr_clf.weights[2]

plt.plot(x_points, y_)
plt.scatter(X[:50, 0], X[:50, 1], label='0')
plt.scatter(X[50:, 0], X[50:, 1], label='1')
plt.legend()

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model
from sklearn import model_selection

def load_data():
    iris=datasets.load_iris() 
    X_train=iris.data
    y_train=iris.target
    return model_selection.train_test_split(X_train, y_train,test_size=0.25,random_state=0,stratify=y_train)

def test_LogisticRegression(*data):
    X_train,X_test,y_train,y_test=data
    regr = linear_model.LogisticRegression()
    regr.fit(X_train, y_train)
    print('Coefficients:%s, intercept %s'%(regr.coef_,regr.intercept_))
    print('Score: %.2f' % regr.score(X_test, y_test))
    
def test_LogisticRegression_multinomial(*data):
    X_train,X_test,y_train,y_test=data
    regr = linear_model.LogisticRegression(multi_class='multinomial',solver='lbfgs')
    regr.fit(X_train, y_train)
    print('Coefficients:%s, intercept %s'%(regr.coef_,regr.intercept_))
    print('Score: %.2f' % regr.score(X_test, y_test))
    
def test_LogisticRegression_C(*data):
    X_train,X_test,y_train,y_test=data
    Cs=np.logspace(-2,4,num=100)
    scores=[]
    for C in Cs:
        regr = linear_model.LogisticRegression(C=C)
        regr.fit(X_train, y_train)
        scores.append(regr.score(X_test, y_test))
    
    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.plot(Cs,scores)
    ax.set_xlabel(r"C")
    ax.set_ylabel(r"score")
    ax.set_xscale('log')
    ax.set_title("LogisticRegression")
    plt.show()

if __name__=='__main__':
    X_train,X_test,y_train,y_test=load_data() 
    test_LogisticRegression(X_train,X_test,y_train,y_test) 
    test_LogisticRegression_multinomial(X_train,X_test,y_train,y_test) 
    test_LogisticRegression_C(X_train,X_test,y_train,y_test)