3. 代码实现

接上一篇文章

3.1 问题背景

应用案例背景

假设我们有一个数据集，包含患者的多种生理指标（如年龄、性别、体重指数（BMI）、血糖水平等）以及他们是否被诊断为糖尿病（是或否）。我们的目标是建立一个逻辑回归模型，根据这些生理指标预测一个患者是否有糖尿病的概率。

数据集表示

假设我们的数据集包含以下特征：

• $X_1$: 年龄
• $X_2$: 性别（0表示女性，1表示男性）
• $X_3$: 体重指数（BMI）
• $X_4$: 血糖水平
• $Y$: 是否有糖尿病（0表示没有，1表示有）

逻辑回归模型

我们将使用逻辑回归模型来预测$P(Y=1|X)$，即给定生理指标X时，患者有糖尿病的概率。模型的形式为：

$P(Y=1|X)=\frac{1}{1+e^{-({\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+{\beta }_4{X}_4)}}$

3.2 数据准备

# 模拟创建数据
import pandas as pd
import numpy as np
# 设置随机种子以保证结果可复现
np.random.seed(0)
# 创建数据
data_size = 1000  # 数据集大小
# 年龄范围从20到80
ages = np.random.randint(20, 80, size=data_size)
# 性别为0或1
genders = np.random.randint(0, 2, size=data_size)
# BMI范围从18到35
bmis = np.random.uniform(18, 35, size=data_size)
# 血糖水平范围从70到140
glucose_levels = np.random.uniform(70, 140, size=data_size)
# 指定结果：是否有糖尿病
labels = np.random.randint(0, 2, size=data_size)

# 定义结果变量——先假定一个模型来生成数据的结果
# 使用逻辑函数生成糖尿病的概率，并依此生成0或1的标签
beta = np.array([-8, 0.05, 0.25, 0.1, 0.05])  # 假定的模型参数
X = np.column_stack((np.ones(data_size), ages, genders, bmis, glucose_levels))
linear_combination = X.dot(beta)
probabilities = 1 / (1 + np.exp(-linear_combination))
labels = np.random.binomial(1, probabilities)# 创建DataFrame
df = pd.DataFrame({'Age': ages,'Gender': genders,'BMI': bmis,'Glucose_Level': glucose_levels,'Diabetes': labels
})# 显示前几行数据
print(df.head())

   Age  Gender        BMI  Glucose_Level  Diabetes
0   64       1  21.150926      74.332242         1
1   67       0  34.928860     127.491243         1
2   73       0  20.199048      96.576651         1
3   20       1  26.014774     110.008514         1
4   23       1  19.157583     138.848879         1

3.3 数据探索

# 数据集中有多少患者被诊断为糖尿病
print(df['Diabetes'].value_counts())

1    880
0    120
Name: Diabetes, dtype: int64

3.4 手动实现逻辑回归

逻辑回归模型的训练过程就是最大化似然函数的过程。我们可以使用梯度下降法来最大化似然函数。下面我们将手动实现逻辑回归模型的训练过程。

使用梯度下降来最大化逻辑回归的似然函数，实际上我们通常是在最小化对数似然函数的负值，这被称为对数损失或交叉熵损失。这种方法称为梯度下降，其目标是找到一组参数（权重和偏置），使得损失函数最小化。
公式为：
$$J(\mathbf{w},b)=-\sum _{i=1}^n\left[y^{(i)}\log (\sigma (z^{(i)}))+(1-y^{(i)})\log (1-\sigma (z^{(i)}))\right]$$

其中${\beta }_0$就是b，w就是${\beta }_1到{\beta }_n$，z就是线性组合${\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+{\beta }_4{X}_4$，$\sigma$就是sigmoid函数。
现在我们的目标就是最小化上面的公式，要求它的最小值对应的$\beta$，我们就又需要用到梯度下降法了。
在梯度下降法中，我们更新权重 $\mathbf{w}$ 和偏置 $b$ 以最小化损失 $J$，更新规则为：

$\mathbf{w}\leftarrow \mathbf{w}-\alpha {\nabla }_{\mathbf{w}}J$

$b\leftarrow b-\alpha \frac{\partial J}{\partial b}$

所以现在我们就先要对上面的公式求导，求出梯度，然后再用梯度下降法来更新$\beta$。

对 ${\beta }_0$ 的偏导数（对于偏置项）

1. 写出损失函数 $J$：
   $J(w,b)=-{\sum }_{i=1}^n\left[y^{(i)}\log (\sigma (z^{(i)}))+(1-y^{(i)})\log (1-\sigma (z^{(i)}))\right]$
2. 写出 $z^{(i)}$：
   $z^{(i)}={\beta }_0+{\beta }_1{x}_1^{(i)}+\dots +{\beta }_n{x}_n^{(i)}$
3. 对 $z^{(i)}$ 求 ${\beta }_0$ 的偏导数：
   $\frac{\partial z^{(i)}}{\partial {\beta }_0}=1$
4. 对 Sigmoid 函数求导：
   $\frac{d\sigma }{dz}=\sigma (z)\cdot (1-\sigma (z))$
   其中，$\sigma (z)=\frac{1}{1+e^{-z}}$。
5. 对损失函数 $J$ 关于 $\sigma (z^{(i)})$ 求偏导数：
   $\frac{\partial J}{\partial \sigma (z^{(i)})}=-\left(\frac{y^{(i)}}{\sigma (z^{(i)})}-\frac{1-y^{(i)}}{1-\sigma (z^{(i)})}\right)$
6. 应用链式法则得到 $J$ 关于 ${\beta }_0$ 的偏导数：
   $\frac{\partial J}{\partial {\beta }_0}={\sum }_{i=1}^n\frac{\partial J}{\partial \sigma (z^{(i)})}\cdot \frac{d\sigma }{d{z}^{(i)}}\cdot \frac{\partial z^{(i)}}{\partial {\beta }_0}$
7. 计算最终的偏导数表达式：
   $\frac{\partial J}{\partial {\beta }_0}={\sum }_{i=1}^n\left(\sigma (z^{(i)})-y^{(i)}\right)$

对 ${\beta }_j$ 的偏导数（对于权重项）

1. 对 $z^{(i)}$ 求 ${\beta }_j$ 的偏导数：
   $\frac{\partial z^{(i)}}{\partial {\beta }_j}=x_j^{(i)}$

2. 使用链式法则得到 $J$ 关于 ${\beta }_j$ 的偏导数：
   $\frac{\partial J}{\partial {\beta }_j}={\sum }_{i=1}^n\frac{\partial J}{\partial \sigma (z^{(i)})}\cdot \frac{d\sigma }{d{z}^{(i)}}\cdot \frac{\partial z^{(i)}}{\partial {\beta }_j}$

3. 计算最终的偏导数表达式：
   $\frac{\partial J}{\partial {\beta }_j}={\sum }_{i=1}^n\left(\sigma (z^{(i)})-y^{(i)}\right)x_j^{(i)}$

   在实际的算法实现中，这些导数会用来更新 ${\beta }_0$ 和 ${\beta }_j$：

因为$\sigma (z^{(i)})$其实就是我们的预测值，所以这个公式的意思就是我们的预测值减去实际值，然后再求和，就是我们的偏导数。

${\beta }_0:={\beta }_0-\alpha \frac{\partial J}{\partial {\beta }_0}$
${\beta }_j:={\beta }_j-\alpha \frac{\partial J}{\partial {\beta }_j}$

1. 定义 Sigmoid 函数： 这是逻辑回归中的核心函数，将线性组合映射到 (0, 1) 区间。
2. 定义损失函数： 这个函数计算当前权重和偏置下的对数损失。
3. 定义梯度计算函数： 计算损失函数关于每个参数的梯度。
4. 定义梯度下降更新规则： 使用计算得到的梯度更新权重和偏置。
5. 执行训练过程： 迭代执行梯度下降步骤，直至满足收敛条件或达到指定的迭代次数。

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

定义损失函数
公式为：J(w, b) = $-{\sum }_{i=1}^n[y^{(i)}\log (\sigma (z^{(i)})+(1-y^{(i)})\log (1-\sigma (z^{(i))})]$

# 定义对数损失函数——我们的目的就是让损失最小，就是让似然函数最大，当求的这些beta能够使得损失最小的时候，就是我们的最优解，beta就是我们的最优参数
def compute_loss(X, y, w, b):# 对应上面的z=(β0+β1X1+β2X2+β3X3+β4X4)z = np.dot(X, w) + b probs = sigmoid(z)# np.mean()求损失的均值loss = -np.mean(y * np.log(probs) + (1 - y) * np.log(1 - probs))return loss

在 compute_gradients 函数中，参数的含义如下：

X: 特征数据矩阵。它是一个二维数组，其中每一行代表一个训练样本，每一列代表一个特征。
y: 标签向量。它是一个一维数组，包含与特征矩阵 X 中的每一行相对应的标签。在二分类逻辑回归中，这些标签通常是 0 或 1。
w: 权重向量。这是一个一维数组，其中每个元素对应于 X 矩阵中某个特征的权重。
b: 偏置项。这是一个标量，它与权重向量 w 一起决定了模型的预测。

# 定义梯度计算函数——这个函数就是对上面的损失函数求导，求出梯度
# w对应beta1到beta4，b对应beta0，X对应X1到X4，y对应Y（实际的值）
def compute_gradients(X, y, w, b):z = np.dot(X, w) + bprobs = sigmoid(z)error = probs - y # 对应sigmoid(z)-y# 梯度计算，就是前面推导的公式grad_w = np.dot(X.T, error) / len(y) # 对应(1/n)*sum((sigmoid(z)-y)*X)grad_b = np.mean(error) # 对应(1/n)*sum(sigmoid(z)-y)return grad_w, grad_b

# 定义梯度下降函数
def gradient_descent(X, y, w, b, alpha, iterations):''':param X: 表示特征数据矩阵:param y: 表示标签向量（实际值）:param w: 权重向量:param b: 偏置项:param alpha: 学习率——就是我们的步长，每一次更新多少，是一个自己定义的参数:param iterations: 迭代次数——就是我们的梯度下降要迭代多少次:return: 返回最终的权重和偏置'''for i in range(iterations):grad_w, grad_b = compute_gradients(X, y, w, b)# 更新权重和偏置w -= alpha * grad_wb -= alpha * grad_b# 每迭代一定次数打印损失if i % 100 == 0:print(f"Iteration {i}: Loss = {compute_loss(X, y, w, b)}")return w, b

准备数据

• $X_1$: 年龄
• $X_2$: 性别（0表示女性，1表示男性）
• $X_3$: 体重指数（BMI）
• $X_4$: 血糖水平
• $Y$: 是否有糖尿病（0表示没有，1表示有）

# 准备数据
# x = df的前四列，y = df的最后一列
X = df.iloc[:, :-1].values # iloc是通过行号来取行数据的，这里取所有行，去掉最后一列
y = df.iloc[:, -1].values # 取所有行，取最后一列
X.shape, y.shape

((1000, 4), (1000,))

# 将X和y分为训练集和测试集
from sklearn.model_selection import train_test_split
# 这里就是把X和y分为训练集和测试集，test_size=0.2表示测试集占20%，random_state=0表示随机种子，保证每次运行结果一样
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

# 特征标准化——目的是让数据的均值为0，方差为1，这样可以加快梯度下降的收敛速度
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

# 初始化权重和偏置
# 这里的w是一个一维数组，长度为X_train的列数，也就是特征的个数，这里是4
w = np.zeros(X_train.shape[1])
b = 0

# 训练模型
# 这里的alpha是学习率，iterations是迭代次数
w, b = gradient_descent(X_train, y_train, w, b, alpha=0.1, iterations=10000)

Iteration 0: Loss = 0.67694848293987
Iteration 100: Loss = 0.31170512579698983
Iteration 200: Loss = 0.28006450564964497
Iteration 300: Loss = 0.26967717122142904
Iteration 400: Loss = 0.2649489155820069
Iteration 500: Loss = 0.26246507021700166
Iteration 600: Loss = 0.2610482148241369
Iteration 700: Loss = 0.26019578953398675
.........................................
Iteration 8300: Loss = 0.2586331617938211
Iteration 8400: Loss = 0.25863316179381984
Iteration 8500: Loss = 0.2586331617938189
Iteration 8600: Loss = 0.25863316179381834
Iteration 8700: Loss = 0.25863316179381785
Iteration 8800: Loss = 0.2586331617938175
Iteration 8900: Loss = 0.2586331617938173
Iteration 9000: Loss = 0.2586331617938171
Iteration 9100: Loss = 0.25863316179381696
Iteration 9200: Loss = 0.25863316179381696
Iteration 9300: Loss = 0.2586331617938169
Iteration 9400: Loss = 0.25863316179381685
Iteration 9500: Loss = 0.25863316179381685
Iteration 9600: Loss = 0.2586331617938168
Iteration 9700: Loss = 0.2586331617938168
Iteration 9800: Loss = 0.2586331617938168
Iteration 9900: Loss = 0.2586331617938168

# 训练好模型后，我们可以使用它来进行预测
# 预测概率
test_probs = sigmoid(np.dot(X_test, w) + b)
# 预测类别
test_preds = np.where(test_probs > 0.5, 1, 0)
test_preds

array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,1, 1])

# 计算准确率
from sklearn.metrics import accuracy_score
accuracy = accuracy_score(y_test, test_preds)
accuracy

0.885

# 比较真实值和预测值
results = pd.DataFrame({'Actual': y_test,'Predicted': test_preds,'Probability': test_probs
})
results.head(100)

	Actual	Predicted	Probability
0	1	1	0.974574
1	1	1	0.961203
2	1	1	0.874061
3	1	1	0.996859
4	1	1	0.832528
...	...	...	...
95	0	1	0.873257
96	1	1	0.944671
97	1	1	0.903943
98	1	1	0.979292
99	1	1	0.998257

100 rows × 3 columns

# 最终的权重和偏置
print(f"Weight: {w}")
print(f"Bias: {b}")
# 最终准确度
print(f"Final Accuracy: {accuracy}")

Weight: [1.41949725 0.42282911 0.70729532 0.90150334]
Bias: 3.0447132565927673
Final Accuracy: 0.885

总结

在这个项目中，我们手动实现了逻辑回归模型的训练过程。我们首先定义了 Sigmoid 函数和对数损失函数，然后使用梯度下降法最小化损失函数。我们将这个模型应用于一个模拟数据集，其中包含患者的生理指标和是否患有糖尿病的标签。最终，我们评估了模型的准确性，并展示了预测结果。
通过最终的结果，我们可以看到我们的模型在测试集上的准确率为 0.885，说明我们的逻辑回归模型在这个数据集上表现还是不错的。

补充——采用sklearn库实现逻辑回归

# 使用sklearn库实现逻辑回归
from sklearn.linear_model import LogisticRegression
# 创建逻辑回归模型
model = LogisticRegression()
# 训练模型
model.fit(X_train, y_train)
# 预测
y_pred = model.predict(X_test)
# 计算准确率
accuracy = accuracy_score(y_test, y_pred)
accuracy

0.885

# 比较调库实现的逻辑回归和手动实现的逻辑回归的权重和偏置
print("Manual Model:")
print(f"Weight: {w}")
print(f"Bias: {b}")
print("\nSklearn Model:")
print(f"Weight: {model.coef_}")
print(f"Bias: {model.intercept_}")

Manual Model:
Weight: [1.41949725 0.42282911 0.70729532 0.90150334]
Bias: 3.0447132565927673Sklearn Model:
Weight: [[1.37051336 0.40978953 0.68512706 0.87338347]]
Bias: [2.988773]

可以发现，两者的权重和偏置非常接近，这说明我们手动实现的逻辑回归模型是正确的。同时，我们也可以看到，使用sklearn库实现的逻辑回归模型的准确率与我们手动实现的模型准确率相同，这说明我们的模型在这个数据集上表现良好。