Logistic回归详解

Logistic回归原理推导及代码实现

Logistic回归为概率型非线性回归模型，是研究二分类观察结果之间关系的一种多变量分析方法，在工业界应用广泛，因此了解其原理及实现较为重要。本篇文章将通过极大似然法对logistic回归进行推导，并利用python进行代码实现。

原理篇

首先考虑线性分类器z = w_0 + w_1x_1 +……+ w_kx_k，为了进行分类任务，利用sigmoid函数将z映射为概率进行分类。Logistic回归通过假设每个事件服从伯努利分布(1,p),而p则受sigmoid函数控制。根据伯努利分布，可写出每个事件的概率分布，再利用极大似然法可求出参数w_1…w_k，下面分别讨论每一个步骤

Sigmoid函数

Logistic回归中的sigmoid函数形式如下

$$\begin{cases} g(z) = \frac{1}{1+e^{-z}} \\ z = w_0 + w_1 x_1+......w_k x_k \end{cases}$$
其中
$$y = \begin{cases} 1 , &g(z)\geq 0.5 \\ 0 , &g(z)<0.5 \end{cases}$$
若利用概率将其写成一个等式来描述y的分布，即 $p(y) = g(z)^y(1-g(z))^{1-y}$
Sigmoid函数有一些特性，其导数如下

$\frac{\partial g(z)}{\partial z} = g(z) - g^2(z)$

极大似然推导

假设每个事件$y_i$ ~ $b(1,g(z_i))$，则根据上文$p(y_i) = g(z_i)^{y_i}(1-g(z_i))^{1-y_i}$，由此可求出似然函数L(w)

$L(w) = \prod_{i=1}^n g(z_i)^{y_i}(1-g(z_i))^{1-y_i}$
取对数求得对数似然函数l(w)
$l(w) = \sum_{i=1}^ny_iln g(z_i) + \sum_{i=1}^n (1-y_i)ln (1-g(z_i))$
对每个 $w_j$求其偏导可得
$\frac{\partial l(w)}{\partial w_j} = \sum_{i=1}^ny_i(1-g(z_i))x_{ij} - \sum_{i=1}^n(1-y_i)g(z_i)x_{ij}= \sum_{i=1}^n (y_i-g(z_i))x_{ij}$
可以发现，l(w)的梯度和线性回归类似，只是logistic回归是g(z)而线性回归是z。但由于这点不同导致logistic回归的似然函数写不出最大值的表达式，因此需要通过牛顿法或梯度上升来求得参数。

梯度上升法

根据上文有似然函数的梯度如下

$\frac{\partial l(w)}{\partial w_j} = \sum_{i=1}^n (y_i-g(z_i))x_{ij}$
写成矩阵形式即
$\nabla_{w_j} = \frac{\partial l(w)}{\partial w_j} = (y-\frac{1}{1+e^{-XW}})'x_j$
那么根据梯度上升每个参数的更新如下，其中 $\alpha$为梯度步长
$w_j := w_j + \alpha \nabla_{w_j}$
如此可求得Logistic回归的权重 $w_1......w_k$
让我们整理一下思路，Logistic的推导思路如下

1. 首先z与自变量x线性相关，为了对z进行分类，利用sigmoid函数将其映射为概率p
2. 假设每个事件服从伯努利分布，则伯努利分布的参数即为p
3. 由此可写出似然函数，利用极大似然估计，根据链式法则求导，再利用梯度上升求得权重$w_1......w_k$

实现篇

首先随便创建一个数据集，利用sklearn.datasets导入一个0，1类别的数据集并进行标准化，代码如下

from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_breast_cancer
cancer = load_breast_cancer()
x = cancer.data
y = cancer.target.reshape(-1,1)

sc_x = StandardScaler()
sc_y = StandardScaler()
X = sc_x.fit_transform(x)
Y = sc_y.fit_transform(y)

下面实现Logistic回归算法部分，其中有四大部分，Sigmoid函数，梯度上升训练，预测以及输出错误率

class logistic():
    def __init__(self):
        self.X = 0 
        self.y = 0
        self.W = 0
        self.y_hat = 0

    def sigmoid(X,W):
            z = np.dot(X,W) 
            g = 1/(1+np.exp(-z))
            return g
    def fit(self,X,y,alpha=0.1,iternum = 100):
        X = np.matrix(X)
        y = np.matrix(y)
        self.X = X #类定义
        self.y = y

        n = X.shape[0]
        k = X.shape[1]
        W = np.random.normal(size=(k,1)) #初始化W权重向量

        for i in range(iternum):
            rand_index = np.random.randint(n) #随机选取一个索引
            x_i = X[rand_index,:]
            y_i = y[rand_index,0]
            y_hat = sigmoid(x_i,W)
            gradient = x_i.T*(y_i-y_hat)
            W += alpha*gradient
        self.W = W

    def predict(self,X):
        g = sigmoid(X,self.W)
        g[np.where(g>0.5)] = 1
        g[np.where(g<=0.5)] = 0
        self.y_hat = g
        return g

    def error(self):
        err = np.sum(np.abs(self.y_hat-self.y))/len(self.y_hat)
        return err

下面测试一下

clf = logistic()
clf.fit(X,y)
y_hat = clf.predict(X)
clf.error()

输出错误率为4.2%