逻辑回归

逻辑回归原理

逻辑回归是一种对数概率回归，虽然名字中带有回归，实际上是一种解决分类问题的二分类算法。他是用sigmoid函数估计出样本属于正样本的概率。
sigmoid的函数表达式为：
$$h(z)=\frac{1}{1+e^{(-z)}}$$其中$z=w_0+w_1\cdot x_1+\dots +w_n\cdot x_n$，即$$h(z)=\frac{1}{1+e^{(-w^{T}x)}}$$$w$表示权重向量，$x$表示输入向量。

逻辑回归是一个线性模型

上面的线性映射实际上就是线性回归，最后通过sigmoid函数并不能改变他是线性模型的本质。样本属于正样本的概率为$$p(y=1|x)=h(x)$$样本属于负样本的概率为$$p(y=1|x)=1-h(x)$$其中$y$表示为类别标签，取值0或1，分别对应负样本和正样本。样本属于正样本和负样本的概率比的对数称为对数似然比：
$$ln\frac{p(y=1|x)}{p(y=0|x)}=ln\frac{\frac{1}{1+e^{(-w^{T}x)}}}{1-\frac{1}{1+e^{(-w^{T}x)}}}=w^{T}x$$分类归则：如果正样本的概率大于负样本的概率，即$$\frac{p(y=1|x)}{p(y=0|x)}>1$$也就是$$ln\frac{p(y=1|x)}{p(y=0|x)}=ln\frac{\frac{1}{1+e^{(-w^{T}x)}}}{1-\frac{1}{1+e^{(-w^{T}x)}}}=w^{T}x>0$$
即$w^{T}x>0$属于正样本，$w^{T}x<0$属于负样本，因此说逻辑回归是一个线性模型。

逻辑回归为什么使用交叉熵作为损失函数

假设训练样本集为$(x_i,y_i),i=1,2,3...l$，其中$x_i$为n维特征向量，$y_i$为类别标签，取值0或1，给定参数$w$和样本特征向量$x$，样本属于每个类的概率可以统一写成如下形式：$$p(y|x,w)=(h(x){)}^y(1-h(x){)}^{1-y}$$证明：令$y$为0或者1，上式分别属于正负样本的概率。逻辑回归输出的是样本属于一个类的概率，而样本的标签属于离散的0或者1，因此不适合直接使用欧式距离来作为损失函数，这里通过最大似然估计来确定参数。由于样本之间相互独立，训练样本集的似然函数为$$L(w)=\prod _{i=1}^{l}p(y|x_i,w)=\prod _{i=1}^l(h(x_i){)}^{y_i}(1-h(x_i){)}^{1-y_i}$$对数函数为：
$$lnL(w)=\sum _{i=1}^l{y}_{i}lnh(x_i)+（1-y_i）ln(1-h(x_i))$$这个函数称为二项式对数似然函数，要求该函数的最大值，等价于下面函数的最小值：
$$f(w)=-lnL(w)=-\sum _{i=1}^l{y}_{i}lnh(x_i)+（1-y_i）ln(1-h(x_i))$$根据求导公式，目标函数的梯度为：$${\nabla }_{w}f(w)=\sum _{i=1}^l（h(x_i)-y_i)x_i$$
目标函数的Hessian矩阵(对分量的二次偏导)为$${\nabla }_w^{2}f(w)=\sum _{i=1}^{l}h(x_i)(1-h(x_i))X_i$$如果单个样本的特征向量为$x_i=[x_{i1},x_{i2},...x_{in}{]}^T$，矩阵$X_i$定义为
$$X_i=\left[\begin{array}{ccc}{x}_{i1}^2& \cdots & x_{i1}{x}_{in}\\ ⋮& ⋮& ⋮\\ x_{in}{x}_{i1}& \cdots & x_{in}^2\end{array}\right]$$
此矩阵可以写成如下的乘积形式$$X_i=x_i{x}_i^T$$,对于任意的不为0的向量$x$有$$x_i{X}_i{x}_i^T\ge 0$$从而确定$X_i$半正定，另外由于$h(x_i)(1-h(x_i))>0$因此Hessian矩阵半正定，目标函数为凸函数。逻辑回归求解优化问题是不带约束的凸优化问题。如果使用欧式距离作为目标函数即$$L(w)=\frac{1}{2l}\sum _{i=1}^l(h(x_i)-y_i{)}^2$$则无法保证目标函数为凸函数。这就是使用交叉熵而不用欧式距离作为优化函数的主要原因之一。

逻辑回归梯度下降更新权重

使用梯度下降法求解梯度，迭代公式为$$w_{k+1}=w_k-\alpha (h_w(x_i)-y_i)x_i$$
$w$的初始值可全部设置为1。

代码实现

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import matplotlib as mpl
mpl.use('TkAgg')

def sigmoid(seta,X):
    # sigmoid 函数
    return 1. / (1 + np.exp(-(seta.dot(X.transpose()))))
#导数
def sigmoid_derivative(x):
    s = 1 / (1 + np.exp(-x))
    ds = s * (1 - s)
    return ds

def cross_entropy_error(y,t):
    # 交叉熵损失函数(平均损失)
    # t表示预测标签
    # y表示真实标签
    delta=1e-7  #添加一个微小值可以防止负无限大(np.log(0))的发生。
    return -np.sum(y*np.log(t+delta)+(1-y)*np.log(1-t+delta))/len(y)

def Optimizer(seta,X,Y):
    # 迭代参数，优化参数
    return seta+0.001*(Y.transpose()-sigmoid(seta,X)).dot(X)

def get_error_rate(T, Y):
    # 训练错误lv
    all_num = len(Y)
    error_num = 0
    for i in range(all_num):
        pred = T[0,i] > 0.5
        if pred != bool(Y[i]):
            error_num += 1
    return error_num * 1.0 / all_num

#画出分类结果
def plot_Decision(data,weights):
    plt.scatter(data[:,0], data[:,1], s=30, c=data[:,3], marker='s')
    x = np.arange(-0.0, 2.0, 0.1);  # 从-3.0 到 3.0 步长为0.1
    y = (-weights[0,2] - weights[0,0] * x) / weights[0,1]  # 决策边界 z=w0*x1+w1*x2+w2   z=0   x1代表横坐标 x2代表纵坐标 和点坐标对应 所以可以求出x2  即y
    plt.plot(x, y)
    plt.show()
    
if __name__ == '__main__':
    ones=np.random.random([1000,4])
    zeros=np.random.random([1000,4])+1
    ones[:,2]=1
    zeros[:,2]=1
    ones[:,3]=1#类别标签
    zeros[:,3]=0#类别标签

    train=np.append(ones,zeros,axis=0)#矩阵合并
    X_train, X_test, y_train, y_test = train_test_split(train[:,0:3], train[:,3], test_size=0.2)

    m, n = np.shape(X_train)
    seta=np.ones([1,n])
    for i in range(1000):
        T_train=sigmoid(seta, X_train)#训练集预测
        seta=Optimizer(seta, X_train, y_train)
        # 测试集合预测
        T_test = sigmoid(seta, X_test)
    print(seta)
    print("训练损失：", round(cross_entropy_error(y_train, T_train), 2), "训练准确率:", 1 - get_error_rate(T_train, y_train))
    print("测试损失：", round(cross_entropy_error(y_test, T_test), 2), "测试准确率:", 1 - get_error_rate(T_test, y_test))
    plot_Decision(train,seta)

结果显示

[[-6.41108114 -6.36417383 12.58099168]]
训练损失： 0.02 训练准确率: 1.0
测试损失： 0.03 测试准确率: 1.0