logistics回归原理推导+纯代码实现泰坦尼克号预测

前言

刚开始接触机器学习的时候，了解比较多的还是回归问题，但其实实际生活中，更多的是分类问题，例如人脸识别、目标识别、过滤垃圾邮件等等。而logistics可以说是最简单的分类算法。

sigmoid函数

sigmoid函数又称logistics函数，如图呈现s形，sigmoid可以将线性回归函数值映射到区间 ( 0,1 ) ，作为概率输出，大于或等于0.5判为正类，小于0.5判为父类。那么又有疑问了？那是不是logistics只能做二分类呢？那肯定不是啦，只要能做二分类，那肯定就能做多分类，因为我们可以在多分类中确定其中的一类，那么其他类就是负类了嘛，只是要做多几次。

logistics内部原理推导

首先我们用线性回归表达出来：
$$y={\theta }_0+{\theta }_{1}x+\dots +{\theta }_{n}x=X\theta$$
其中$X=\left[1,x_1,\dots ,x_n\right]$，$\theta =\left[{\theta }_1,\dots ,{\theta }_n\right]$。
利用logistics函数进行映射：
$$g\left(x\right)=g\left(X\theta \right)=\frac{1}{1+e^{-X\theta }}$$
那可以得到我们的目标函数：
$$\min J\left(\theta \right)=\frac{1}{2m}\sum _{i=1}^m\left(g\left(x^{\left(i\right)}\right)-y^{\left(i\right)}\right)=\frac{1}{2m}\sum _{i=1}^m\left(\frac{1}{1+e^{-x^{\left(i\right)}\theta }}-y^{\left(i\right)}\right)$$其中$m$为样本数，$i$代表第几个样本。
类似线性回归，我们希望能对目标函数求导，得到梯度函数，但是logistics模型的目标函数比线性回归复杂得很，所以求导是一个非常复杂的过程。于是转换表达方式，使用统计学上的最大似然法（MLE），当样本为正类时，其概率可以表达为：$p\left(y=1|x^{\left(i\right)};\theta \right)=g\left(x^{\left(i\right)}\right)$，当样本为反类时：$p\left(y=0|x^{\left(i\right)};\theta \right)=1-g\left(x^{\left(i\right)}\right)$。那么样本的概率可以表达为：
$p\left(y|x^{\left(i\right)};\theta \right)={\left(g\left(x^{\left(i\right)}\right)\right)}^{y^{\left(i\right)}}{\left(1-g\left(x^{\left(i\right)}\right)\right)}^{1-y\left(i\right)}$
构造似然函数：
$$l\left(\theta \right)=\prod _{i=1}^m{\left(g\left(x^{\left(i\right)}\right)\right)}^{y^{\left(i\right)}}{\left(1-g\left(x^{\left(i\right)}\right)\right)}^{1-y\left(i\right)}$$我们的目标就是使得$l(\theta )$尽可能大，因为值越大，对于是否为正类或负类更加确定。
对等式两边同时取对数：
$$\max \ln \left(l\left(\theta \right)\right)=\sum _{i=1}^m\left[y^{\left(i\right)}\ln g\left(x^{\left(i\right)}\right)+\left(1-y^{\left(i\right)}\right)\ln \left(1-g\left(x^{\left(i\right)}\right)\right)\right]$$那么接下来就使用梯度下降法进行优化：
$$\frac{\partial J\left(\theta \right)}{\partial \theta }=-\frac{1}{m}\sum _{i=1}^m\left[y^{\left(i\right)}\left(\frac{1}{g\left(x^{\left(i\right)}\right)}\right)-\left(1-y^{\left(i\right)}\right)\frac{1}{1-g\left(x^{\left(i\right)}\right)}\right]\frac{\partial g\left(x^{\left(i\right)}\right)}{\partial \theta }$$
结合：$g\left(z\right)=\frac{1}{1+e^{-z}}$,$g^{\prime }\left(z\right)=g\left(z\right)\cdot \left(1-g\left(z\right)\right)$，可以将其式子转化为：
$$\frac{\partial J\left(\theta \right)}{\partial \theta }=-\frac{1}{m}\sum _{i=1}^m\left[y^{\left(i\right)}\left(1-g\left(x^{\left(i\right)}\right)\right)-\left(1-y^{\left(i\right)}\right)g\left(x^{\left(i\right)}\right)\right]x^{{\left(i\right)}^T}=-\frac{1}{m}\sum _{i=1}^m\left[y^{\left(i\right)}-g\left(x^{\left(i\right)}\right)\right]x^{{\left(i\right)}^T}$$
接着我们使用梯度下降法对目标函数进行优化，算法步骤如下：
step1 $k=0$，$\epsilon =10^{-5}$，初始化${\theta }_k=[1,1,\dots ,1]$
step2 计算$g_k={\frac{\partial J\left(\theta \right)}{\partial \theta }}_{\theta ={\theta }_k}$，判断$\Vert g_k\Vert <\epsilon$，跳出循环，输出${\theta }_k$。
step3 取$d_k=-g_k$，armijo准则求$\alpha$。
step4 ${\theta }_k={\theta }_k+{\alpha }_k{d}_k$，$k=k+1$，转step 2。

纯代码实现

使用logistics回归预测泰坦尼克号数据集，预处理部分就不展式了，我们重点关注实现logistics的代码：

m = len(X)  #训练样本数量
b = np.ones(m)
X = np.insert(X,0,b,axis=1) #初始化
n = len(X[0])

def fun(theta):
    return -1/m*np.sum([y[i]*np.log(1/(1+np.e**(-X[i]@theta)))+(1-y[i])*np.log(1-1/(1+np.e**(-X[i]@theta))) for i in range(m)])

def gfun(theta):
    return -1/m*np.sum([(y[i]-1/(1+np.e**(-X[i]@theta)))*X[i].reshape(-1,1) for i in range(m)],axis=0)

#非精确线搜索
def Armijo(xk,dk):
    beta=0.5;sigma=0.2;m=0;maxm=20;mk=0;
    while m<maxm:
        if fun(xk+beta**m*dk)<=fun(xk)+sigma*beta**m*gfun(xk).T@dk:
            mk=m
            break
        m +=1
    newxk=xk+beta**mk*dk
    return newxk,mk

theta = np.ones(n).reshape(-1,1)
epsilon = 10**(-5)
maxn = 0

while maxn<10000:
    gk = gfun(theta)
    dk = -gk
    [theta,mk] = Armijo(theta,dk)
    if np.linalg.norm(gk)<epsilon:
        break
    maxn += 1
    print('第%s次迭代，梯度范数：'%maxn,np.linalg.norm(gfun(theta)),'极小点是：',theta,'函数值：',fun(theta))

我们使用精准率、召回率和F1值作为模型评价指标

#模型评价
TP=0; FP=0; FN=0; TN=0
for i in range(m):
    if X[i]@theta <0:
        if y[i] == 0:
            TP += 1
        else:
            FP +=1
    else:
        if y[i] == 0:
            FN += 1
        else:
            TN +=1    
P = TP/(TP+FP)  
print('模型对标签0的精准率为：',P)
R = TP/(TP+FN)
print('模型对标签0的召回率为：',R)
F1Score = 2*P*R/(P+R)
print('模型对标签0的F1-Score为：',F1Score)
#标签为1的
TP=0; FP=0; FN=0; TN=0
for i in range(m):
    if X[i]@theta >0:
        if y[i] == 1:
            TP += 1
        else:
            FP +=1
    else:
        if y[i] == 1:
            FN += 1
        else:
            TN +=1    
P = TP/(TP+FP)  
print('模型对标签1的精准率为：',P)
R = TP/(TP+FN)
print('模型对标签1的召回率为：',R)
F1Score = 2*P*R/(P+R)
print('模型对标签1的F1-Score为：',F1Score)

在迭代一万次后，可以得到在负类上F1值为0.83，在正类上达到0.72。总的来说，效果还是可以的，但是梯度下降法迭代的速度和效率都太慢了。