《统计学习方法》之逻辑斯蒂回归代码实现

最新推荐文章于 2022-10-17 20:47:58 发布

原创最新推荐文章于 2022-10-17 20:47:58 发布 · 1.2k 阅读

2 ·

CC 4.0 BY-SA版权

机器学习专栏收录该内容

12 篇文章

订阅专栏

本文深入解析逻辑斯蒂回归模型，解释对数几率回归的概念，介绍如何使用极大似然估计求解参数，并提供手写数字识别的代码实现。

逻辑斯蒂回归，又被称为对数几率回归。下面根据两个名字对模型进行解释。
先说为什么被称为对数几率回归，什么是几率？假设一个事件发生的概率为p则事件发生的几率为 $p1−p\frac {p}{1-p}$ ，那么对数几率也就表示为 $log(p1−p)log(\frac {p}{1-p})$ ，现在考虑二分类问题，假设标签为0，1，p表示预测标签为1的概率，也就是 $p = P (y = 1 ∣ X)$ ，我们用线性模型来拟合事件的对数几率，也就是
$log(p1−p)=log(P(y=1∣X)1−P(y=1∣X))=wx+blog(\frac {p}{1-p})=log(\frac {P(y=1|X)}{1-P(y=1|X)})=wx+b$
通过上面的式子就得到了：
$P(y=1∣X)=ewx+b1+ewx+bP(y=1|X)=\frac {e^{wx+b}} {1+e^{wx+b}}$

$P(y=0∣X)=11+ewx+bP(y=0|X)=\frac {1} {1+e^{wx+b}}$
有了上面的式子就可以根据极大似然估计求解参数w了。
再说逻辑斯蒂回归，先计算出线性模型的值， $w x + b$ ，由于我们是要做分类问题，当然可直接选择一个阈值，通过将线性模型得到的值与阈值比较大小进行类别判断，但这样略显粗糙？(这样难以定义损失函数进行优化)。于是我们找来了逻辑函数
$P(t)=11+e−tP(t)=\frac {1} {1+e^{-t}}$

将线性模型的输出值作为t输入到逻辑函数函数中，将输出作为标签为1的概率，即：
$P(y=1∣X)=11+e−(wx+b)=ewx+b1+ewx+bP(y=1|X)=\frac {1}{1+e^{-(wx+b)}}=\frac {e^{wx+b}} {1+e^{wx+b}}$

实验注意事项：实验中使用的数据为手写数字数据集，每个像素点的取值范围为[0，255]，如果不对数据进行预处理，则在进行梯度更新时很容易出现溢出，针对这个问题有三种解决办法，具体见代码中的注释。

下面是代码实现部分：

#!/usr/bin/env python
# -*- encoding: utf-8 -*-
"""
Created on 2020/1/31 17:22
@author: phil
"""

from dataloader import load_data        # 自定义
import numpy as np
from random import shuffle
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.preprocessing import StandardScaler


def preprocess_data(images, labels):
    # 只保留数字为0，1的图片
    mask = (labels[:, 0] == 0)+(labels[:, 0] == 1)
    images = images[mask]
    labels = labels[mask]

    # 在输入数据中新增一列1，作为wx+b中的常数项，后面直接使用Wx表示Wx+b
    one_column = np.ones((images.shape[0], 1))
    images = np.c_[images, one_column]

    # 数据集分割，将20%数据作为测试集
    X_train, X_test, y_train, y_test = train_test_split(images, labels, test_size=0.2, random_state=42)

    # 如果直接使用原始数据，每个样本中的数据值取值范围为[0,255]，计算过程会使得exp(wx)过大导致结算结果出错
    # 有三种解决办法
    # 方法1：将p(y=1|x)=exp(wx)/(1+exp(wx))改写为p(y=1|x)=1/(1+exp(-wx))，见fit函数
    # 方法2：将输入数据归一化 X = X / 255，使得属性取值范围缩到[0,1]
    # X_train = X_train / 255
    # X_test = X_test / 255
    # 方法3：将输入数据标准化，代码如下：注意测试集的标准化参数来自训练集
    # stdSca = StandardScaler()
    # X_train = stdSca.fit_transform(X_train)
    # X_test = stdSca.transform(X_test)

    return X_train, X_test, y_train, y_test


class LogisticRegression():
    def __init__(self):
        self.W = None

    def fit(self, X, y, learning_rate=0.01, epoch=10):
        # X每一行是一个样本
        n, dim = X.shape
        # 初始化参数
        self.W = np.zeros((dim, 1))
        # 这里使用随机梯度下降算法，每次选择一个样本对参数进行更新
        for i in range(epoch):
            # 打乱训练集顺序
            index = list(range(n))
            shuffle(index)
            X, y = X[index], y[index]
            for j in range(n):
                Xj = X[j].reshape(1, -1)
                yj = y[j]
                # 原始计算公式
                # e_Wx = np.exp(Xj.dot(self.W))
                # p1 = e_Wx / (1+e_Wx)
                # 避免数值溢出
                p1 = 1 / (1 + np.exp(-Xj.dot(self.W)))
                self.W += learning_rate * (yj - p1) * Xj.T

    def predict(self, X):
        preds = []
        for i in range(len(X)):
            Xi = X[i]
            # 见fit函数
            # e_Wx = np.exp(Xi.dot(self.W))
            # p1 = e_Wx / (1+e_Wx)
            p1 = 1 / (1 + np.exp(-Xi.dot(self.W)))
            if p1 >= 0.5:
                preds.append(1)
            else:
                preds.append(0)
        return np.array(preds)


if __name__ == "__main__":
    images, labels = load_data()
    X_train, X_test, y_train, y_test = preprocess_data(images, labels)
    print("train size", len(X_train))
    print("test size", len(X_test))
    model = LogisticRegression()
    model.fit(X_train, y_train)
    print("accuracy score:", accuracy_score(model.predict(X_test), y_test))

"""
Output:
    train size 10132
    test size 2533
    accuracy score: 0.9976312672720095
"""