案例 2：乳腺肿块特征的逻辑回归#
学习目标与记号#
本案例使用由细针穿刺图像计算得到的 30 个数值特征，建立逻辑回归模型区分数据标签 M 与 B。重点是理解二分类概率模型、损失函数和参数更新过程，而不是给出任何临床建议。
本案例用 \(m\) 表示样本量，用 \(n\) 表示特征数量；设计矩阵记为 \(\boldsymbol{X}\in\mathbb{R}^{m\times n}\)，权重向量记为 \(\boldsymbol{w}\in\mathbb{R}^{n\times1}\)，偏置记为 \(b\)，观测标签向量记为 \(\boldsymbol{y}\in\{0,1\}^{m\times1}\)。主要学习目标如下：
从线性运算结果 \(\boldsymbol{z}=\boldsymbol{X}\cdot\boldsymbol{w}+b\) 得到 sigmoid 函数输出的概率，并写出数值稳定的二元交叉熵；
推导批量梯度，理解为什么逻辑回归的优化目标不是分类准确率；
比较梯度下降与牛顿法每一步的代价和收敛速度；
把概率、阈值与混淆矩阵分开，避免把 0.5 当成自然法则；
说明本案例中医疗数据的适用边界、类别含义和两类误判可能带来的不同后果。
对应正文： “线性回归与逻辑回归”中的逻辑回归、交叉熵、梯度下降和 Newton-Raphson 算法。
数据来源与伦理边界#
说明页：https://archive.ics.uci.edu/dataset/17/breast-cancer-wisconsin-diagnostic
直链：https://archive.ics.uci.edu/static/public/17/breast+cancer+wisconsin+diagnostic.zip
许可：CC BY 4.0
规模：569 条记录、30 个连续特征，标签 M/B。
数据适合算法教学，但样本来源、采集设备和群体代表性有限。本模型未经过临床验证，不能用于诊断、分诊或治疗。为了避免措辞造成误导，下文只报告“数据标签 M/B”。
# 统一环境、随机种子与下载缓存
from pathlib import Path
import hashlib
import os
import shutil
import tempfile
import urllib.request
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
SEED = 42
rng = np.random.default_rng(SEED)
CACHE_ROOT = Path(os.environ.get(
"AI_COURSE_DATA_DIR",
Path.home() / ".cache" / "ai-course-cases",
)).expanduser()
CACHE_ROOT.mkdir(parents=True, exist_ok=True)
def download(url, filename, sha256=None):
"""匿名下载到独立缓存；已有文件仍会计算摘要，避免在没有提示的情况下使用错误数据。"""
destination = CACHE_ROOT / filename
if not destination.exists():
request = urllib.request.Request(url, headers={"User-Agent": "ai-course-case/1.0"})
with urllib.request.urlopen(request, timeout=120) as response:
with tempfile.NamedTemporaryFile(dir=CACHE_ROOT, delete=False) as tmp:
shutil.copyfileobj(response, tmp)
temporary = Path(tmp.name)
temporary.replace(destination)
digest = hashlib.sha256(destination.read_bytes()).hexdigest()
if sha256 is not None and digest != sha256:
destination.unlink(missing_ok=True)
raise ValueError(f"SHA-256 不匹配：{digest}")
print(f"{destination.name}：{destination.stat().st_size / 1024**2:.2f} 兆字节，SHA-256={digest}")
return destination
plt.rcParams["figure.dpi"] = 120
print("缓存目录：", CACHE_ROOT)
缓存目录： /private/tmp/ai-course-case-data
第 1 步：读取原始格式并核对标签#
原始 wdbc.data 没有表头：第 1 列是标识符，第 2 列是标签，之后依次是 10 种细胞核测量的 mean、standard error、worst，共 30 列。标识符不是生物学特征，应立即删除。
import io
import zipfile
URL = "https://archive.ics.uci.edu/static/public/17/breast+cancer+wisconsin+diagnostic.zip"
archive = download(URL, "uci_wdbc.zip")
base = ["radius", "texture", "perimeter", "area", "smoothness",
"compactness", "concavity", "concave_points", "symmetry", "fractal_dimension"]
base_labels = ["半径", "纹理", "周长", "面积", "平滑度",
"紧凑度", "凹度", "凹点数", "对称性", "分形维数"]
group_labels = {"mean": "均值", "se": "标准误", "worst": "最大（最差）"}
feature_names = [f"{name}_{group}" for group in ("mean", "se", "worst") for name in base]
feature_labels = {
f"{name}_{group}": f"{label}（{group_labels[group]}）"
for group in ("mean", "se", "worst")
for name, label in zip(base, base_labels)
}
columns = ["id", "label", *feature_names]
with zipfile.ZipFile(archive) as zf:
candidates = [n for n in zf.namelist() if n.lower().endswith("wdbc.data")]
if len(candidates) != 1:
raise RuntimeError(candidates)
raw = pd.read_csv(io.BytesIO(zf.read(candidates[0])), header=None, names=columns)
print("数据维度：", raw.shape, "；标签计数：", raw["label"].value_counts().to_dict())
assert raw.shape[1] == 32
assert set(raw["label"]) == {"M", "B"}
display(raw.head(3).rename(columns={"id": "标识符", "label": "标签", **feature_labels}))
uci_wdbc.zip：0.05 兆字节，SHA-256=bc154869ef13f753f9e2b5a17e248cfe1ba4b6721db7c4da9f4880e40b05d3af
数据维度： (569, 32) ；标签计数： {'B': 357, 'M': 212}
| 标识符 | 标签 | 半径（均值） | 纹理（均值） | 周长（均值） | 面积（均值） | 平滑度（均值） | 紧凑度（均值） | 凹度（均值） | 凹点数（均值） | ... | 半径（最大（最差）） | 纹理（最大（最差）） | 周长（最大（最差）） | 面积（最大（最差）） | 平滑度（最大（最差）） | 紧凑度（最大（最差）） | 凹度（最大（最差）） | 凹点数（最大（最差）） | 对称性（最大（最差）） | 分形维数（最大（最差）） | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 842302 | M | 17.99 | 10.38 | 122.8 | 1001.0 | 0.11840 | 0.27760 | 0.3001 | 0.14710 | ... | 25.38 | 17.33 | 184.6 | 2019.0 | 0.1622 | 0.6656 | 0.7119 | 0.2654 | 0.4601 | 0.11890 |
| 1 | 842517 | M | 20.57 | 17.77 | 132.9 | 1326.0 | 0.08474 | 0.07864 | 0.0869 | 0.07017 | ... | 24.99 | 23.41 | 158.8 | 1956.0 | 0.1238 | 0.1866 | 0.2416 | 0.1860 | 0.2750 | 0.08902 |
| 2 | 84300903 | M | 19.69 | 21.25 | 130.0 | 1203.0 | 0.10960 | 0.15990 | 0.1974 | 0.12790 | ... | 23.57 | 25.53 | 152.5 | 1709.0 | 0.1444 | 0.4245 | 0.4504 | 0.2430 | 0.3613 | 0.08758 |
3 rows × 32 columns
第 2 步：检查数据#
先查看类别比例和各个特征的数值尺度。本案例不根据全部数据删除“异常值”：极端值可能正是重要样本，而且任何删除规则都只能根据训练数据确定。许多特征之间高度相关，这可能使 Hessian 矩阵对应的线性方程难以稳定求解，因此后面的 Newton-Raphson 更新会加入一个很小的正数以改善数值稳定性。
summary = raw[feature_names].describe().T[["mean", "std", "min", "max"]]
summary = summary.rename(index=feature_labels, columns={
"mean": "均值", "std": "标准差", "min": "最小值", "max": "最大值",
})
display(summary.head(10))
raw["label"].value_counts(normalize=True).rename("类别占比").plot.bar()
plt.ylabel("比例")
plt.title("标签分布")
plt.show()
print("缺失值总数：", int(raw.isna().sum().sum()))
| 均值 | 标准差 | 最小值 | 最大值 | |
|---|---|---|---|---|
| 半径（均值） | 14.127292 | 3.524049 | 6.98100 | 28.11000 |
| 纹理（均值） | 19.289649 | 4.301036 | 9.71000 | 39.28000 |
| 周长（均值） | 91.969033 | 24.298981 | 43.79000 | 188.50000 |
| 面积（均值） | 654.889104 | 351.914129 | 143.50000 | 2501.00000 |
| 平滑度（均值） | 0.096360 | 0.014064 | 0.05263 | 0.16340 |
| 紧凑度（均值） | 0.104341 | 0.052813 | 0.01938 | 0.34540 |
| 凹度（均值） | 0.088799 | 0.079720 | 0.00000 | 0.42680 |
| 凹点数（均值） | 0.048919 | 0.038803 | 0.00000 | 0.20120 |
| 对称性（均值） | 0.181162 | 0.027414 | 0.10600 | 0.30400 |
| 分形维数（均值） | 0.062798 | 0.007060 | 0.04996 | 0.09744 |
缺失值总数： 0
第 3 步：分层划分与训练集标准化#
\(y=1\) 表示 M，\(y=0\) 表示 B。分层划分使训练集与测试集中的类别比例大致相同。均值和标准差只根据训练集计算，测试集只使用已经确定的均值和标准差进行转换。观测标签 \(\boldsymbol{y}\) 保留为 \((m,1)\) 的二维列向量，便于核对梯度的维度。
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
X = raw[feature_names].to_numpy(np.float64)
y = (raw[["label"]].to_numpy() == "M").astype(np.float64)
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=.25, stratify=y, random_state=SEED,
)
scaler = StandardScaler().fit(X_train)
X_train_z = scaler.transform(X_train)
X_test_z = scaler.transform(X_test)
assert X_train_z.shape[0] == y_train.shape[0]
assert y_train.shape[1] == 1
print("训练特征维度：", X_train_z.shape, "；测试特征维度：", X_test_z.shape)
print("训练集和测试集的 M 标签占比：", y_train.mean(), y_test.mean())
训练特征维度： (426, 30) ；测试特征维度： (143, 30)
训练集和测试集的 M 标签占比： 0.3732394366197183 0.3706293706293706
第 4 步：数值稳定的 sigmoid 函数与交叉熵#
模型先计算线性运算结果 \(\boldsymbol{z}=\boldsymbol{X}\cdot\boldsymbol{w}+b\)，再令 \(\boldsymbol{p}=\sigma(\boldsymbol{z})\)。当 \(z\) 的绝对值很大时，直接计算 \(\log(p)\) 和 \(\log(1-p)\) 可能出现 \(\log(0)\)。利用恒等式，单样本损失可以写为 \(\log(1+\exp(z))-yz\)。NumPy 的 logaddexp(0, z) 能以更加稳定的方式完成其中的对数与指数计算。
批量平均梯度为
代码同时返回各个结果，便于在训练前用断言检查其维度。
def sigmoid(z):
out = np.empty_like(z, dtype=np.float64)
positive = z >= 0
out[positive] = 1.0 / (1.0 + np.exp(-z[positive]))
exp_z = np.exp(z[~positive])
out[~positive] = exp_z / (1.0 + exp_z)
return out
def logistic_loss_grad(X_batch, y_batch, w, b, l2=0.0):
z = X_batch @ w + b
p = sigmoid(z)
data_loss = np.mean(np.logaddexp(0.0, z) - y_batch * z)
loss = float(data_loss + .5 * l2 * np.sum(w ** 2))
error = p - y_batch
grad_w = X_batch.T @ error / len(X_batch) + l2 * w
grad_b = float(error.mean())
return loss, grad_w, grad_b, p
w0 = np.zeros((X_train_z.shape[1], 1))
loss0, gw0, gb0, p0 = logistic_loss_grad(X_train_z, y_train, w0, 0.0)
print("初始损失：", loss0, "；权重梯度维度：", gw0.shape,
"；偏置梯度类型：", type(gb0), "；概率维度：", p0.shape)
assert np.isclose(loss0, np.log(2))
初始损失： 0.6931471805599453 ；权重梯度维度： (30, 1) ；偏置梯度类型： <class 'float'> ；概率维度： (426, 1)
第 5 步：用中心差分检查一个权重与偏置#
即使梯度公式写错，逻辑回归的损失有时仍可能下降，因此不能只根据损失曲线判断梯度是否正确。数值梯度用于检查解析梯度，而不是用于正式训练。本例使用 float64 数值和较小的 \(\epsilon\)，并通过相对误差比较两种梯度。
epsilon = 1e-6
j = 7
w_plus, w_minus = w0.copy(), w0.copy()
w_plus[j, 0] += epsilon
w_minus[j, 0] -= epsilon
lp = logistic_loss_grad(X_train_z[:64], y_train[:64], w_plus, 0.0)[0]
lm = logistic_loss_grad(X_train_z[:64], y_train[:64], w_minus, 0.0)[0]
analytic = logistic_loss_grad(X_train_z[:64], y_train[:64], w0, 0.0)[1][j, 0]
numeric = (lp - lm) / (2 * epsilon)
relative_error = abs(analytic - numeric) / max(1.0, abs(analytic), abs(numeric))
print({"解析梯度": analytic, "数值梯度": numeric, "相对误差": relative_error})
assert relative_error < 1e-7
{'解析梯度': np.float64(-0.37804498587415436), '数值梯度': -0.37804498587146185, '相对误差': np.float64(2.6925128793209296e-12)}
第 6 步：批量梯度下降#
逻辑回归损失是可微函数，而准确率需要先根据阈值把概率转换为类别，其数值会随着阈值或参数呈阶梯式变化，因此不能直接用普通梯度优化准确率。\(L_2\) 惩罚项只作用于 \(\boldsymbol{w}\)，不作用于偏置 \(b\)。程序每 20 步保存一次训练损失，测试集不参与模型或训练方案的选择。
w_gd = np.zeros((X_train_z.shape[1], 1))
b_gd = 0.0
gd_history = []
for step in range(1200):
loss, grad_w, grad_b, _ = logistic_loss_grad(
X_train_z, y_train, w_gd, b_gd, l2=1e-3,
)
w_gd -= 0.08 * grad_w
b_gd -= 0.08 * grad_b
if step % 20 == 0:
gd_history.append((step, loss))
gd_history = np.asarray(gd_history)
plt.semilogy(gd_history[:, 0], gd_history[:, 1], label="梯度下降")
plt.xlabel("更新步")
plt.ylabel("训练交叉熵 + L2")
plt.legend()
plt.show()
第 7 步：Newton-Raphson 算法为何更新次数少，但每步计算更多#
把权重和偏置合并为参数向量 \(\boldsymbol{\theta}\) 后，Hessian 矩阵为 \(\widetilde{\boldsymbol{X}}^{\mathsf T}\cdot\boldsymbol{R}\cdot\widetilde{\boldsymbol{X}}/m\)，其中 \(\boldsymbol{R}\) 是对角矩阵，其第 \(i\) 个对角元素为 \(p_i(1-p_i)\)。Newton-Raphson 算法先求解 \(\boldsymbol{H}\cdot\boldsymbol{\Delta}=\boldsymbol{g}\)，再更新 \(\boldsymbol{\theta}\leftarrow\boldsymbol{\theta}-\boldsymbol{\Delta}\)。每一步都要构造并求解一个 \((n+1)\) 维线性方程组，因此当特征很多时，一步 Newton-Raphson 更新通常比一步梯度下降需要更多计算。本例只有 \(n=30\) 个特征，适合用来比较两种方法。
代码中加入的很小正数是数值稳定处理，不改变这里要讲解的基本公式；它用于降低特征高度相关时线性方程求解不稳定的风险。
Xb = np.c_[X_train_z, np.ones((len(X_train_z), 1))]
theta = np.zeros((Xb.shape[1], 1))
newton_history = []
for step in range(12):
z = Xb @ theta
p = sigmoid(z)
loss = float(np.mean(np.logaddexp(0.0, z) - y_train * z))
gradient = Xb.T @ (p - y_train) / len(Xb)
curvature = (p * (1 - p)).ravel()
hessian = Xb.T @ (Xb * curvature[:, None]) / len(Xb)
hessian[:-1, :-1] += 1e-3 * np.eye(Xb.shape[1] - 1)
hessian += 1e-8 * np.eye(Xb.shape[1])
step_direction = np.linalg.solve(hessian, gradient)
theta -= step_direction
newton_history.append(loss)
print("牛顿法损失：", np.round(newton_history, 5))
w_newton, b_newton = theta[:-1], float(theta[-1, 0])
牛顿法损失： [0.69315 0.23737 0.13534 0.08711 0.06131 0.04577 0.03649 0.03119 0.028
0.02583 0.02421 0.02293]
第 8 步：概率评价与阈值评价分开#
ROC 曲线下面积（ROC-AUC）衡量模型把正类排在负类之前的能力，对数损失衡量完整概率预测的质量，而混淆矩阵必须在给定分类阈值后才能计算。这里固定阈值 0.5，只用于比较两种优化算法是否得到接近的结果。真实医疗决策必须在独立验证集上结合漏报与误报的实际代价选择阈值，并经过临床和监管验证。
from sklearn.metrics import (
ConfusionMatrixDisplay, accuracy_score, log_loss, roc_auc_score,
)
models = {
"梯度下降": sigmoid(X_test_z @ w_gd + b_gd).ravel(),
"牛顿法": sigmoid(X_test_z @ w_newton + b_newton).ravel(),
}
rows = []
for name, probability in models.items():
prediction = (probability >= .5).astype(int)
rows.append({
"模型": name,
"对数损失": log_loss(y_test.ravel(), probability),
"ROC 曲线下面积": roc_auc_score(y_test.ravel(), probability),
"阈值 0.5 的准确率": accuracy_score(y_test.ravel(), prediction),
})
display(pd.DataFrame(rows))
matrix_display = ConfusionMatrixDisplay.from_predictions(
y_test.ravel(), (models["梯度下降"] >= .5).astype(int),
display_labels=["B", "M"], cmap="Blues",
)
matrix_display.ax_.set_xlabel("预测标签")
matrix_display.ax_.set_ylabel("真实标签")
plt.title("固定阈值 0.5（仅教学比较）")
plt.show()
| 模型 | 对数损失 | ROC 曲线下面积 | 阈值 0.5 的准确率 | |
|---|---|---|---|---|
| 0 | 梯度下降 | 0.066350 | 0.997694 | 0.979021 |
| 1 | 牛顿法 | 0.258992 | 0.980503 | 0.958042 |
第 9 步：阈值改变的不是模型概率，而是行动规则#
扫描阈值并画出召回率与精确率。为了保持测试集纯净，正式项目应在验证集上选阈值；这里的曲线只用于展示机制，不据此宣称某阈值可用于临床。
from sklearn.metrics import precision_score, recall_score
probability = models["梯度下降"]
threshold_rows = []
for threshold in np.linspace(.05, .95, 19):
prediction = (probability >= threshold).astype(int)
threshold_rows.append({
"阈值": threshold,
"精确率": precision_score(y_test.ravel(), prediction, zero_division=0),
"召回率": recall_score(y_test.ravel(), prediction, zero_division=0),
})
threshold_table = pd.DataFrame(threshold_rows)
threshold_table.plot(x="阈值", y=["精确率", "召回率"], ylim=(0, 1.05))
plt.title("阈值—指标关系（测试集上的教学可视化）")
plt.show()
简单基线与结果解释规则#
在训练较复杂的算法前，应先报告“始终预测多数类”这一简单方法的准确率，以及训练集中 M 标签所占的比例。多数类基线没有排序能力，也不能给出有用的个体概率，但它能揭示一个常见问题：类别不平衡时，很高的准确率并不一定说明模型学到了有效规律。逻辑回归至少应在交叉熵、ROC-AUC 以及与任务有关的召回率和精确率上提供更多信息。
阅读结果时先确认标签含义：代码把 M 编码为 1，因此召回率表示真实 M 标签中被模型识别出的比例。接着查看对数损失，判断模型是否给出了过于极端却错误的概率；再查看 AUC 判断排序能力；最后在预先确定的阈值下查看混淆矩阵。若两种优化算法的准确率相同但对数损失不同，说明它们给出的类别可能相同，但概率预测质量不同，不能据此认为两个模型完全等价。
正式研究还应从训练数据中划分验证集。阈值、正则化强度、特征选择和停止轮数只能根据训练集与验证集确定，测试集只在模型与训练方案确定后评价一次。由于样本量较小，还应使用分层交叉验证说明结果的不确定性，不能把单次随机划分中很小的数值差异解释为稳定改进。
进一步诊断：系数不是“重要性真相”#
标准化后的权重绝对值可以作为局部线性关联的线索，但高度相关的半径、周长和面积会相互分摊系数；系数符号也可能因加入相关变量而改变。因此不要按单次拟合的绝对值给特征做生物学排名。更稳妥的教学检查是重复交叉验证，观察系数方向和性能区间，并把这种不稳定性写进结论。
还应检查极端概率。若模型把某样本预测为接近 0 或 1 却判断错误，它对 log loss 的影响会很大；这通常提示分布外样本、异常量纲、标签噪声或模型错设。查看这类样本可以诊断算法，但展示时不要泄露标识符，也不能根据测试集个案反复调整模型。
结论、局限与常见错误#
sigmoid 函数给出模型概率，但概率是否校准还需单独验证；0.5 不是普遍最优阈值。
牛顿法步数少不代表总成本低。它的内存和线性求解成本随特征数快速增加。
切勿把 id 当特征；也不要先在全数据上标准化或按全数据挑阈值。
类别标签来自特定数据与流程，不能外推到新医院、新设备或新人群。
本案例只检查算法实现。任何真实医疗用途都需要前瞻性验证、公平性检查、人工监督和合规审批。
最小检查清单#
每次重跑都应确认：标签编码未反转；训练与测试样本不重叠；标准化器只拟合训练集；损失为有限值；预测概率全部位于 [0,1]；指标函数的正类仍是 M。
综合练习#
新划出验证集，在“漏报代价是误报 5 倍”的假设下选择阈值，再只评价一次测试集。
把 L2 强度改为 0、10^-4、10^-2、1，比较权重范数、log loss 与 AUC。
记录梯度下降和 Newton-Raphson 算法的实际运行时间、更新次数与最终损失。比较时使用相同的数据、停止条件和运行环境，并说明两种算法在单步计算量与所需更新次数上的差别。
在 z=±1000 处比较稳定损失与直接 log(sigmoid(z))，解释溢出来源。
用 scikit-learn 的
LogisticRegression作为经过广泛使用的实现进行对照，并核对正类标签和正则化项的定义。