机器学习Python实战-第三章-分类问题-2.逻辑回归算法-EW帮帮网

3.2.1 原理简介

1.普通线性回归

2.广义线性模型

3.逻辑回归

前面都是算法公式原理，代码实操部分在实战和实验。

3.2.1 原理简介

逻辑回归是一种常见的广义线性模型，线性模型中的“线性”就是一系列一次特征的线性组合，在二维空间中是一条直线，在三维空间中是一个平面，如果推广到n维空间，就可以理解为广义线性模型。线性模型（linear model）的形式为：

$f(x)=w\cdot x+b$

其中， $x=(x^{(1)},x^{(2)},\dots ,x^{(n)})^T$ 是用列向量表示的样本，该样本有 $n$ 种特征，用 $x^{(i)}$ 表示样本 $x$ 的第 $i$ 个特征。 $w=(w^{(1)},w^{(2)},\dots ,w^{(n)})^T$ 为每个特征对应的权重生成的权重向量，权重向量直观的表达了各个特征在预测中的重要性。

1.普通线性回归

线性回归是一种回归分析技术，回归分析本质上就是一个函数估计的问题（函数估计包括参数估计和非参数估计两类），就是找出因变量和自变量之间的因果关系。回归分析的因变量应该是连续变量，若因变量为离散变量，则问题转换为分类问题，回归分析是一个有监督学习的问题。

给定数据集 $T=\{(x_1, y_1),(x_2, y_2),\dots,(x_N, y_N)\},x_i\in \chi \subseteq R^n,y_i\in Y \subseteq R,i=1,2,\dots ,N$ ，其中 $x=[x^{(1)},x^{(2)},\dots ,x^{(n)}]^T$ ，学习的模型为：

$f(x)=w\cdot x+b$

下面根据已知的数据集 $T$ 来计算参数 $w$ 和 $b$ 。

对于给定的样本 $x_i$ ，其预测值为 $\hat{y}_i=f(x_i)=w\cdot x_i+b$ 。采用平方损失函数，则在训练集 $T$ 上，模型的损失函数为:

$L(f)=\sum_{i=1}^{N}(\hat{y}_i-y_i)^2=\sum_{i=1}^{N}(w\cdot x_i+b-y_i)^2$

要使损失函数最小，即：

$(w^{*}, b^{*})=arg\underset{w,b}{min}\sum_{i=1}^{N}(w\cdot x_i+b-y_i)^2$

可以用梯度下降法来求解上述最优化问题的数值解，同时要对特征进行归一化处理，并利用最小二乘法来求解解析解。

令

$\tilde{w} =(w^{(1)},w^{(2)},\dots ,w^{(n)},b)^T=(w^T,b)^T$

$\tilde{x} =(x^{(1)},x^{(2)},\dots ,x^{(n)},1)^T=(x^T,1)^T$

$y=(y_1,y_2,\dots ,y_N)^T$

则有

${\textstyle \sum_{i=1}^{N}}(w\cdot x_i+b-y_i)^2=(y-(\tilde{x}_1,\tilde{x}_2,\dots,\tilde{x}_N)^T\tilde{w})^T(y-(\tilde{x}_1,\tilde{x}_2,\dots,\tilde{x}_N)^T\tilde{w})$

令

$x=(\tilde{x}_1,\tilde{x}_2,\dots,\tilde{x}_N)^T= \begin{bmatrix} \tilde{x}_1^T\\ \tilde{x}_2^T\\ \vdots\\ \tilde{x}_N^T \end{bmatrix}= \begin{bmatrix} x_1^{(1)} & x_1^{(2)} & \dots & x_1^{(n)} & 1\\ x_2^{(1)} & x_2^{(2)} & \dots & x_2^{(n)} & 1\\ \vdots & \vdots & \ddots & \vdots & 1\\ x_N^{(1)} & x_N^{(2)} & \dots & x_N^{(n)} & 1 \end{bmatrix}$

则有

$\tilde{w}^*=arg\underset{\tilde{w}}{min}(y-x\tilde{w})^T(y-x\tilde{w})$

令 $E_{\tilde{w}}=(y-x\tilde{w})^T(y-x\tilde{w})$ ，求它的极小值。对 $\tilde{w}$ 求导，令导数为0，得到以下解析解:

$\frac{\partial E_{\tilde{w}}}{\partial \tilde{w}}=2x^T(x\tilde{w}-y)=0\Rightarrow x^Tx\tilde{w}=x^Ty$

(1) 当 $x^Tx$ 为满秩矩阵或者正定矩阵时，可得

$\tilde{w}^*=(x^Tx)^{-1}x^Ty$

其中， $(x^Tx)^{-1}$ 为 $x^Tx$ 的逆矩阵。于是得到多元线性回归模型为

$f(\hat{x}_i)=\hat{x}_i^T\hat{w}^*$

(2) 当 $x^Tx$ 不是满秩矩阵时，比如 $N<n$ （样本数量小于特征种类的数量），根据 $x$ 的秩小于或等于 $(N,n)$ 中的最小值，即小于或等于 $N$ （矩阵的秩一定小于或等于矩阵的行数和列数）；而矩阵 $x^Tx$ 是 $n\times n$ 大小的，它的秩一定小于或等于 $N$ ，因此不是满秩矩阵。此时存在多个解析解。常见的作法是引入正则化项，如 $L_1$ 正则化或者 $L_2$ 正则化。以 $L_2$ 正则化为例:

$\tilde{w}^*=arg\underset{\tilde{w}}{min}[(y-x\tilde{w})^T(y-x\tilde{w})+\lambda\|\tilde{w}\|_2^2]$

其中， $\lambda >0$ 时调整正则化项与均方误差的比例， $\|\dots \|_2$ 为 $L_2$ 范数。

2.广义线性模型

考虑单调可导函数 $h(\cdot)$ ，令 $h(y)=w^Tx+b$ ，这样得到的模型称为广义线性模型（generalized linear model）。广义线性模型的一个典型的例子就是对数线性回归。当 $h(\cdot)=\ln(\cdot)$ 时的广义线性模型就是对数线性回归，即

$\ln y=w^Tx+b$

它是通过 $exp(w^Tx+b)$ 来拟合 $y$ 的。它虽然被称为广义线性回归，但实质上是非线性的。

3.逻辑回归

上述的学习方法都是使用线性模型进行回归学习的，而线性模型也可以用作分类。考虑二分类问题，给定数据集 $T=\{(x_1, y_1),(x_2, y_2),\dots,(x_N, y_N)\},x_i\in \chi \subseteq R^n,y_i\in \{0,1\},i=1,2,\dots ,N$ ，其中 $x_i=(x_i^{(1)},x_i^{(2)},\dots ,x_i^{(n)})^T$ 。需要知道 $P(y|x)$ ，这里用条件概率的原因是：预测时都已知 $x$ ，然后需要判断此时对应的 $y$ 值。

考虑到 $w\cdot x+b$ 取值是连续的，因此它不能拟合离散变量。可以考虑用它来拟合条件概率 $P(y=1|x)$ ，因为概率的取值也是连续的。但是对于 $w\ne 0$ （若等于零向量则没有求解的价值）， $w\cdot x+b$ 的取值范围是 $-\infty \sim +\infty$ ，不符合概率取值在范围 $[0,1]$ 的要求，因此考虑采用广义线性模型，最理想的是单位阶跃函数：

$P(y=1|x)= \left\{\begin{matrix} 0, & z<0\\ 0.5, & z=0\\ 1, & z>0 \end{matrix}\right.,\ z=w\cdot x+b$

但是阶跃函数不满足单调可导的性质。因此，需要寻找一个可导的、与阶跃函数相似的函数。对数概率函数（logistic function）就是这样一个替代函数：

$P(y=1|x)=\frac{1}{1+e^{-z}},z=w\cdot x+b$

由于 $P(y=0|x)=1-P(y=1|x)$ ，则有

$\ln \frac{P(y=1|x)}{P(y=0|x)}=z=w\cdot x+b$

$\frac{P(y=1|x)}{P(y=0|x)}$ 表示样本为正例的可能性与为反例的可能性之比，称为概率（odds），反映了样本作为正例的相对可能性。概率的对数称为对数概率（log odds，又称logit）。

下面给出逻辑回归模型参数估计。给定训练数据集 $T=\{(x_1, y_1),(x_2, y_2),\dots,(x_N, y_N)\}$ ，其中 $x_i\in R^n,y_i\in \{0,1\}$ 。模型估计的原理使用极大似然法估计模型参数，

为了便于讨论，将参数吸收进中，即令

$\tilde{w} =(w^{(1)},w^{(2)},\dots ,w^{(n)},b)^T\in R^{n+1}$

$\tilde{x} =(x^{(1)},x^{(2)},\dots ,x^{(n)},1)^T\in R^{n+1}$

$P(Y=1|\tilde{x})=\pi(\tilde{x})=\frac{exp(\tilde{w}\cdot \tilde{x})}{1+exp(\tilde{w}\cdot \tilde{x})}$

$P(Y=0|\tilde{x})=1-\pi(\tilde{x})$

则似然函数为

$\prod_{i=1}^{N}[\pi(\tilde{x_i})]^{y_i}[1-\pi(\tilde{x_i})]^{1-y_i}$

对数似然函数为

$L(\tilde{w})=\sum_{i=1}^N[y_i\log\pi(\tilde{x})+(1-y_i)\log(1-\pi(\tilde{x}))]$

$=\sum_{i=1}^N[y_i\log\frac{\pi(\tilde{x})}{1-\pi(\tilde{x})}+\log(1-\pi(\tilde{x}))]$

又由于

$\pi(\tilde{x})=\frac{exp(\tilde{w}\cdot \tilde{x})}{1+exp(\tilde{w}\cdot \tilde{x})}$

因此有

$L(\tilde{w})=\sum_{i=1}^N[y_i(\tilde{w}\cdot \tilde{x})-\log(1+exp(\tilde{w}\cdot \tilde{x}))]$

对 $L(\tilde{w})$ 求极大值，得到 $\tilde{w}$ 的估计值。设估计值为 $\tilde{w}^*$ ，则逻辑回归模型为

$P(Y=1|X=\tilde{x})=\frac{exp(\tilde{w}^*\cdot \tilde{x})}{1+exp(\tilde{w}^*\cdot \tilde{x})}$

$P(Y=0|X=\tilde{x})=\frac{1}{1+exp(\tilde{w}^*\cdot \tilde{x})}$

以上讨论的都是二分类的逻辑回归模型，可以推广到多分类逻辑回归模型。设离散型随机变量 $Y$ 的取值集合为 $\{1,2,\dots,K\}$ ，则多分类逻辑回归模型为

$P(Y=k|\tilde{x})=\frac{exp(\tilde{w}^*\cdot \tilde{x})}{1+\sum_{k=1}^{K-1}exp(\tilde{w}^*\cdot \tilde{x})},\ k=1,2,\dots,K-1$

$P(Y=K|\tilde{x})=\frac{1}{1+\sum_{k=1}^{K-1}exp(\tilde{w}^*\cdot \tilde{x})},\ \tilde{x}\in R^{n+1},\tilde{w}_k\in R^{n+1}$

其参数估计方法与二分类逻辑回归模型类似。

3.2.2 算法步骤

输入：数据集 $T=\{(x_1, y_1),(x_2, y_2),\dots,(x_N, y_N)\},x_i\in \chi \subseteq R^n,y_i\in Y \subseteq R,i=1,2,\dots ,N$ ，正则化项系数 $\lambda >0$ 。

输出：

$f(x)=w\cdot x+b$

算法步骤：

令

$\tilde{w} =(w^{(1)},w^{(2)},\dots ,w^{(n)},b)^T=(w^T,b)^T$

$\tilde{x} =(x^{(1)},x^{(2)},\dots ,x^{(n)},1)^T=(x^T,1)^T$

$y=(y_1,y_2,\dots ,y_N)^T$

计算

$x=(\tilde{x}_1,\tilde{x}_2,\dots,\tilde{x}_N)^T= \begin{bmatrix} \tilde{x}_1^T\\ \tilde{x}_2^T\\ \vdots\\ \tilde{x}_N^T \end{bmatrix}= \begin{bmatrix} x_1^{(1)} & x_1^{(2)} & \dots & x_1^{(n)} & 1\\ x_2^{(1)} & x_2^{(2)} & \dots & x_2^{(n)} & 1\\ \vdots & \vdots & \ddots & \vdots & 1\\ x_N^{(1)} & x_N^{(2)} & \dots & x_N^{(n)} & 1 \end{bmatrix}$

优化求解

$\tilde{w}^*=arg\underset{\tilde{w}}{min}[(y-x\tilde{w})^T(y-x\tilde{w})+\lambda\|\tilde{w}\|_2^2]$

最终得到模型

$f(\hat{x}_i)=\hat{x}_i^T\tilde{w}^*$

3.2.3 实战

1.数据集

1.线性回归
在线性回归问题中，使用的数据集是sklearn自带的一个糖尿病病人的数据集。该数据集从糖尿病病人采样并整理后，特点如下：

数据集有442个样本。

每个样本有10个特征。

每个特征都是浮点数，数据的范围是 -0.2 ~ 0.2。

样本的目标为25~346的整数。

这里给出加载数据集的函数：
from sklearn import datasets
from sklearn.model_selection import train_test_split


def load_data():
    diabetes = datasets.load_diabetes()
    return train_test_split(diabetes.data, diabetes.target, test_size=0.25, random_state=0)
使用该数据集返回值是一个元组，元组依次是：训练样本集、测试样本集、训练样本集对应的标签值、测试样本集对应的标签值。

2.逻辑回归

为了测试逻辑回归模型的分类性能，此处选用经典的数据集：鸢尾花数据集。

2.sklearn实现

1.线性回归

测试线性回归模型，代码如下：

from sklearn import datasets, linear_model
from sklearn.model_selection import train_test_split
import numpy as np
import pytest


def load_data():
    diabetes = datasets.load_diabetes()
    return train_test_split(diabetes.data, diabetes.target, test_size=0.25, random_state=0)


@pytest.fixture
def data():
    """Fixture to load data."""
    return load_data()


def test_LinearRegression(data):
    X_train, X_test, y_train, y_test = data
    regression = linear_model.LinearRegression()
    regression.fit(X_train, y_train)
    print('\nCoefficients:%s, intercept:%.2f' % (regression.coef_, regression.intercept_))
    print("Residual sum of squares:%.2f" % np.mean((regression.predict(X_test) - y_test) ** 2))
    print('Score:%.2f' % regression.score(X_test, y_test))

代码结果如下：

测试集中预测结果的均方误差为3180.16，预测性能得分为0.36（该值越大越好，最大为1.0）。

2.逻辑回归

测试逻辑回归模型（使用的数据集是鸢尾花数据集，不是糖尿病数据集），代码如下：

def load_iris():
    iris = datasets.load_iris()
    x = iris.data
    y = iris.target
    return model_selection.train_test_split(x, y, test_size=0.3, random_state=1, shuffle=True, stratify=y)


@pytest.fixture
def data2():
    """Fixture to load data."""
    return load_iris()


def test_LogisticRegression(data2):
    X_train, X_test, y_train, y_test = data2
    regression = linear_model.LogisticRegression()
    regression.fit(X_train, y_train)
    print('\nCoefficients:%s, intercept:%s' % (regression.coef_, regression.intercept_))
    print('Score:%.2f' % regression.score(X_test, y_test))

代码结果如下

测试集中的预测结果性能得分为0.98，即预测准确率为98%。

下面考察multi_class参数对分类结果的影响。默认采用的是one-vs-rest策略，但是逻辑回归模型的原型就支持多分类，给出的测试函数如下：

def test_LogisticRegression_multinomial(data2):
    X_train, X_test, y_train, y_test = data2
    regression = linear_model.LogisticRegression(multi_class='multinomial', solver='lbfgs')
    regression.fit(X_train, y_train)
    print('\nCoefficients:%s, intercept:%s' % (regression.coef_, regression.intercept_))
    print('Score:%.2f' % regression.score(X_test, y_test))

代码结果如下

最后，考察参数C对分类模型的预测性能的影响。C是正则化项系数的倒数，它越小则正则化项的权重越大。测试函数如下：

def test_LogisticRegression_C(data2):
    X_train, X_test, y_train, y_test = data2
    Cs = np.logspace(-2, 4, num=100)
    scores = []
    for C in Cs:
        regression = linear_model.LogisticRegression(C=C)
        regression.fit(X_train, y_train)
        scores.append(regression.score(X_test, y_test))
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(Cs, scores)
    ax.set_xlabel(r"C")
    ax.set_ylabel(r"score")
    ax.set_xscale('log')
    ax.set_title(r"Logistic Regression")
    plt.show()

代码结果如下

测试结果如下图。可以看到随着C的增大（即正则化项减小），LogisticRegression的预测准确率上升。当C增大到一定程度（即正则化项减小到一定程度）时，LogisticRegression的预测准确率维持在较高的水准保持不变。

3.算法实现

为了使用逻辑回归模型对鸢尾花进行分类，此处选用经典的鸢尾花数据集。

现只取数据集Iris中的两个特征Sepal.length（花萼长度）和Petal.length（花瓣长度），定义为 $X(X1,X2)$ ，对应 y 分类中的两个类别（0，1），将根据 $X(X1,X2)$ 的值对鸢尾花进行分类。首先绘制这两个特征的散点图，代码如下。

from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
import numpy as np


iris = load_iris()
data = iris.data
target = iris.target
X = data[0:100, [0, 2]]
y = target[0:100]
label = np.array(y)
index_0 = np.where(label == 0)
plt.scatter(X[index_0, 0], X[index_0, 1], marker='x', color='b', label='0', s=15)
index_1 = np.where(label == 1)
plt.scatter(X[index_1, 0], X[index_1, 1], marker='o', color='r', label='1', s=15)
plt.xlabel('X1')
plt.ylabel('X2')
plt.legend(loc='upper left')
plt.show()

代码结果如下

接着编写一个逻辑回归模型的类，然后训练测试，计算损失函数（损失函数的本质是衡量“模型预估值”到“实际值”的距离）。注意损失函数值越小，模型越好，而且损失函数尽量是一个凸函数，便于收敛计算。逻辑回归模型预估的是样本属于某个分类的概率，其损失函数可以采用均方差、对数、概率等方法。计算损失函数的代码如下：

import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_iris


# Iris 数据集中的目标变量 y 是整数向量（类标签），
# 但逻辑回归模型期望 y 为二进制格式（0 和 1）以进行二进制分类。需要相应地转换标签。
class LogisticRegressionBinary(object):
    def __init__(self):
        self.W = None

    def train(self, X, y, lr=0.01, num_iters=5000):
        num_train, num_feature = X.shape
        self.W = 0.001 * np.random.randn(num_feature, 1).reshape((-1, 1))
        loss = []
        for i in range(num_iters):
            error, dW = self.compute_loss(X, y)
            self.W += - lr * dW
            loss.append(error)
            if i % 200 == 0:
                print('i= %d, error= %f' % (i, error))
        return loss

    def compute_loss(self, X, y):
        num_train = X.shape[0]
        h = self.output(X)
        loss = - np.sum((y * np.log(h) + (1 - y) * np.log(1 - h))) / num_train
        dW = X.T.dot(h - y) / num_train
        return loss, dW

    def output(self, X):
        g = np.dot(X, self.W)
        return self.sigmoid(g)

    def sigmoid(self, X):
        return 1 / (1 + np.exp(-X))

    def predict(self, X_test):
        h = self.output(X_test)
        return np.where(h >= 0.5, 1, 0)


# 加载 Iris 数据集
iris = load_iris()
X = iris.data
y = iris.target

# 鸢尾花数据集有三个类（0，1，2），筛选出类 0 和 1 的数据
# 对于二元分类，筛选类 0 和 1 的数据
binary_filter = y < 2
X = X[binary_filter]
y = y[binary_filter].reshape((-1, 1))

# 在 X 矩阵左侧添加全 1 的列，说明模型中的截距项
one = np.ones((X.shape[0], 1))
X_train = np.hstack((one, X))

# 训练 Logistic 回归模型，使用 Logistic 回归进行二进制分类。
classify = LogisticRegressionBinary()
loss = classify.train(X_train, y)

# 输出学习的权重
print("Learned weights:\n", classify.W)

# 绘制迭代的损失曲线
plt.plot(loss)
plt.xlabel('Iteration number')
plt.ylabel('Loss value')
plt.title('Loss curve for Logistic Regression')
plt.show()

（书上代码达不到给出的效果，具体原因以及修改部分我在上述代码中添加了注释）

训练之后，损失值图将显示误差随着迭代次数的增加而减少。

以绘图的方式对决策边界进行可视化处理，代码如下：

import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_iris


class logistic(object):
    def __init__(self):
        self.W = None

    def train(self, X, y, lr=0.01, num_iters=5000):
        num_train, num_feature = X.shape
        self.W = 0.001 * np.random.randn(num_feature, 1).reshape((-1, 1))
        loss = []
        for i in range(num_iters):
            error, dW = self.compute_loss(X, y)
            self.W += - lr * dW
            loss.append(error)
            if i % 200 == 0:
                print('i= %d, error= %f' % (i, error))
        return loss

    def compute_loss(self, X, y):
        num_train = X.shape[0]
        h = self.output(X)
        loss = - np.sum((y * np.log(h) + (1 - y) * np.log(1 - h))) / num_train
        dW = X.T.dot(h - y) / num_train
        return loss, dW

    def output(self, X):
        g = np.dot(X, self.W)
        return self.sigmoid(g)

    def sigmoid(self, X):
        return 1 / (1 + np.exp(-X))

    def predict(self, X_test):
        h = self.output(X_test)
        return np.where(h >= 0.5, 1, 0)


iris = load_iris()
data = iris.data
target = iris.target
X = data[0:100, [0, 2]]
y = target[0:100]
y = y.reshape((-1, 1))
one = np.ones((X.shape[0], 1))
X_train = np.hstack((one, X))
classify = logistic()
loss = classify.train(X_train, y)
label = np.array(y)
index_0 = np.where(label == 0)
plt.scatter(X[index_0, 0], X[index_0, 1], marker='x', c='b', label='0', s=15)
index_1 = np.where(label == 1)
plt.scatter(X[index_1, 0], X[index_1, 1], marker='o', c='r', label='1', s=15)
# 绘制分类边界线
x1 = np.arange(4, 7.5, 0.5)
x2 = (- classify.W[0] - classify.W[1] * x1) / classify.W[2]
plt.plot(x1, x2, color='black')
plt.xlabel('X1')
plt.ylabel('X2')
plt.legend(loc='upper left')
plt.show()

运行结果如图所示，可以看出，最后学习得到的决策边界成功隔开两个类别。

3.2.4 实验

1.实验目的

理解逻辑回归的算法原理，掌握损失函数的优化方法，并分别利用sklearn的相关包、python语言编程来实现该算法。

2.实验数据

数据集选用经典的鸢尾花数据集。

3.实验要求

实现数据可视化：通过数据集文件导入数据，并使用Matplotlib工具建立对应散点图。

将线性回归参数初始化为0，计算损失函数（cost function）的初始值，根据算法基本原理中的损失函数计算公式来计算。

选择以下两种优化方法分别求解逻辑回归参数。

梯度下降法

牛顿迭代法

对验证集进行验证。

画出分类边界。

4.实验代码

1. 导入鸢尾花数据集，绘制散点图。代码如下：

import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
import pandas as pd

# 加载鸢尾花数据集
iris = load_iris()
iris_df = pd.DataFrame(data=iris.data, columns=iris.feature_names)
iris_df['species'] = iris.target

# 将数字标签转换为类别标签
iris_df['species'] = iris_df['species'].map({0: 'setosa', 1: 'versicolor', 2: 'virginica'})

# 绘制散点图
sns.scatterplot(data=iris_df, x='sepal length (cm)', y='sepal width (cm)', hue='species')
plt.title('Iris Dataset Scatter Plot')
plt.xlabel('Sepal Length (cm)')
plt.ylabel('Sepal Width (cm)')
plt.legend(title='Species')
plt.show()

代码结果如下：

2.线性回归参数初始化为0，根据算法基本原理中的损失函数计算损失函数的初始值。在原理部分提到的损失函数可以进一步变成均方误差（MSE）的形式

$L(f)=\frac{1}{N}\sum_{i=1}^{N}(\hat{y}_i-y_i)^2=\frac{1}{N}\sum_{i=1}^{N}(w\cdot x_i+b-y_i)^2$

代码如下（数据集使用糖尿病病人数据集）：

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split

# 加载糖尿病数据集
diabetes = load_diabetes()
X = diabetes.data  # 特征数据
y = diabetes.target  # 目标数据

# 将y调整成列向量
y = y.reshape(-1, 1)

# 将数据集拆分为训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=0)

# 初始化权重 W 和偏置 b
W = np.zeros((X_train.shape[1], 1))  # 初始化为0，X_train.shape[1]是特征数
b = 0

# 初始预测值 h(x) 全为 0
y_pred = np.dot(X_train, W) + b  # 由于 W 和 b 都是 0，所以预测值为 0

# 计算初始损失值 (MSE)
initial_loss = (1 / (2 * len(y_train))) * np.sum((y_pred - y_train) ** 2)

print("糖尿病数据集的初始损失值:", initial_loss)

代码结果如图：

3.使用梯度下降法和牛顿迭代法分别求解逻辑回归参数。

# 加载鸢尾花数据集
iris = load_iris()
X = iris.data
y = iris.target

# 只保留 Setosa 和 Versicolor 两类 (target == 0 是 Setosa, target == 1 是 Versicolor)
mask = y < 2
X = X[mask]
y = y[mask]

# 目标变量转换为 0 和 1
y = (y == 0).astype(int)

# 标准化数据
scaler = StandardScaler()
X = scaler.fit_transform(X)


# 定义 Sigmoid 函数
def sigmoid(z):
    return 1 / (1 + np.exp(-z))


# 定义梯度计算函数
def compute_gradient(X, y, theta):
    m = X.shape[0]
    h = sigmoid(X @ theta)
    gradient = (X.T @ (h - y)) / m
    return gradient


# 梯度下降法实现
def gradient_descent(X, y, theta, learning_rate, num_iterations):
    for i in range(num_iterations):
        gradient = compute_gradient(X, y, theta)
        theta = theta - learning_rate * gradient
    return theta


# 定义 Hessian 计算函数
def compute_hessian(X, theta):
    m = X.shape[0]
    h = sigmoid(X @ theta)
    R = np.diag(h * (1 - h))
    H = X.T @ R @ X / m
    return H


# 牛顿法实现
def newton_method(X, y, theta, num_iterations):
    for i in range(num_iterations):
        gradient = compute_gradient(X, y, theta)
        hessian = compute_hessian(X, theta)
        theta = theta - np.linalg.inv(hessian) @ gradient
    return theta

设置预测函数进行交叉验证。

# 预测函数
def predict(X, theta):
    return (sigmoid(X @ theta) >= 0.5).astype(int)


# 交叉验证
def cross_validate_model(X, y, method, num_folds=5, learning_rate=0.01, num_iterations=1000):
    kf = KFold(n_splits=num_folds)
    accuracies = []

    for train_index, val_index in kf.split(X):
        X_train, X_val = X[train_index], X[val_index]
        y_train, y_val = y[train_index], y[val_index]

        theta = np.zeros(X_train.shape[1])

        # 根据选择的方法优化参数
        if method == 'gd':
            theta_optimal = gradient_descent(X_train, y_train, theta, learning_rate, num_iterations)
        elif method == 'newton':
            theta_optimal = newton_method(X_train, y_train, theta, num_iterations=10)

        # 在验证集上进行预测
        y_pred = predict(X_val, theta_optimal)

        # 计算准确率
        accuracy = np.mean(y_pred == y_val)
        accuracies.append(accuracy)

    avg_accuracy = np.mean(accuracies)
    return avg_accuracy

而后绘制分类边界。

# 绘制分类边界
def plot_decision_boundary(X, y, theta, title):
    # 创建网格范围
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01),
                         np.arange(y_min, y_max, 0.01))

    # 计算网格点的预测值
    Z = sigmoid(np.c_[xx.ravel(), yy.ravel()] @ theta)
    Z = Z.reshape(xx.shape)

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 设置字体
    plt.rcParams['axes.unicode_minus'] = False

    # 绘制决策边界
    plt.contourf(xx, yy, Z, levels=[0, 0.5, 1], alpha=0.5, cmap='coolwarm')
    plt.scatter(X[:, 0], X[:, 1], c=y, edgecolor='k', cmap='coolwarm')
    plt.xlabel('花萼长度')
    plt.ylabel('花萼宽度')
    plt.title(title)
    plt.show()

主函数

# 主程序l
if __name__ == "__main__":
    # 使用交叉验证对梯度下降和牛顿法进行验证
    accuracy_gd = cross_validate_model(X, y, method='gd', num_folds=5)
    accuracy_newton = cross_validate_model(X, y, method='newton', num_folds=5)

    print(f'梯度下降法在交叉验证中的平均准确率: {accuracy_gd:.9f}')
    print(f'牛顿迭代法在交叉验证中的平均准确率: {accuracy_newton:.9f}')

    # 选择前两个特征进行可视化
    X_selected = X[:, :2]  # 取前两个特征
    X_train, X_test, y_train, y_test = train_test_split(X_selected, y, test_size=0.3, random_state=42)

    # 使用梯度下降法优化
    theta = np.zeros(X_train.shape[1])
    theta_optimal_gd = gradient_descent(X_train, y_train, theta, learning_rate=0.01, num_iterations=1000)

    # 使用牛顿法优化
    theta_optimal_newton = newton_method(X_train, y_train, theta, num_iterations=10)

    # 绘制分类边界
    plot_decision_boundary(X_train, y_train, theta_optimal_gd, '梯度下降法分类边界')
    plot_decision_boundary(X_train, y_train, theta_optimal_newton, '牛顿迭代法分类边界')

结果：

分类结果非常好，没有错误。由于数据较少，并且只选择了花萼宽度和花萼长度作为二分类的类别，对于大规模的数据或许达不到这么好的效果。如果代码存在任何问题，欢迎大家指出。需要源码请访问git仓库：机器学习pytho实战源码

机器学习Python实战-第三章-分类问题-2.逻辑回归算法

3.2.1 原理简介

1.普通线性回归

2.广义线性模型

3.逻辑回归

3.2.2 算法步骤

3.2.3 实战

1.数据集

2.sklearn实现

3.算法实现

3.2.4 实验

1.实验目的

2.实验数据

3.实验要求

4.实验代码

网站公告

今日签到

热门文章

最新发布