各个模型解读专题(2) 逻辑回归

1. 回顾线性回归

$$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\dots {\theta }_n{x}_n={\theta }^{T}x$$

$$x=\left(\begin{array}{c}1\\ x_1\\ x_2\\ ⋮\\ x_n\end{array}\right),\theta =\left(\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right)$$

$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m{\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)}^2$$

$$Y=X\theta$$

2. 改进用于分类

• sigmoid

$$g(z)=\frac{1}{1+e^{-z}}$$

$$g^{\prime }(z)=g(z)(1-g(z))$$

结合线性模型,$z=x\theta$

$$h(x)=\frac{1}{1+e^{-x\theta }}$$

其中x为样本输入，h(x)为模型输出，可以理解为某一分类的概率大小。而θ为分类模型的要求出的模型参数。h(x)的值越小，而分类为0的概率越高，反之，值越大分类为1的概率越高。如果靠近临界点，则分类准确率会下降.

$$\begin{array}{c}P(y=1\mid x,\theta )=h(x)\\ P(y=0\mid x,\theta )=1-h(x)\end{array}$$

y的概率分布函数表达式用矩阵法表示，即为：

$$P(Y\mid X,\theta )=h(x{)}^y(1-h(x){)}^{1-y}$$

得到了y的概率分布函数表达式，就可以用似然函数最大化来求解需要的模型系数θ。为了方便求解，这里用对数似然函数最大化。其中

似然函数的代数表达式为

$$L(\theta )=\prod _{i=1}^m{\left(h\left(x^{(i)}\right)\right)}^{y^{(i)}}{\left(1-h\left(x^{(i)}\right)\right)}^{1-y^{(i)}}$$

$$\begin{array}{rlr}l(\theta )& =\log L(\theta )& =\sum _{i=1}^m{y}^{(i)}\log h\left(x^{(i)}\right)+\left(1-y^{(i)}\right)\log \left(1-h\left(x^{(i)}\right)\right)\end{array}$$

$$\begin{array}{rl}\frac{\partial }{\partial {\theta }_j}l(\theta )& =\left(y\frac{1}{g\left({\theta }^{T}x\right)}-(1-y)\frac{1}{1-g\left({\theta }^{T}x\right)}\right)\frac{\partial }{\partial {\theta }_j}g\left({\theta }^{T}x\right)\\ & =\left(y\frac{1}{g\left({\theta }^{T}x\right)}-(1-y)\frac{1}{1-g\left({\theta }^{T}x\right)}\right)g\left({\theta }^{T}x\right)\left(1-g\left({\theta }^{T}x\right)\right)\frac{\partial }{\partial {\theta }_j}{\theta }^{T}x\\ & =(y\left(1-g\left({\theta }^{T}x\right)-(1-y)g\left({\theta }^{T}x\right)\right)x_j\\ & =(y-h(x))x_j\end{array}$$

根据上面的推导设立恰当的损失函数来方便逻辑回归的算法求解:

$$J(\theta )=\frac{1}{m}l(\theta )$$

此时训练时的梯度下降方向即为：

$${\theta }_j:={\theta }_j-\alpha \frac{\partial }{\partial {\theta }_j}J(\theta )={\theta }_j+\alpha \frac{1}{m}\sum _{i=1}^m\left(y-h\left(x^{(i)}\right)\right)x_j^{(i)}$$

以上是不加正则项的推导，现实情况往往是需要加正则项的，具体加什么样的正则项需要根据具体任务判断。

3. 实践

3.1. 自定义

# %%
import numpy as np
import math


class Sigmoid:
    def function(self, x):
        return 1 / (1 + np.exp(-x))

    def derivative(self, x):
        return self.function(x) * (1 - self.function(x))


class LogisticRegression(object):
    def __init__(self, learning_rate=0.1):
        self.w = None
        self.sigmoid = Sigmoid()
        self.lr = learning_rate

    def fit(self, X, y, iterations=10000):
        lr = self.lr
        n_samples, n_features = np.shape(X)
        self.w = np.random.uniform(
            -1 / math.sqrt(n_features), 1 / math.sqrt(n_features), (n_features, )
        )
        print(f"w's shape : {self.w.shape}")
        for i in range(iterations):
            y_pred = self.sigmoid.function(X.dot(self.w))

            delta = y - y_pred
            grad = self.sigmoid.function(X.dot(self.w)) * (
                1 - self.sigmoid.function(X.dot(self.w))
            )
            # print(f"delta.shape:{delta.shape},gard.shape:{grad.shape} || gard:{grad}")
            self.w -= (lr * X.T.dot(-delta * grad))

    def predict(self, X):
        pred = self.sigmoid.function(X.dot(self.w))
        y_prediction = np.round(pred).astype("int")

        return y_prediction


# %%
from sklearn.datasets import load_iris

iris = load_iris()

X = iris.data
y = iris.target
# y = y.reshape(-1,1)
names = iris.feature_names


from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X[0:100], y[0:100], test_size=0.2, random_state=2021
)


# %%
model = LogisticRegression()
model.fit(X_train, y_train)


# %%

from sklearn.metrics import accuracy_score

y_pred = model.predict(X_test)

accuracy_score(y_test,y_pred)

3.2. 调用sklearn库

from sklearn.linear_model import LogisticRegression

lgr = LogisticRegression()
lgr.fit(X_train,y_train)
accuracy_score(y_test,lgr.predict(X_test))

其中的重要参数如下：

• penalty

默认是L2的正则化，也可以选择L1正则化。如果在实际项目中发现过拟合情况严重或者由于模型特征过多，从而想让模型系数稀疏化可以选择L1正则化。

• solver

优化算法，一共有下面4种选择：

1）linlinear：即坐标轴下降法；
2）lbfgs：拟牛顿法的一种；
3）newton-cg：牛顿法的一种；
4）sag：随机平均梯度下降，该算法每次仅使用了部分样本计算梯度迭代，适合样本数据量很大的情况，收敛速度优于SGD和其他3种算法。

需要注意的是正则化项的选择会影响损失函数优化算法的选择。4种优化算法newton-cg、lbfgs、liblinear和sag都适用于L2正则化，但是如果选择L1正则化，则只能选择liblinear。这是由于L1正则化的损失函数不是连续可导的，而另外3种优化算法时要求损失函数的一阶可导或者二阶可导。

• multi_class

分类方式，multi_class参数即指分类方式的选择，可以选择one-vs-rest（OvR）和many-vs-many（multinomial）。需要注意的是liblinear无法用于multinomial这种分类方式。

• 权重类型 class_weight

class_weight参数是指分类模型中各种类型的权重。如果在class_weight选择balanced，则会根据训练样本量来计算权重。某种类型样本量越多，则权重越低，样本量越少，则权重越高。提高了某种分类的权重，相比不考虑权重，会有更多的样本分类划分到高权重的类别中，从而可以解决样本不平衡的问题。sample_weight也可以用于解决这类问题。

4. plus 逻辑回归模型如何用于多分类

多元逻辑回归常见的有one-vs-rest（OvR）和many-vs-many（MvM）两种，OvR指每次将一个类的样例作为正例，所有其他类的样例作为反例来训练N个分类器，在测试时若仅有一个分类器预测为正类，则对应的类别标记作为分类结果。MvM指每次将若干个类作为正类，若干个其他类作为反类。

比较简单的多分类算法是，假设是n元分类模型，把所有第N类样本作为正例，其余作为反例， 在此基础上进行二元回归，得到第N类的分类模型，以此类推。

$$\begin{array}{c}P(Y=n\mid x)=\frac{e^{\left({\theta }_{n}x\right)}}{1+{\sum }_{n=1}^{N-1}{e}^{\left({\theta }_{n}x\right)}}n=1,2,\cdots ,N-1\\ P(Y=N\mid x)=\frac{e^{\left({\theta }_{n}x\right)}}{1+{\sum }_{n=1}^{N-1}{e}^{\left({\theta }_{n}x\right)}}\end{array}$$

5. 逻辑回归和线性回归的比较

逻辑回归的本质是线性回归，只不过套了一层壳用于分类。套的这层壳就是$sigmoid$函数。

• 优化函数

再者，线性回归的优化函数是MSE,最小二乘，而逻辑回归的优化目标是一个似然函数。

加粗样式

• 代码来源: gihub

---

6. 项目实践

导入

import matplotlib.pyplot as plt
import numpy as np
from io import BytesIO

import torch
import torch.nn.functional as F

import pandas as pd

print(torch.__version__)

数据集准备阶段

df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data')

print("dataset preview:\n{}".format(df.head()))
print(df.tail())
ds = np.lib.DataSource()
fp = ds.open('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data')

x = np.genfromtxt(BytesIO(fp.read().encode()), delimiter=',', usecols=range(2), max_rows=100)
y = np.zeros(100)
y[50:] = 1

np.random.seed(1)
idx = np.arange(y.shape[0])
np.random.shuffle(idx)
X_test, y_test = x[idx[:25]], y[idx[:25]]
X_train, y_train = x[idx[25:]], y[idx[25:]]
mu, std = np.mean(X_train, axis=0), np.std(X_train, axis=0)
X_train, X_test = (X_train - mu) / std, (X_test - mu) / std

fig, ax = plt.subplots(1, 2, figsize=(7, 2.5))
ax[0].scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1])
ax[0].scatter(X_train[y_train == 0, 0], X_train[y_train == 0, 1])
ax[1].scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1])
ax[1].scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1])
plt.show()

数据集准备的思路如下 :

• 下载现有的数据

我们看一下数据的构成

可以看到每个数据有四个特征, 一个分类结果. 这就是著名的鸢尾花数据集

• 为了方便可视化展示,选取两个特征作为分类依据.

x = np.genfromtxt(BytesIO(fp.read().encode()), delimiter=',', usecols=range(2), max_rows=100)
y = np.zeros(100)
y[50:] = 1

该数据集一共100项,前50项为 Iris-setosa花朵的数据, 后面50项为Iris-setosa的数据.

利用numpy从中读取数据,usecols = range(2)就是只取前两个特征作为X的值.

读取完成之后,X的效果如下:

对于y,我们将前50设置为0,后50设置为1就好.

之后,将数据进行乱序处理:

np.random.seed(1) # 将这个乱序记录下来,之后再次训练这个代码的时候,还是按照这个顺序排列
idx = np.arange(y.shape[0]) # 计数器
np.random.shuffle(idx) # 将100个数字乱序形成数组 idx 
X_test, y_test = x[idx[:25]], y[idx[:25]]  # 取乱序数组的前25个值作为索引,索引X的值作为测试集合
X_train, y_train = x[idx[25:]], y[idx[25:]] # 同理作为训练集

• 将数据正则化

为了让计算量变小,让不同的衡量指标在同一个数据维度内可以比较 ,需要将数据正则化.

(例如: 在金融预测中, 对于金融数据 几百万到几百亿, 和哪一天这种个位数需要一起处理,那么需要将数据都限制到足够小才行.一般采用的方案就是正则化)

mu, std = np.mean(X_train, axis=0), np.std(X_train, axis=0)
X_train, X_test = (X_train - mu) / std, (X_test - mu) / std

$\frac{X-\mu }{\sigma }$ 就是将数据正则化的常用公式.((概率论常识))

• 画图表示数据

  fig, ax = plt.subplots(1, 2, figsize=(7, 2.5))
  ax[0].scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1]) # 在第一幅图中画散点图,将训练集中所有的第一类数据,f1和f2分别作为横轴和纵轴,画出来.
  ax[0].scatter(X_train[y_train == 0, 0], X_train[y_train == 0, 1])
  ax[1].scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1])
  ax[1].scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1])
  plt.show()
  

模型构建阶段

选择设备

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")

是利用GPU还是CPU训练, 一般来说, 尤其对于大型模型来说,GPU是CPU训练速度的几十倍.所以,有钱真好.

设立模型

def custom_where(cond, x_1, x_2):
    return (cond * x_1) + ((1-cond) * x_2)

class LogisticRegression1():
    def __init__(self, num_features):
        self.num_features = num_features
        self.weights = torch.zeros(num_features, 1,
                                   dtype=torch.float32, device=device)
        self.bias = torch.zeros(1, dtype=torch.float32, device=device)

    def forward(self, x):
        linear = torch.add(torch.mm(x, self.weights), self.bias)
        probas = self._sigmoid(linear)
        return probas

    def backward(self, probas, y):
        errors = y - probas.view(-1) # X.view(-1)中的-1本意是根据另外一个数来自动调整维度。
        return errors

    def predict_labels(self, x):
        probas = self.forward(x)
        labels = custom_where(probas >= .5, 1, 0)
        return labels

    def evaluate(self, x, y):
        labels = self.predict_labels(x).float()
        accuracy = torch.sum(labels.view(-1) == y) / y.size()[0]
        return accuracy

    def _sigmoid(self, z):
        return 1. / (1. + torch.exp(-z))

    def _logit_cost(self, y, proba):
        tmp1 = torch.mm(-y.view(1, -1), torch.log(proba))
        tmp2 = torch.mm((1 - y).view(1, -1), torch.log(1 - proba))
        return tmp1 - tmp2

    def train(self, x, y, num_epochs, learning_rate=0.01):
        for e in range(num_epochs):
            #### Compute outputs ####
            probas = self.forward(x)

            #### Compute gradients ####
            errors = self.backward(probas, y)
            
            neg_grad = torch.mm(x.transpose(0, 1), errors.view(-1, 1))

            #### Update weights ####
            self.weights += learning_rate * neg_grad
            self.bias += learning_rate * torch.sum(errors)

            #### Logging ####
            print('Epoch: %03d' % (e + 1), end="")
            print(' | Train ACC: %.3f' % self.evaluate(x, y), end="")
            print(' | Cost: %.3f' % self._logit_cost(y, self.forward(x)))

设立模型的思路如下:

• 从上述数据中 可以看出,这是一个明显的线性分类任务

• 设置线性模型

  $y=XW+b$

  W就是一个权重矩阵. 我们数据中选取的特征是2个,那么X的形状就是$X_{75\times 2}$ , 那么$W_{2\times 1}$ . 关于权重参数的数量,一般是有几个特征 ,就 选取几个参数.

• 模型结构:

• 初始化

  • 特征数量
  • 初始权重值
  • 偏置

• 前向传播

  • $prob=sigmoid(XW+b)$

    $y=\{\begin{array}{ll}1& prob\ge 0.5\\ 0& prob<0.5\end{array}$

• 向后传播

  后向传播就是更新权值而已

  • 计算梯度

    $\nabla =X^T(y-prob)=X^{T}errors$

    其实其含义就是:

    假如我们需要用以下数据做预测
    $$\begin{array}{ccc}{f}_1& f_2& y\\ x_{1,1}& x_{1,2}& 1\\ x_{2,1}& x_{2,2}& 1\\ x_{3,1}& x_{3,2}& 0\end{array}$$

    我们对于第1项数据的预测值是0.7, 第二项数据的预测值是0.6, 第三项的预测值是0.8.那么

    $erro{r}_1=1-0.7=0.3$

    $erro{r}_2=1-0.6=0.4$

    $erro{r}_3=0-0.8=-0.8$

    我们计算的梯度就是如下公式
    $$\nabla =\left(\begin{array}{c}0.3\\ 0.4\\ -0.8\end{array}\right)$$

    这就是梯度,也是下一步参数更新的方向.

  • 参数更新
    $$\begin{gathered} W=W+\eta \nabla \\ b=b+\eta (erro{r}_1+erro{r}_2+erro{r}_3) \end{gathered}$$

    偏置的更新与参数的更新不同,此处的偏置是一个常数, 更新的时候需要将所有的error值加起来用于更新.

    $\eta$是更新步长, 是一个可以自由设定的,控制更新步伐的一个超参数

训练评估阶段

训练,就是将前向传播和后向传播的操作组合起来就OK,评估就是计算分类正确率

X_train_tensor = torch.tensor(X_train, dtype=torch.float32, device=device)
y_train_tensor = torch.tensor(y_train, dtype=torch.float32, device=device)

logr = LogisticRegression1(num_features=2)
logr.train(X_train_tensor, y_train_tensor, num_epochs=10, learning_rate=0.1)

print('\n  Model parameters:')
print('  Weights: %s' % logr.weights)
print('  Bias: %s' % logr.bias)

---

结果可视化显示

w, b = logr.weights, logr.bias

x_min = -2
y_min = ((-(w[0] * x_min) - b[0]) / w[1])

x_max = 2
y_max = ((-(w[0] * x_max) - b[0]) / w[1])

fig, ax = plt.subplots(1, 2, sharex=True, figsize=(7, 3))
ax[0].plot([x_min, x_max], [y_min, y_max])
ax[1].plot([x_min, x_max], [y_min, y_max])

ax[0].scatter(X_train[y_train == 0, 0], X_train[y_train == 0, 1], label='class 0', marker='o')
ax[0].scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1], label='class 1', marker='s')
ax[1].scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1], label='class 0', marker='o')
ax[1].scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1], label='class 1', marker='s')

ax[1].legend(loc='upper left')
plt.show()