逻辑回归实现糖尿病预测（GD、SGD、K折）

本实验使用逻辑回归，根据数据集中包含的某些诊断测量值来预测对象是否患有糖尿病。将数据预处理后，先使用两种方法实现逻辑回归：梯度下降和调用sklearn库中的函数LogisticRegression。作为延伸，还使用了随机梯度下降和K折交叉验证。

基础理论

考虑二分类问题，给定数据集：
$$D=(x_1,y_1),(x_2,y_2),\cdots ,(x_N,y_N),x_i\subseteq R^n,y_i\in 0,1,i=1,2,\cdots ,N$$在线性可分情况下，决策边界可以表示为: $w^{T}x+b=0$，对于样本点若有 $w^{T}x+b>0$，那么可以判定其类别为 1，否则类别为0。

逻辑回归表达式

$$\begin{gathered} y=\mathrm{f}(\mathbf{z})=\frac{1}{1+{\mathrm{e}}^{-\mathbf{z}}}, \\ \mathbf{z}={\mathbf{w}}^{\mathrm{T}}\mathbf{x}+\mathrm{b} \end{gathered}$$其中 $f(\cdot )$ 为sigmoid函数。

使用了sigmoid函数映射后，分界面在 $z=0$ 处的间隔将变大，有利于优化$w$。

同时，模型计算的值 $y$ 可以理解为将样本划分为1类的概率，若其大于0.5，视为划分为1类，否则划分为0类。
$$\begin{gathered} P(y=1|x,w)=\frac{1}{1+e^{-w^{T}x+b}} \\ P(y=0|x,w)=\frac{e^{-w^{T}x+b}}{1+e^{-w^{T}x+b}} \end{gathered}$$

模型优化过程

令$\mathbf{p}=P(y=1|x,w)$， 则$P(y=0|x,w)=1-\mathbf{p}$。给定样本 $x_i$，其分类为$y_i$的概率为：
$$P\left({\mathrm{y}}_{\mathrm{i}}|{\mathbf{x}}_{\mathbf{i}}\right)={\mathrm{p}}^{{\mathrm{y}}_{\mathrm{i}}}{\left(1-\mathrm{p}\right)}^{1-{\mathrm{y}}_{\mathrm{i}}}$$优化模型参数$w$满足所有样本的情况，则要：
$$\max \mathrm{P}=\prod _{\mathrm{i}=1}^{\mathrm{n}}P\left({\mathrm{y}}_{\mathrm{i}}|{\mathrm{x}}_{\mathrm{i}}\right)=\prod _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}^{{\mathrm{y}}_{\mathrm{i}}}{\left(1-\mathrm{p}\right)}^{\left(1-{\mathrm{y}}_{\mathrm{i}}\right)}$$两边取对数后：
$$h(w)=\ln \mathrm{P}=\sum _{i=1}^n(y_i\ln p+(1-y_i)\ln (1-p))$$使用梯度下降算法优化参数。

计算梯度：
$$\begin{array}{rl}\nabla h\left(\mathbf{w}\right)& =\nabla (\sum _{i=1}^n({\mathbf{y}}_{\mathbf{i}}\ln \mathbf{p}+(1-{\mathbf{y}}_{\mathbf{i}})\ln (1-\mathbf{p})))\\ & =\sum _{i=1}^n({\mathbf{y}}_{\mathbf{i}}\mathbf{p}+\left(1-{\mathbf{y}}_{\mathbf{i}}\right)\frac{-{\mathbf{p}}^{\prime }}{1-\mathbf{p}})\end{array}$$因为sigmoid函数$f(\cdot )$的梯度为 $f(\cdot )\ast (1-f(\cdot ))$，故$p=\frac{1}{1+e^{-w^{T}x+b}}$的梯度为 ${\mathbf{p}}^{\prime }=\mathbf{p}(1-\mathbf{p})\mathbf{x}$。

代入上式得：
$$\begin{array}{rl}\nabla \mathrm{h}\left(\mathbf{w}\right)& =\sum _{i=1}^n({\mathbf{y}}_{\mathbf{i}}(1-\mathbf{p}){\mathbf{x}}_{\mathbf{i}}+(1-{\mathbf{y}}_{\mathbf{i}})-\mathbf{p}{\mathbf{x}}_{\mathbf{i}})\\ & =\sum _{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{y}}_{\mathrm{i}}-\mathrm{p}){\mathbf{x}}_{\mathbf{i}}\end{array}$$其中 $\mathbf{p}=\frac{1}{1+e^{-w^{T}x+b}}$.

优化参数：
$${\mathbf{w}}_{\mathrm{t}+1}={\mathbf{w}}_{\mathbf{t}}-\alpha \nabla h\left({\mathbf{w}}_{\mathbf{t}}\right)$$

数据集说明

此数据集来自kaggle，意在根据数据集中包含的某些诊断测量值来预测对象是否患有糖尿病。数据集包含多个医学预测变量（‘Pregnancies’,'Glucose’等）和一个目标变量（‘Outcome’）。

变量：
0列为Pregnancies(怀孕次数);
1列为Glucose(口服葡萄糖耐量试验中2小时后的血浆葡萄糖浓度);
2列为BloodPressure(舒张压,单位:mmHg);
3列为SkinThickness(三头肌皮褶厚度,单位:mm)
4列为Insulin(餐后血清胰岛素,单位:mm) 
5列为BMI,体重指数(体重(kg)/身高(m)^2)
6列为DiabetesPedigreeFunction(糖尿病家系作用)
7列为Age(年龄)
目标：
8列为Outcome(是否患有糖尿病,0或1)

库的导入

train_test_split函数用于将数据集分为训练集和测试集。StandardScaler函数用于对数据进行标准化处理。print_coefs函数打印某种逻辑回归方法的超参数、偏置项与准确率。

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

def print_coefs(name,w,b,acc):
    print("\n"+name+":")
    print("变量参数 w =",w.reshape(-1))
    print("偏置项 b = %.8f"%b)
    print("准确率:",acc)

数据预处理

1. 处理缺省值

   由于数据集中的数据存在异常0值，需要对这些缺省值进行处理。对于缺省值，可以使用均值或者中位数进行填充。本实验使用中位数进行填充。

2. 划分训练集与测试集

   这里使用了sklearn库中的train_test_split函数将数据集分为训练集和测试集，并设定训练集占总数据集的80%，测试集占20%。该函数的random_state参数用于设置随机数种子，以得到对数据集随机的划分。

3. 数据标准化

   使用sklearn库中的StandardScaler函数对数据进行标准化处理，使得数据集中的每个特征服从标准正态分布。这里，先用fit_transform方法处理训练集，再用transform方法处理测试集。这时，在归一化测试集时，使用的是训练集的统计量，使用训练集和测试集更相似，使算法在两者上的表现尽可能相同。

def Preprocess(data,seed=5932):
    medians=data.median()
    names = data.columns[1:-1]
    data[names] = data[names].replace(0, np.nan)    # 将除Pregnancies和Outcome两列外所有列中的缺省值标记为NaN
    data = data.fillna(medians)                     # 使用各列的中位数填充缺省值

    X,y=data.iloc[:,:-1],data.iloc[:,-1]            # 将数据划分为X和y
    # 使用train_test_split函数将数据划分为训练集和测试集,设定训练数据占80%
    X_train,X_test,y_train,y_test=train_test_split(X,y,train_size=0.8,random_state=seed)
    
    ss = StandardScaler()
    X_train = pd.DataFrame(ss.fit_transform(X_train),columns=X.columns)
    X_test = pd.DataFrame(ss.transform(X_test),columns=X.columns)   #分别对训练和测试数据的特征进行标准化处理
    return X_train,X_test,y_train,y_test

梯度下降(GD)实现逻辑回归

根据逻辑回归的原理，
$${\mathbf{w}}_{\mathrm{t}+1}={\mathbf{w}}_{\mathrm{t}}-\alpha \nabla h\left(\mathbf{w}\right)={\mathbf{w}}_{\mathrm{t}}-\alpha {\mathbf{X}}^{\mathrm{T}}\left(\frac{1}{1+\exp \left(-\mathbf{Xw}\right)}-\mathbf{y}\right)$$其中，$\mathbf{X}$和$\mathbf{w}$均经过增广。

于是，使用numpy中的矩阵处理，增广X_train和X_test，并将y_train和y_test写成列向量的形式。初始化$\mathbf{w}$为全1列向量，学习率$\alpha$为0.001，迭代次数为1000次，得到分类参数$\mathbf{w}$。将$\mathbf{w}$代入sigmoid函数，对测试集中的每个样本计算得到y的预测值。sigmoid(X_test@w)大于等于0.5时，预测为1，否则预测为0。将预测值与真实值进行比较，得到梯度下降法的准确率。

def sigmoid(x):
    return 1/(1+np.exp(-x))
    
def GradientDescent(X_train,X_test,y_train,y_test,alpha=0.001):
    X_train=np.hstack((X_train.values,np.ones((X_train.shape[0],1))))
    X_test=np.hstack((X_test.values,np.ones((X_test.shape[0],1))))          # 将X_train和X_test增加一列1
    y_train,y_test=y_train.values.reshape(-1,1),y_test.values.reshape(-1,1) # 将y_train和y_test转换为列向量
    w=np.ones([X_train.shape[1],1])                     # 初始化参数w
    for i in range(1000):
        w-=alpha*X_train.T@(sigmoid(X_train@w)-y_train) # 梯度下降
    y_predict=(sigmoid(X_test@w)>=0.5)                  # 预测结果二值化
    acc=np.mean(y_predict==y_test)                      # 计算准确率,即预测正确的样本数占总样本数的比例
    print_coefs("Gradient Descent",w[:-1],w[-1],acc)

调用库函数 LogisticRegression 实现逻辑回归

使用sklearn库中的LogisticRegression函数实现逻辑回归。其中，model.coef_对应于未增广的$\mathbf{w}$，model.intercept_对应于偏置项。

def Sklearn(X_train,X_test,y_train,y_test):
    from sklearn.linear_model import LogisticRegression
    model=LogisticRegression()          # 初始化逻辑回归模型
    model.fit(X_train,y_train)          # fit方法使用训练集进行训练,得到模型超参数
    y_predict=model.predict(X_test)     # predict方法使用测试集进行预测,得到预测的y
    acc=np.mean(y_predict==y_test)
    print_coefs("Sklearn",model.coef_,model.intercept_,acc)

随机梯度下降(SGD)实现逻辑回归

与梯度下降法相比，随机梯度下降法每次迭代只使用一个样本，因此速度更快。迭代公式为
$${\mathbf{w}}_{\mathrm{t}+1}={\mathbf{w}}_{\mathrm{t}}-\alpha {\mathbf{x}}_{\mathrm{t}}\left(\frac{1}{1+\exp \left(-{\mathbf{x}}_{\mathrm{t}}\mathbf{w}\right)}-y_{\mathrm{t}}\right)$$其中$\mathrm{t}$的取值为$1,2,...,n$，$n$为训练集的样本数，在这里为X_train.shape[0]。

随机梯度下降虽然迭代速度较快，但是可能会陷入局部最优解，准确率通常较梯度下降法低。

def SGD(X_train,X_test,y_train,y_test,alpha=0.001):
    X_train=np.hstack((X_train.values,np.ones((X_train.shape[0],1))))
    X_test=np.hstack((X_test.values,np.ones((X_test.shape[0],1))))
    y_train,y_test=y_train.values.reshape(-1,1),y_test.values.reshape(-1,1)
    w=np.ones([X_train.shape[1],1])
    for j in range(X_train.shape[0]):
        w-=alpha*(sigmoid(X_train[j]@w)-y_train[j])*X_train[j].reshape(-1,1)    # 随机梯度下降
    y_predict=(sigmoid(X_test@w)>=0.5)
    acc=np.mean(y_predict==y_test)
    print_coefs("Stochastic Gradient Descent",w[:-1],w[-1],acc)

K折划分实现逻辑回归

K折交叉验证将数据集划分为K个大小相似的互斥子集，每次使用K-1个子集的并集作为训练集，剩下的一个子集作为测试集，重复K次，可以更加充分地利用数据集。计算出多个模型的准确率后，从中找出一个泛化性能最好、超参数组合最佳的模型，使用这组超参数对整个训练集进行训练，得到最终的模型。

这里使用sklearn库中的KFold函数将数据集划分为K个子集，KFold函数的n_splits参数用于设置划分的份数（这里设为5），random_state参数用于设置随机数种子，以得到对数据集随机的划分。对每一种划分，使用LogisticRegression函数实现逻辑回归，得到每种模型的准确率，并取准确率最高的模型为最优模型，使用最优模型的超参数对完整训练集进行训练。

def KFoldSklearn(X_train,X_test,y_train,y_test,k=5,seed=5932):
    from sklearn.linear_model import LogisticRegression
    from sklearn.model_selection import KFold
    kf = KFold(n_splits=k,shuffle=True,random_state=seed)   # 初始化KFold对象,设定K=5
    max_acc=0                                               # max_acc用于记录几个模型中最高的准确率
    for train_index,test_index in kf.split(X_train):        # split方法用于产生每次划分的训练集的索引和测试集的索引
        model=LogisticRegression()
        model.fit(X_train.iloc[train_index],y_train.iloc[train_index])                  # iloc方法用于按索引进行切片
        acc=np.mean(model.predict(X_train.iloc[test_index])==y_train.iloc[test_index])  # 计算每个模型的准确率
        if max_acc<acc:
            max_acc=acc                                     # 找寻最高的准确率
            best_model=model                                # 记录最高准确率对应的模型,即最优模型
    y_predict=best_model.predict(X_test)                    # 使用最优模型的参数对整个测试集进行预测
    acc=np.mean(y_predict==y_test)                          # 在整个测试集上计算准确率
    print_coefs("KFold",best_model.coef_,best_model.intercept_,acc)

主函数

主函数中，使用库函数read_csv读取数据文件diabetes.csv，Preprocess函数将数据进行预处理与划分，得到训练集和测试集。GradientDescent、SGD、Sklearn和KFoldSklearn函数分别使用梯度下降、随机梯度下降、调用库函数LogisticRegression和K折划分实现逻辑回归，并打印出各自的超参数、偏置项和准确率。

if __name__=='__main__':
    data=pd.read_csv(r"diabetes.csv")
    datasets=Preprocess(data)
    GradientDescent(*datasets)      # 梯度下降
    SGD(*datasets)                  # 随机梯度下降
    Sklearn(*datasets)              # 调用库函数
    KFoldSklearn(*datasets)         # K折交叉验证

Gradient Descent:
变量参数 w = [ 0.32284836  1.06512242 -0.12022423 -0.06247633 -0.13979151  0.7101504  0.31527763  0.19169785]
偏置项 b = -0.72258697
准确率: 0.8636363636363636

Stochastic Gradient Descent:
变量参数 w = [ 0.94821771  1.00063203  0.93099904  0.93718738  0.96261441  0.96416734  0.98095596  0.93800946]
偏置项 b = 0.88192861
准确率: 0.6428571428571429

Sklearn:
变量参数 w = [ 0.31805227  1.0466346  -0.11284161 -0.05621546 -0.13238095  0.69474348  0.31040481  0.19211869]
偏置项 b = -0.71923251
准确率: 0.8636363636363636

KFold:
变量参数 w = [ 0.2517363   1.03985257 -0.07014447 -0.02436512 -0.1031414   0.65096068  0.27245426  0.28586422]
偏置项 b = -0.61968836
准确率: 0.8701298701298701

实验结果分析

可以看到，在当前随机数种子设定下，梯度下降和调用库函数进行逻辑回归的准确率达到86.36%。使用五折划分得到的准确率略有提升，为87.01%。相比之下，随机梯度下降得到的准确率较低，仅为64.28%。随机梯度下降每次迭代只使用一个样本，速度更快，但是准确率较低，更适合大体量数据集的训练。