机器学习之Logistic&Sofmax回归模型--python实现

——《 Python机器学习算法:原理，实现与案例》读书笔记

相关数据集：
链接：https://pan.baidu.com/s/1fBNk7mbqCczKawmOyk491Q
提取码：niz2

概述：

Logistic回归虽然名为回归，实为分类算法，用于处理二分类问题


Sotfmax回归是在Logistic回归的基础上推广得到的，用于处理多元分类问题

Logistic回归

---

线性分类器

先从简单的二维平面的二分类问题开始


如下图所示：
在二维平面内存在两种类别的样本，我们可以找到一条直线，将这两种样本尽可能好得划分开
$$假设直线方程为：$$
$$w^{T}x+b=0w,b\in R$$

$$直线上方的点都满足：w^{T}x+b>k$$
$$而直线下方的点都满足：w^{T}x+b<k其中k为常数$$

这意味着，我们可以根据找出的这条直线，判断样本的类别
$$设z为样本的类别(0,1)$$
$$z=\{\begin{array}{l}0,w^{T}x+b=<k\\ 1,w^{T}x+b=>k\end{array}$$

对于多维空间的分类问题，类似上述问题，我们只需要找到一个超平面，可以将所有样本尽可能好得划分开（类似SVM）

$$假设超平面方程为：$$
$$w^{T}x+b=0w\in R^n,b\in R$$
类似的，根据此函数的输出函数值，我们可以实现对样本的分类

logistic函数

对于二维平面的二分类问题，最理想的分类函数显然为单位阶跃函数：
$$H(z)=\{\begin{array}{l}1,z\ge 0\\ 0,z\le 0\end{array}$$
但此函数有一个严重的缺点：不连续–>不是处处可微(有些算法就不可以运用)


logistic函数便是一直常用的替代函数
$$\sigma (z)=\frac{1}{1+e^{-z}}$$

概率意义

logistic是一种sigmoid函数，其值域在（0,1）之间连续

函数的输出可以视为x条件下样本分类为y=1的条件概率
$$P(y=1|x)=\sigma (g(x))=\frac{1}{1+e^{w^T+b}}$$

对数概率函数

在统计学中，几率(odds) 定义为：事件发生的概率与事件不发生概率的比值

几率表示的是样本作为正例的相对可能性

对数几率：
$$log\frac{p}{1-p}$$对数几率大于0表明正例的概率大，大于0表明负例的概率大

logistic回归模型假设一个实例为正例的对数几率是 输入(x)的线性函数

$$log\frac{p}{1-p}=w^{T}x+b$$
反求p得：
$$p=\sigma (g(x))=\frac{1}{1+e^{-(w^T+b)}}$$

一阶导数

$$\frac{\mathrm{d}\sigma (z)}{\mathrm{d}z}=\sigma (z)(1-\sigma (z))$$

Logistic回归模型

Logistic回归模型假设函数为：
$$h_{w,b}(x)=\sigma (g(x))=\frac{1}{1+e^{-(w^{T}x+b)}}$$

通常将b作为w0纳入权向量w,同时为输入向量添加常数1作为x0:
$$\begin{gathered} w=(b,w_1,w_2,...,w_n{)}^T \\ x=(1,x_1,x_2,...,x_n{)}^T \end{gathered}$$
此时，logistic函数为：
$$h_w(x)=\sigma (g(x))=\frac{1}{1+e^{-(w^{T}x)}}$$

$$h_w(x)的输出为预测x为正例的概率,通过训练可以确定模型参数w,构建二元分类函数：$$
$$H(h_w(x))=\{\begin{array}{l}1,h_w(x)\ge 0.5\\ 0,h_w(x)\le 0.5\end{array}$$

极大似然法估计参数

对于给定的包含m个样本的数据集D,可以使用极大似然法估计w
$$易知，模型将输入实例x_i预测为y_i的概率为：$$

$$P(y=y_i|x_i;w)=h_w(x_i{)}^{y_i}(1-h_w(x_i{)}^{1-y_i})y_i=0,1$$

$$定义似然函数为:$$

$$L(w)=\prod _{i=1}^{m}P(y=y_i|x_i;w)$$

$$比较Low的写法：L(w)=h_w(x_i{)}^k(1-h_w(x_i){)}^{n-k}n，抽样总次数，k，抽样为y=1的次数$$

极大似然法估计参数w的核心思想：已知x,y(样本分布),选择w,使得观测的数据(样本分布)出现的概率最大

(个人理解：不同的w,直接导致样本的分布不同，根据已知的分布，选择合适的参数（w），使得对任意样本集的分布 近似为 已知样本集的分布)
$$参考--极大似然估计$$
$$参考--极大似然估计$$

$$\hat{w}=\underset{w}{argmax}L(w)$$

为了简化模型，我们将目标函数转化为：
$$\hat{w}=\underset{w}{argmax}ln(L(w))$$对数似然函数（对数函数：积变和，单调递增函数）
$$l(w)=ln(L(w))=\sum _{i=1}^m{y}_{i}ln(h_w(x_i))+(1-y_i)ln(1-h_w(x_i))$$

梯度下降更新

对于Logistic回归模型，定义其损失函数为：
$$J(w)=-\frac{1}{m}\sum _{i=1}^m{y}_{i}ln(h_w(x_i))+(1-y_i)ln(1-h_w(x_i))$$
损失函数最小值与对数似然函数最大值等价

---

计算损失函数梯度：
$$计算J(w)对分量w_i的偏导数$$
$$\frac{\partial J(w)}{\partial w_j}=-\frac{1}{m}\frac{\partial }{\partial w_j}\sum _{i=1}^m{y}_{i}ln(h_w(x_i))+(1-y_i)ln(1-h_w(x_i))$$
$$=-\frac{1}{m}\sum _{i=1}^m{y}_i\frac{\partial }{\partial w_j}ln(h_w(x_i))+(1-y_i)\frac{\partial }{\partial w_j}ln(1-h_w(x_i))$$
$$=-\frac{1}{m}\sum _{i=1}^m{y}_i\frac{1}{h_w(x_i)}\frac{\partial h_w(x_i)}{z_i}\frac{\partial z_i}{w_j}+(1-y_i)\frac{1}{1-h_w(x_i)}\frac{-\partial h_w(x_i)}{z_i}\frac{\partial z_i}{w_j}$$
$$=-\frac{1}{m}\sum _{i=1}^m(y_i-h_w(x_i))\frac{\partial z_i}{\partial w_j}$$
$$=\frac{1}{m}\sum _{i=1}^m(h_w(x_i)-y_i)x_i$$

$$h_w(x_i)-y_i可解释为模型预测x_i为正例的概率与实际类别之间的误差$$

对于随机梯度下降算法，相应的梯度计算公式为：

$$▽J(w)=(h_w(x_i)-y_i)x_i$$

$$模型参数w的更新公式为：w:=w-\eta ▽J(w),\eta 为学习率$$

算法流程：

数据处理

循环

计算内积
概率转换
计算损失
求解梯度
更新参数

判断条件

import numpy as np

class LogisticRegression:
    def __init__(self, n_iter=200, eta=1e-3, tol=None):
        # 训练迭代次数
        self.n_iter = n_iter
        #学习率
        self.eta = eta
        #误差变化阈值
        self.tol = tol
        #模型参数
        self.w = None
    
    def _preprocess_data_X(self, X):
        '''数据预处理 '''
        
        m, n = X.shape
        X_ = np.empty((m, n+1))
        X_[:, 0] = 1 
        X_[:, 1:] = X
        
        return X_
    
    def _z(self, X, w):
        '''g(x)函数，计算内积(x与w)'''
        return np.dot(X, w)
    
    def _sigmod(self, z):
        '''Logistic函数'''
        return 1./(1. + np.exp(-z))
    
    def _predict_proba(self, X, w):
        '''h(x)函数，预测y=1的概率'''
        z = self._z(X, w)    
        return self._sigmod(z)
    
    def _loss(self, y, y_proba):
        '''计算损失'''
        m = y.size
        p = y_proba * (2 * y - 1) + (1 - y) 
        return -np.sum(np.log(p)) / m 
    
    def _gradient(self, X, y, y_proba):
        '''计算梯度'''  
        return np.matmul(y_proba-y, X) / y.size
    
    def _gradient_descent(self, w, X, y):
        '''梯度下降算法'''
        if self.tol is not None:
            loss_old = np.inf 
            
        for step_i in range(self.n_iter):
            y_proba = self._predict_proba(X, w) 
            
            loss = self._loss(y, y_proba)
            print('%4i Loss: %s' % (step_i, loss))
            
            if self.tol is not None:
                if loss_old - loss < self.tol:
                    break
                
            loss_old = loss 
        
            grad = self._gradient(X, y, y_proba)

            w -= self.eta * grad
            
    def train(self, X_train, y_train):
        '''训练'''
        #预处理
        X_train = self._preprocess_data_X(X_train)
        
        #初始化参数向量
        _, n = X_train.shape
        self.w = np.random.random(n) * 0.5
        
        #梯度下降训练
        self._gradient_descent(self.w, X_train, y_train)
    
    def predict(self, X): 
        '''预测''' 
        X = self._preprocess_data_X(X) 
        
        y_pred = self._predict_proba(X, self.w)
        
        # p>0.5 y=1
        return np.where(y_pred >= 0.5, 1, 0)

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# 数据导入
X = np.genfromtxt('wine.data', delimiter=',', usecols=range(1, 14))
y = np.genfromtxt('wine.data', delimiter=',', usecols=0)
idx = (y != 3)
X = X[idx]
y = y[idx]
y -= 1

#划分
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

#标准化
ss = StandardScaler()
ss.fit(X_train)
StandardScaler(copy=True, with_mean=True, with_std=True)
X_train_std = ss.transform(X_train)
X_test_std = ss.transform(X_test)

# 模型处理
clf = LogisticRegression(n_iter=2000, eta=0.01, tol=0.0001)
clf.train(X_train_std, y_train)

# 模型检验
from sklearn.metrics import accuracy_score

y_pred = clf.predict(X_test_std)
accuracy = accuracy_score(y_test, y_pred)
accuracy

0.9743589743589743

Softmax回归

Logistic回归只能处理二元分类问题，在此基础上推广得到的Sotfmax回归可以处理多元分类问题。Sotfmax回归也被称为多元Logistic回归。

Sotfmax函数

假设有k个分类，Sotfmax对实例x的类别进行预测时，需分别计算x属于每一个类别的概率，因此每个类别拥有各自独立的线性函数：
$$z_j=g_j(x)=w_j^{T}x$$

$$W=\left[\begin{array}{c}{w}_1^T\\ w_2^T\\ ⋮\\ w_k^T\end{array}\right]$$

---

$$定义Sotfmax回归的g(x)函数为：$$
$$z=g(x)=Wx=\left[\begin{array}{c}{z}_1\\ z_2\\ ⋮\\ z_k\end{array}\right]$$

定义Softmax函数为：
$$\sigma (z{)}_j=\frac{e^{z_j}}{{\sum }_{k=1}^K{e}^{z_k}}$$
Softmax函数的输出为：
$$\sigma (z)=\left[\begin{array}{c}\sigma (z{)}_1\\ \sigma (z{)}_2\\ ⋮\\ \sigma (z{)}_k\end{array}\right]$$

---

Sotfmax回归模型假设函数为：
$$h_w(x)=\sigma (g(x))=\frac{1}{{\sum }_{k=1}^K{e}^{w_k^{T}x}}\left[\begin{array}{c}{e}^{w_1^{T}x}\\ e^{w_2^{T}x}\\ ⋮\\ e^{w_K^{T}x}\end{array}\right]$$

Sotfmax函数的输出是预测x为各类别的概率，通过训练确定了模型参数W,便可构建多元分类函数：
$$H(h_w(x))=\underset{k}{argmax}{h}_W(x{)}_k=\underset{k}{argmax}(W_k^{T}x)$$

梯度下降更新公式

Sotfmax回归模型的损失函数为：
$$J(w)=-\frac{1}{m}\sum _{i=1}^m\sum _{j=1}^{K}I(y_i=j)ln{h}_W(x_i)$$

$$I为指示函数，当y_i=j时，I=1，否则I=0$$

这里直接给出梯度的计算公式：
$${▽}_{w_j}J(W)=\frac{1}{m}\sum _{i=1}^m(h_W(x_i{)}_j-I(y_i=j))x_i$$

$$h_W(x_i{)}_j-I(y_i=j)可以解释为模型预测x_i为第j类别的概率与其是否为第j类别之间的误差$$

对于随机梯度下降算法，每次只使用一个样本来计算梯度，相应梯度计算公式为：
$${▽}_{w_j}J(W)=\sum _{i=1}^m(h_W(x_i{)}_j-I(y_i=j))x_i$$

$$w_j的更新公式为：w_j:=w_j-\eta {▽}_{w_j}J(W)$$
即：
$$W:=W-\eta \left[\begin{array}{c}{▽}_{w_1}J(w{)}^T\\ {▽}_{w_2}J(w{)}^T\\ ⋮\\ {▽}_{w_k}J(w{)}^T\end{array}\right]$$

import numpy as np
import random

class SoftmaxRegression:
    def __init__(self, n_iter=200, eta=1e-3, tol=None):
        # 训练迭代次数
        self.n_iter = n_iter
        # 学习率
        self.eta = eta
        # 误差变化阈值
        self.tol = tol
        # 模型参数W
        self.W = None
        
    def _z(self, X, W):
        '''g(x)函数：计算X与W内积'''
        if X.ndim == 1 :
            return np.dot(W, X) 
        return np.matmul(X, W.T)
    
    def _softmax(self, Z):
        '''softmax函数'''
        E = np.exp(Z) 
        if Z.ndim == 1:
            return E / np.sum(E) 
        return E / np.sum(E, axis=1, keepdims=True)
    
    def _predict_proba(self, X, W): 
        '''h(x)函数，预测y为各类别的概率'''
        Z = self._z(X, W)
        return self._softmax(Z)
    
    def _loss(self, y, y_proba):
        '''计算损失'''
        m = y.size 
        # 获取每个样本对应种类的概率，如果这个样本为1类，就选取其对应1类的概率
        p = y_proba[range(m), y]
        return -np.sum(np.log(p)) / m 
    
    def _gradient(self, xi, yi, yi_proba):
        '''计算梯度'''
        K = yi_proba.size
        y_bin = np.zeros(K)
        # 对应样本类型为1
        y_bin[yi] = 1 
        
        return (yi_proba - y_bin)[:,None] * xi
    
    def _stochastic_gradient_descent(self, W, X, y): 
        '''随机梯度下降算法''' 
        if self.tol is not None:
            loss_old = np.inf
            end_count = 0
            
        # 核心代码
        m = y.size
        idx = np.arange(m)
        for step_i in range(self.n_iter):
            # 计算损失
            y_proba = self._predict_proba(X, W) 
            loss = self._loss(y, y_proba)
            print('%4i Loss: %s '%(step_i, loss)) 
            
            if self.tol is not None:
                # 随机梯度下降的loss曲线起伏较大，连续多次低于阈值，停止算法
                if loss_old - loss < self.tol:
                    end_count += 1 
                    if end_count == 5:
                        break
                else:
                    end_count = 0
                
                loss_old = loss 
            
            # 每一轮迭代前，打乱数据集
            np.random.shuffle(idx) 
            for i in idx: 
                yi_proba = self._predict_proba(X[i], W)
                
                grad = self._gradient(X[i], y[i], yi_proba)
                # 根据每一个样本更新梯度
                W -= self.eta * grad
                
    def _preprocess_data_X(self, X): 
        m, n = X.shape
        X_ = np.empty((m, n+1)) 
        X_[:, 0] = 1 
        X_[:, 1:] = X
        
        return X_ 
    
    def train(self, X_train, y_train):
        
        X_train = self._preprocess_data_X(X_train) 
        
        k = np.unique(y_train).size
        
        _, n = X_train.shape
        self.W = np.random.random((k, n)) * 0.05
        
        self._stochastic_gradient_descent(self.W, X_train, y_train)
        
    def predict(self, X): 
        X = self._preprocess_data_X(X) 
        Z = self._z(X, self.W)
        return np.argmax(Z, axis=1) 

X = np.genfromtxt('wine.data', delimiter=',', usecols=range(1, 14))
y = np.genfromtxt('wine.data', delimiter=',', usecols=0)
y -= 1 
y=y.astype(np.int32)

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
ss = StandardScaler()
ss.fit(X_train)
StandardScaler(copy=True, with_std=True, with_mean=True)
X_train_std = ss.transform(X_train)
X_test_std = ss.transform(X_test)

clf = SoftmaxRegression(n_iter=2000, eta=0.01, tol=0.0001)
clf.train(X_train_std, y_train)

# 模型检验
from sklearn.metrics import accuracy_score
y_pred = clf.predict(X_test_std)
accuracy_score(y_test, y_pred)

Python机器学习算法:原理，实现与案例 – 刘硕