05.逻辑斯蒂回归

算法思路

线性回归可以对数据进行拟合，但是对于分类问题，就不是怎么$work$。
我们可以借用同样的思路，取拟合他的其他特征，比如概率，但是概率$\in [0,1]$，所以这里拟合数据的对数几率函数，即：
$$\ln \frac{p^{\prime }}{1-p^{\prime }}=w^{T}x+b$$
这就是逻辑斯蒂回归，同时也被人称为对数几率回归。
然后我们采用最大似然估计，来估计出$w,b$的值即可。

模型定义

输入： X = { x ( 1 ) , x ( 2 ) , … , x ( K ) } , y = { y ( 1 ) , y ( 2 ) , … , y ( K ) } X=\{x^{(1)}, x^{(2)}, \dots, x^{(K)}\}, y=\{y^{(1)}, y^{(2)}, \dots, y^{(K)}\} X={x(1),x(2),…,x(K)},y={y(1),y(2),…,y(K)}，其中 x ( i ) ∈ R D , y ( i ) ∈ { 0 , 1 } x^{(i)}\in R^D,y^{(i)}\in\{0,1\} x(i)∈RD,y(i)∈{0,1}。
输出： $w,b$
模型： ${\prod }_{k=1}^K({p^{\prime (k)}}^{y^{(k)}}{(1-p^{\prime (k)})}^{1-y^{(k)}})$

模型推导

$$\begin{array}{rl}L(w,b|X,y)& =\prod _{k=1}^K({p^{\prime (k)}}^{y^{(k)}}{(1-p^{\prime (k)})}^{1-y^{(k)}})\\ \ln L(w,b|X,y)& =\sum _{k=1}^K{y}^{(k)}\ln p^{\prime (k)}+(1-y^{(k)})\ln (1-p^{\prime (k)})\\ & =\sum _{k=1}^K{y}^{(k)}\ln \frac{p^{\prime (k)}}{1-p^{\prime (k)}}+\ln (1-p^{\prime (k)})\\ & =\sum _{k=1}^K{y}^{(k)}(w^T{x}^{(k)}+b)-\ln (1+\exp (w^T{x}^{(k)}+b))\end{array}$$
对其求偏导得：
$$\begin{array}{rl}\frac{\partial \ln L(w,b|X,y)}{\partial w}& =\sum _{k=1}^K{y}^{(k)}{x}^{(k)}-\frac{\exp (w^T{x}^{(k)}+b)x^{(k)}}{1+\exp (w^T{x}^{(k)}+b)}\\ & =\sum _{k=1}^K(y^{(k)}-p^{\prime (k)})x^{(k)}\end{array}$$
$$\begin{array}{rl}\frac{\partial \ln L(w,b|X,y)}{\partial b}& =\sum _{k=1}^K{y}^{(k)}-\frac{\exp (w^T{x}^{(k)}+b)}{1+\exp (w^T{x}^{(k)}+b)}\end{array}$$

代码实现

% acc: 99.16%
import numpy as np


def load_data(file_name):
    data, label = [], []
    fr = open(file_name, 'r')
    for line in fr.readlines():
        line = line.strip().split(',')
        if int(line[0]) == 0:
            label.append(1)
        else:
            label.append(0)
        data.append([int(num)/255 for num in line[1:]])
    return data, label


def logistic(data, label, epochs=300, lr=0.01):
    data, label = np.array(data).T, np.array(label)
    dims, nums = data.shape
    weights, bias = np.zeros((dims, )), 0
    for epoch in range(epochs):
        for k in range(nums):
            weights_mul_x = np.dot(weights.T, data[:, k]) + bias
            p = np.exp(weights_mul_x) / (1 + np.exp(weights_mul_x))
            weights += lr * (label[k] - p) * data[:, k]
            bias += lr * (label[k] - p)
    return weights, bias


def predict(data, weights, bias):
    data = np.array(data).T
    label = []
    for k in range(data.shape[1]):
        weights_mul_x = np.dot(weights.T, data[:, k]) + bias
        p = np.exp(weights_mul_x) / (1 + np.exp(weights_mul_x))
        label.append(1 if p > 0.5 else 0)
    return label


if __name__ == '__main__':
    train_data, train_label = load_data('../data/mnist_train.csv')
    test_data, test_label = load_data('../data/mnist_test.csv')
    weights, bias = logistic(train_data, train_label)
    y_hat = predict(test_data, weights, bias)
    error = 0
    for y1, y2 in zip(test_label, y_hat):
        if y1 != y2:
            error += 1
    print('acc: {}%'.format(100 - error / len(y_hat) * 100))