音乐播放器
My Brain
 
文章 标签
8

Powered by Gridea | Theme: Fog
载入天数...
载入时分秒...

机器学习实战(七):逻辑回归

逻辑回归

一、基本思想

逻辑回归是一种用于分类任务的经典机器学习算法。

我们之前介绍过机器学习算法大致可以分为两类:

P(x,y)=P(xy)P(y)P(yx)\begin{aligned} &①生成学习:估计P(x,y)=P(x|y)P(y)\\ &②判别学习:估计P(y|x) \end{aligned}

上节学习的朴素贝叶斯属于生成学习算法,而逻辑回归则属于判别学习算法,它与高斯朴素贝叶斯相对应,即判别式的高斯朴素贝叶斯。

逻辑回归的核心函数为:sigmoidsigmoid函数,也成为激活函数:

sigmoid(x)=σ(x)=11+exsigmoid(x)=\sigma(x)=\frac{1}{1+e^{-x}}

我们对逻辑回归模型的建模如下:

P(yx)=11+ey(wTx+b)P(y|x)=\frac{1}{1+e^{-y(w^Tx+b)}}

与感知机原理相同,我们利用升维将偏差项bb处理到ww中,则模型简化为:

P(yx)=11+ey(wTx)P(y|x)=\frac{1}{1+e^{-y(w^Tx)}}

二、逻辑回归的参数

在逻辑回归模型中,重要的参数也为ww,因此我们需要用到之前学过的几个概率估计方法对其进行估计。

2.1 极大似然估计(MLE)

利用极大似然估计逻辑回归的参数:在MLEMLE中,要极大化的条件数据为P(YX,w)P(Y|X,w)X,YX,Y为训练数据,我们要找到一个合适的ww使得特征向量集为XX时,观察到标签YY的概率最大

Xd×nX=[x1,x2,...,xn]Rd×nYnY=[y1,y2,...,yn]\begin{aligned} &X:d\times n维矩阵,即X=[\overrightarrow {x_1},\overrightarrow{x_2},...,\overrightarrow{x_n}]\in R^{d\times n}\\ &Y:n维向量,即Y=[y_1,y_2,...,y_n] \end{aligned}

MLEMLE的假设为:

P(YX,w)=i=1nP(yixi,w)P(Y|X,w)=\prod_{i=1}^nP(y_i|x_i,w)

由此,我们对参数ww进行的极大似然估计过程为:

w^MLE=argmaxw P(YX,w)=argmaxw i=1nP(yixi,w)           =argmaxw i=1nlogP(yixi,w)=argmaxw i=1nlog11+eyi(wTxi)           =argmaxw log(1+eyi(wTxi))=argminw log(1+eyi(wTxi))\begin{aligned} &\hat{w}_{MLE}=\underset{w}{argmax}~P(Y|X,w)=\underset{w}{argmax}~\prod_{i=1}^nP(y_i|x_i,w)\\ &~~~~~~~~~~~=\underset{w}{argmax}~\sum_{i=1}^n\log P(y_i|x_i,w)=\underset{w}{argmax}~\sum_{i=1}^n\log \frac{1}{1+e^{-y_i(w^Tx_i)}}\\ &~~~~~~~~~~~=\underset{w}{argmax}~-\log(1+e^{-y_i(w^Tx_i)})=\underset{w}{argmin}~\log(1+e^{-y_i(w^Tx_i)}) \end{aligned}

为了求解ww,我们引入函数l(w)l(w):在负数域上求解l(w)-l(w)的最大值:

l(w)=i=1nlog(1+eyi(wTxi))l(w)=\sum_{i=1}^n\log(1+e^{-y_i(w^Tx_i)})

我们要寻找max{l(w)}\max\{|-l(w)|\}所对应的ww

2.1.1 1-D Example

我们考虑一个1维问题,如下图所示:++表示正类,oo表示负类:

如上图所示,x0x_0恒为1,数据点特征的不同仅为x1x_1的不同,本质上是一个1维的问题:

x=[x0,x1],w=[w0,w1],l(w)=l([w0,w1])\begin{aligned} &x=[x_0,x_1],w=[w_0,w_1],l(w)=l([w_0,w_1]) \end{aligned}

右图为l(w)l(w)所对应的曲面,我们可以确定:w0=1,w1=0.7w_0=1,w_1=0.7时,l(w)|-l(w)|取得最大值,因此我们计算出w=[1,0.7]w=[1,0.7]

我们也可以用热力图更直观得得到ww的最佳取值:

以下是部分样例对上述图像的贡献图:

2.1.2 2-D Example

我们考虑一个2维问题,如下图所示:++表示正类,oo表示负类:

根据热力图,我们得到l(w)|-l(w)|w=[0.81.0.81]w=[-0.81.0.81]时取得最大值。

MLEMLE计算结果如下所示,其中红色表示正类别的概率很高。黑线表示MLEMLE学习的决策边界。

P(y=1x)=P(y=+1x)11+ewTx=11+ewTxwTx=wTxwTx=0w=[w1,w2]=[0.81,0.81],x=[x1,x2]wTx=w1x1+w2x2=0.81x1+0.81x2=0x1=x2\begin{aligned} &决策边界:P(y=-1|x)=P(y=+1|x)\\ &即:\frac{1}{1+e^{w^Tx}}=\frac{1}{1+e^{-w^Tx}}\rightarrow w^Tx=-w^Tx\rightarrow w^Tx=0\\ &w=[w_1,w_2]=[-0.81,0.81],x=[x_1,x_2]\\ &w^Tx=w_1\cdot x_1+w_2\cdot x_2=-0.81x_1+0.81x_2=0\rightarrow x_1=x_2 \end{aligned}

2.2 极大后验估计(MAP)

在极大后验估计中,我们将ww视为一个随机变量,并且可以预先指定它的先验分布。我们先预先指定ww的先验分布:

wN(0,σ2)P(w)=12πσ2ewTw2σ2w\sim \mathcal{N}(0,\sigma^2)\\ P(w)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{w^Tw}{2\sigma^2}}

这是逻辑回归的高斯近似。

以下计算过程中我们要用到链式法则:

P(A,BC)=P(AB,C)P(BC)P(A,B|C)=P(A|B,C)P(B|C)

同时需要注意XXww是无关的。

我们在MAPMAP中的目标是找到给定数据的最可能的模型参数,即使后验值最大化的参数:

w^MAP=argmaxw P(wD)=argmaxw P(Dw)P(w)P(D)          =argmaxw P(Dw)P(w)=argmaxw P(Y,Xw)P(w)          =argmaxw P(YX,w)P(Xw)P(w)=argmaxw P(YX,w)P(w)          =argmaxw logP(YX,w)+logP(w)          =argmaxw (i=1nlog(1+eyi(wTxi))+log12πσ212σ2wTw)          =argminw (i=1nlog(1+eyi(wTxi))+λwTw)\begin{aligned} &\hat{w}_{MAP}=\underset{w}{argmax}~P(w|D)=\underset{w}{argmax}~\frac{P(D|w)P(w)}{P(D)}\\ &~~~~~~~~~~=\underset{w}{argmax}~P(D|w)P(w)=\underset{w}{argmax}~P(Y,X|w)P(w)\\ &~~~~~~~~~~=\underset{w}{argmax}~P(Y|X,w)P(X|w)P(w)=\underset{w}{argmax}~P(Y|X,w)P(w)\\ &~~~~~~~~~~=\underset{w}{argmax}~\log P(Y|X,w)+\log P(w)\\ &~~~~~~~~~~=\underset{w}{argmax}~\big(-\sum_{i=1}^n\log (1+e^{-y_i(w^Tx_i)})+\log\frac{1}{\sqrt{2\pi\sigma^2}}-\frac1{2\sigma^2}w^Tw\big)\\ &~~~~~~~~~~=\underset{w}{argmin}~\big(\sum_{i=1}^n\log (1+e^{-y_i(w^Tx_i)})+\lambda w^Tw\big) \end{aligned}

其中:λ=12σ2\lambda=\frac1{2\sigma^2},与MLEMLE同理,我们引入函数l(w)l(w),在负数域上求解l(w)|-l(w)|最大值以得到ww

l(w)=i=1nlog(1+ey1(wTxi))+λwTwl(w)=\sum_{i=1}^n\log(1+e^{-y_1(w^Tx_i)})+\lambda w^Tw

三、更新参数

除了上述利用函数求极值的方法求解参数ww的方法,我们也可以用梯度下降的方式更新ww

3.1 损失函数

根据模型预测的概率为:

P(y=1x)=h(x)=σ(wTx)=11+e(wTx)        (1)P(y=0x)=1h(x)                                             (2)\begin{aligned} &P(y=1|x)=h(x)=\sigma(w^Tx)=\frac{1}{1+e^{-(w^Tx)}}~~~~~~~~(1)\\ &P(y=0|x)=1-h(x)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(2) \end{aligned}

合并(1)(1)(2)(2)式得:

P(yx,w)=h(x)y(1h(x))1y,y{0,1}\begin{aligned} &P(y|x,w)=h(x)^y(1-h(x))^{1-y},y\in\{0,1\} \end{aligned}

我们进行极大似然估计可得:

w^MLE=argmaxw P(YX,w)=argmaxw i=1nP(yixi,w)          =argmaxw i=1nh(xi)yi(1h(xi))1yi=argmaxw i=1n(yilnh(xi)+(1yi)ln(1h(xi)))\begin{aligned} &\hat{w}_{MLE}=\underset{w}{argmax}~P(Y|X,w)=\underset{w}{argmax}~\prod_{i=1}^nP(y_i|x_i,w)\\ &~~~~~~~~~~=\underset{w}{argmax}~\prod_{i=1}^nh(x_i)^{y_i}(1-h(x_i))^{1-y_i}=\underset{w}{argmax}~\sum_{i=1}^n\big(y_i\ln h(x_i)+(1-y_i)\ln(1-h(x_i))\big)\\ \end{aligned}

由此我们可以定义损失函数:

L(w)=i=1n(yilnh(xi)+(1yi)ln(1h(xi))L(w)=-\sum_{i=1}^n\big(y_i\ln h(x_i)+(1-y_i)\ln(1-h(x_i)\big)

我们希望损失函数最小化以得到最优的参数ww,与感知机的更新过程相同,我们求导以得到更新项。

3.2 梯度下降

梯度下降的过程为遍历数据集,利用每个样本对参数ww进行更新:令z=wTxz=w^Tx

每一个样本所对应的损失函数为:

l(w)=ylnh(x)(1y)ln(1h(x))l(w)=-y\ln h(x)-(1-y)\ln(1-h(x))

ww进行求导,求出更新项:

l(w)w=l(w)h(x)h(x)zzw,h(x)=11+e(wTx)l(w)h(x)=yh(x)+1y1h(x)=h(x)yh(x)(1h(x))h(x)z=ex(1+ex)=h(x)(1h(x))  zw=xl(w)w=h(x)yh(x)(1h(x))h(x)(1h(x))x=(h(x)y)x\begin{aligned} &\frac{\partial l(w)}{\partial w}=\frac{\partial l(w)}{\partial h(x)}\cdot \frac{\partial h(x)}{\partial z}\cdot \frac{\partial z}{\partial w},h(x)=\frac{1}{1+e^{-(w^Tx)}}\\ &①\frac{\partial l(w)}{\partial h(x)}=-\frac{y}{h(x)}+\frac{1-y}{1-h(x)}=\frac{h(x)-y}{h(x)(1-h(x))}\\ &②\frac{\partial h(x)}{\partial z}=\frac{e^{-x}}{(1+e^{-x})}=h(x)(1-h(x))~~③\frac{\partial z}{\partial w}=x\\ &因此:\frac{\partial l(w)}{\partial w}=\frac{h(x)-y}{h(x)(1-h(x))}\cdot h(x)(1-h(x))\cdot x=(h(x)-y)\cdot x \end{aligned}

则梯度下降的更新过程为:

wwηi=1n(h(xi)yi)x=wηi=1n(σ(wTxi)yi)xiw\rightarrow w-\eta\sum_{i=1}^n(h(x_i)-y_i)\cdot x=w-\eta\sum_{i=1}^n(\sigma(w^Tx_i)-y_i)\cdot x_i

四、模型实现

我们将通过手动实现与调库的方式去构造逻辑回归模型。

4.1 数据集

本次我们仍使用乳腺癌数据集,进行逻辑回归二分类器的实现,数据集详细内容可通过如下代码查看,不再过多赘述。

'''屏蔽warning'''
import warnings
warnings.filterwarnings("ignore")
'''导入重要的库'''
import numpy as np
import pandas as pd
from sklearn.datasets import load_breast_cancer

cancer=load_breast_cancer()
print(cancer.data)#输出数据集内容
print(cancer.target)#输出每个数据集的标签
print(cancer.feature_names)#输出数据集每一列对应的特征
print(cancer.target_names)#输出数据集的标签种类
'''将数据集存储为csv文件'''
df=pd.DataFrame(data=cancer.data,columns=cancer.feature_names)#构建二维表格
df['label']=cancer.target
df.to_csv("data/cancer_data.csv",index=False)

4.2 手动实现模型

我们利用梯度下降的方式对模型进行更新:

'''屏蔽warning'''
import warnings
warnings.filterwarnings("ignore")
'''导入重要的库'''
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
'''激活函数'''
def sigmoid(x):
    return 1/(1+np.exp(-x))
'''定义模型的类'''
class Logistic_Regression(object):
    '''初始化'''
    def __init__(self,eta,n_iter):
        self.eta=eta        #学习率
        self.n_iter=n_iter  #最大训练轮数
    '''模型训练'''
    def fit(self,X_train,y_train):
        feature_num=len(X_train[0]) #特征数
        w=np.zeros(feature_num)
        for i in range(self.n_iter):
            gradient=np.dot(sigmoid(np.dot(X_train,w))-y_train,X_train)
            w=w-self.eta*gradient
        self.w=w
    '''模型预测'''
    def predict(self,X_test):
        y_pred=[]
        for x in X_test:
            if(sigmoid(np.dot(self.w,x))>0.5):
                y_pred.append(1)
            else:
                y_pred.append(0)
        return np.array(y_pred)

'''读取数据'''
df=pd.read_csv("data/cancer_data.csv")
y=df['label']
df.drop("label",axis=1,inplace=True)
X=df.values
'''划分数据'''
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.2)#选取20%的数据作为测试集
'''初始化模型'''
log_reg=Logistic_Regression(0.001,10000)
log_reg.fit(X_train,y_train)
y_pred=log_reg.predict(X_test)
print(accuracy_score(y_test,y_pred))    #0.8947368421052632

4.3 调库实现模型

调库用到的函数为sklearn所提供的LogisticRegression函数,其函数原型如下:

LogisticRegression(penalty='l2', *, dual=False, tol=0.0001, C=1.0, fit_intercept=True, intercept_scaling=1, class_weight=None, random_state=None, solver='lbfgs', max_iter=100, multi_class='auto', verbose=0, warm_start=False, n_jobs=None, l1_ratio=None)

可以发现该模型并不是用梯度下降的方式去更新模型的,而是通过参数的估计,要用到优化求解器,其相关参数如下表所示:

模型参数 Parameter含义 备注
penalty 正则化项 广义线性模型的正则项,可选值包括L1正则项'l1'、L2正则项'l2'、复合正则'elasticnet'和无正则项None,默认值为'l2'。值得注意的是,正则项的选择应与优化求解器相匹配(详见solver参数)。
l1_ratio 正则权重系数 L1正则和L2正则的权重系数,取值空间为0-1。若为0,则相当于为L2正则;若为1,则为L1正则,否则为Elasticnet正则。
dual 对偶问题 默认值False,可设为True将问题转换为对偶问题(详见本博客SVM问题中原问题-对偶问题的推导),仅适用于采用L2正则化且求解器为'liblinear’的情况。当样本数少于特征数时,推荐为True。
tol 迭代阈值 求解器迭代求解时,停止迭代的目标函数改变阈值。
C 正则化系数倒数 注意,其值与正则化强度相反,即C值越小,正则化程度越大。其值必须为正,且默认值为1
fit_intercept 是否预设偏置 控制广义线性模型中是否预设偏置值
intercept_scaling 预设偏置值 广义线性模型中预设的偏置值,仅当求解器为'liblinear'同时fit_intercept时生效。注意:该偏置值会作为新的特征计算其系数,因此也会计入L1和L2正则。
class_weight 样本权重 用于处理样本不均衡问题。默认值为None,即各类别样本权重一样,可通过设置字典定义权重系数,或设为'balanced',即根据样本数自动计算权重。若在fit函数中设置sample_weight参数,两者作用会叠加。
random_state 随机状态 LR模型中的随机性主要体现在求解器迭代时对样本的随机选取,适用于当求解器为'liblinear'或'sag'时。
solver 求解器 sklearn中共提供了5种优化求解器,分别为'liblinear'、'sag'、'saga'、'newton-cg'和'lbfgs'。各求解器的适用条件不同,具体见后文。默认值为'liblinear'。
max_iter 最大迭代步数 ‘newton-cg’、'sag'和'lbfgs' 求解器所需要的最大迭代步数
multi_class 多分类策略 取值可为'ovr'、'multinomial'和'auto'。'ovr'即采用'one vs rest'策略对二分类模型进行集成;'multinomial'即采用'multinomial loss'直接求解多分类问题。默认为'auto',其会在两分类问题或求解器为'liblinear'时选择'ovr',而其它情况下选择'multinomial'。默认值为'auto'

在此我们不对参数作过多解释,具体调参过程就问题而论。

'''屏蔽warning'''
import warnings
warnings.filterwarnings("ignore")
'''导入重要的库'''
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.linear_model import LogisticRegression#导入逻辑回归模型

'''读取数据'''
df=pd.read_csv("data/cancer_data.csv")
y=df['label']
df.drop("label",axis=1,inplace=True)
X=df.values
'''划分数据'''
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.2)#选取20%的数据作为测试集
'''初始化模型'''
log_reg=LogisticRegression()
log_reg.fit(X_train,y_train)
y_pred=log_reg.predict(X_test)
print(accuracy_score(y_test,y_pred))    #0.9385964912280702