ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

西瓜书线性模型Python手实现:对率回归与LDA从推导到可视化

西瓜书线性模型Python手实现:对率回归与LDA从推导到可视化 简介本资源是《机器学习》周志华著俗称“西瓜书”第三章“线性模型”的配套Python实践代码包面向机器学习初学者与高校课程实践者聚焦对率回归与线性判别分析两大核心算法的动手实现与验证。资源完整覆盖教材3.3–3.5节全部编程任务包括在西瓜数据集3.0α上实现并评估对率回归与LDA以及在Diabetes、breast_cancer两个UCI数据集上对比10折交叉验证与留一法的泛化误差估计效果。压缩包共8个文件5个.py主程序脚本、1个.xls、1个.csv、1个.txt总大小仅78KB轻量紧凑便于快速运行与调试其中logistics.py、3.5.py等模块结构清晰注释详实辅以watermelon_3a.txt和breast_cancer.csv等标准数据文件开箱即用。目前已有3195人学习下载是理解线性分类模型原理、掌握sklearn底层实现逻辑及交叉验证实践方法的优质入门级代码参考。1. 西瓜书第三章线性模型Python实现不是抄公式而是把周志华书里那张手绘图跑通的实操包你翻过《机器学习》西瓜书第三章看到“对率回归”“线性判别分析”“最小二乘法”这些词时是不是一边划重点一边怀疑这真能用几行Python跑出来不是调sklearn.fit()就完事而是从零推导梯度、手动写Sigmoid、在watermelon_3a.txt上画出决策边界——这才是本资源的真实定位。它不是课后习题答案集而是一套可调试、可打断点、可改参数、可对比理论推导与数值结果偏差的工程化代码包。里面包含原始西瓜数据集3.0αwatermelon_3a.txt、UCI经典数据集Diabetes.xls、breast_cancer.csv、5个核心脚本3.4.1.py到3.5.py全部基于纯NumPyMatplotlib实现不依赖任何黑盒封装。适合刚学完第三章推导、想验证自己是否真懂“为什么损失函数要取对数似然”“为什么LDA投影方向是Sw⁻¹(Sb)w”的本科生也适合想给新人讲清线性模型底层逻辑的带教工程师。别被“课后答案”误导——这里每个.py文件都留了debug断点位置、每份数据都标注了字段含义、每次plot都标清坐标物理意义。你运行一次就能亲手看见那个课本里抽象的“超平面”怎么在二维散点图上一刀切开好瓜坏瓜。2. 从西瓜数据集3.0α出发手写对率回归全流程拆解2.1 数据加载与预处理为什么watermelon_3a.txt必须手动解析西瓜书第三章使用的watermelon_3a.txt是典型的教学精简数据集共17个样本含密度、含糖率两个特征以及“好瓜/坏瓜”标签。但注意它不是CSV格式而是制表符分隔中文标签无表头。直接用pandas.read_csv会出错必须手动解析import numpy as np def load_watermelon_3a(): data [] labels [] with open(watermelon_3a.txt, r, encodingutf-8) as f: for line in f: if not line.strip(): # 跳过空行 continue parts line.strip().split(\t) # 前两列是连续特征密度、含糖率第三列是中文标签 density float(parts[0]) sugar float(parts[1]) label 1 if parts[2] 是 else 0 # 是→正例(好瓜)否→负例(坏瓜) data.append([density, sugar]) labels.append(label) return np.array(data), np.array(labels) X, y load_watermelon_3a() print(f数据形状: {X.shape}, 标签形状: {y.shape}) # 输出: 数据形状: (17, 2), 标签形状: (17,)提示load_watermelon_3a()函数的关键在于三处硬编码处理——split(\t)应对制表符分隔、parts[2] 是处理中文标签、float()强制转浮点。这是西瓜书配套数据的固有特性不是bug。若你用其他数据集此处需按实际分隔符和标签格式调整。2.2 对率回归建模从数学推导到梯度下降实现对率回归Logistic Regression的本质是求解最大似然估计。西瓜书公式(3.27)给出对数似然函数$$ \ell(\boldsymbol{w}, b) \sum_{i1}^m \left[ y_i (\boldsymbol{w}^\top \boldsymbol{x}_i b) - \ln\left(1 \exp(\boldsymbol{w}^\top \boldsymbol{x}_i b)\right) \right] $$最大化该函数等价于最小化其负值。我们用梯度下降法迭代更新参数def sigmoid(z): # 防止溢出z过大时exp(z)→infz过小时exp(z)→0 z np.clip(z, -500, 500) # 限制z范围避免nan return 1 / (1 np.exp(-z)) def logistic_regression(X, y, lr0.1, max_iter1000, tol1e-6): m, n X.shape w np.zeros(n) # 权重向量 b 0.0 # 偏置项 losses [] for i in range(max_iter): # 计算线性组合 z X w b z X w b # 计算预测概率 p sigmoid(z) p sigmoid(z) # 计算损失负对数似然向量化 loss -np.mean(y * np.log(p 1e-15) (1 - y) * np.log(1 - p 1e-15)) losses.append(loss) # 计算梯度∂L/∂w (1/m) * X.T (p - y), ∂L/∂b (1/m) * sum(p - y) dw (1/m) * X.T (p - y) db (1/m) * np.sum(p - y) # 更新参数 w - lr * dw b - lr * db # 收敛判断梯度模长小于tol if np.linalg.norm(dw) tol and abs(db) tol: print(f梯度下降在第{i1}轮收敛) break return w, b, losses # 执行训练 w, b, losses logistic_regression(X, y, lr0.5, max_iter2000) print(f训练完成权重w{w}, 偏置b{b:.4f})参数说明lr0.5学习率设为0.5而非默认0.1因西瓜数据集样本少、特征尺度小过小学习率导致收敛极慢z np.clip(z, -500, 500)关键防溢出操作sigmoid输入超出[-500,500]时exp计算会溢出这是新手最常翻车点np.log(p 1e-15)防止p0或p1时log(0)报错1e-15是经验性极小值收敛判断用梯度模长而非损失变化量因损失本身震荡大梯度更稳定。2.3 决策边界可视化用等高线画出西瓜书图3.3的复刻版西瓜书图3.3展示的是二维特征空间中的决策边界即p0.5的直线。我们用plt.contour绘制等高线并叠加原始数据点import matplotlib.pyplot as plt def plot_decision_boundary(X, y, w, b): # 创建网格点 h 0.01 x_min, x_max X[:, 0].min() - 0.1, X[:, 0].max() 0.1 y_min, y_max X[:, 1].min() - 0.1, X[:, 1].max() 0.1 xx, yy np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h)) # 计算网格点上的预测概率 Z sigmoid(np.c_[xx.ravel(), yy.ravel()] w b) Z Z.reshape(xx.shape) # 绘制等高线p0.5即决策边界 plt.figure(figsize(8, 6)) plt.contour(xx, yy, Z, levels[0.5], colorsred, linewidths2, linestyles--) plt.contourf(xx, yy, Z, levelsnp.linspace(0, 1, 11), cmapRdYlBu_r, alpha0.6) # 绘制原始数据点 pos X[y 1] neg X[y 0] plt.scatter(pos[:, 0], pos[:, 1], cgreen, markero, s80, label好瓜) plt.scatter(neg[:, 0], neg[:, 1], cred, markerx, s80, label坏瓜) plt.xlabel(密度) plt.ylabel(含糖率) plt.title(西瓜数据集3.0α对率回归决策边界) plt.legend() plt.grid(True, alpha0.3) plt.show() plot_decision_boundary(X, y, w, b)关键细节np.c_[xx.ravel(), yy.ravel()]将网格展平为(N,2)矩阵适配向量化计算contourf填充概率热力图contour(..., levels[0.5])单独画出p0.5的红色虚线——这就是西瓜书图3.3中那条分割线绿色圆圈代表正例好瓜红色叉号代表负例坏瓜直观验证模型是否学到了“高密度高含糖率→好瓜”的物理规律。3. UCI数据集交叉验证实战10折vs留一法的误差对比实验3.1 Diabetes数据集加载与标准化为什么必须做Z-score归一化UCI Diabetes数据集Diabetes.xls包含442个糖尿病患者样本10个生理特征如年龄、BMI、血压等目标变量是疾病进展指标连续值。但注意原书3.4节要求比较“对率回归的错误率”而Diabetes是回归任务这里存在一个关键教学陷阱——实际应将其转化为二分类问题。常见做法是以目标变量中位数为阈值大于中位数记为1病情严重否则为0。import pandas as pd from sklearn.preprocessing import StandardScaler def load_diabetes_binary(): # 加载Excel文件需安装openpyxl df pd.read_excel(Diabetes.xls, headerNone) X df.iloc[:, :-1].values # 前10列为特征 y_continuous df.iloc[:, -1].values # 最后一列为目标 # 转为二分类以中位数为界 threshold np.median(y_continuous) y (y_continuous threshold).astype(int) # Z-score标准化消除特征量纲差异避免梯度爆炸 scaler StandardScaler() X_scaled scaler.fit_transform(X) return X_scaled, y X_dia, y_dia load_diabetes_binary() print(fDiabetes数据集{X_dia.shape[0]}样本{X_dia.shape[1]}特征正例比例{y_dia.mean():.3f})注意StandardScaler是必须步骤。Diabetes各特征量纲差异极大年龄≈几十BMI≈二十多血压≈百位不标准化会导致梯度下降方向严重偏移10折CV结果完全不可信。3.2 10折交叉验证实现手写K-Fold而非调用sklearn为彻底理解CV原理我们手动实现10折划分非随机打乱保持原始顺序def k_fold_split(X, y, k10): 手动实现k折划分返回k组(train_idx, test_idx) n len(X) fold_size n // k indices np.arange(n) folds [] for i in range(k): start i * fold_size end start fold_size if i k-1 else n test_idx indices[start:end] train_idx np.concatenate([indices[:start], indices[end:]]) folds.append((train_idx, test_idx)) return folds def evaluate_logistic_cv(X, y, k10, lr0.1, max_iter500): 执行k折CV返回每折的错误率 folds k_fold_split(X, y, k) errors [] for i, (train_idx, test_idx) in enumerate(folds): X_train, y_train X[train_idx], y[train_idx] X_test, y_test X[test_idx], y[test_idx] # 训练 w, b, _ logistic_regression(X_train, y_train, lrlr, max_itermax_iter) # 测试计算错误率 z_test X_test w b p_test sigmoid(z_test) y_pred (p_test 0.5).astype(int) error_rate np.mean(y_pred ! y_test) errors.append(error_rate) print(f第{i1}折训练{len(train_idx)}样本测试{len(test_idx)}样本错误率{error_rate:.4f}) return np.array(errors) # 执行10折CV errors_10fold evaluate_logistic_cv(X_dia, y_dia, k10, lr0.2) print(f\n10折CV平均错误率{errors_10fold.mean():.4f} ± {errors_10fold.std():.4f})参数选择依据lr0.2Diabetes特征已标准化学习率可比西瓜数据集略大max_iter500样本量大需更多迭代错误率计算用y_pred ! y_test而非损失函数值严格对应“错误率”定义。3.3 留一法LOO实现当kn时的极端情况留一法是k折CV中kn的特例。对442样本的Diabetes需训练442次——计算量巨大但能暴露模型稳定性def leave_one_out(X, y, lr0.1, max_iter500): 留一法每次留1个样本作测试 n len(X) errors [] for i in range(n): # 构造训练集除第i个样本外所有样本 X_train np.vstack([X[:i], X[i1:]]) y_train np.hstack([y[:i], y[i1:]]) X_test, y_test X[i:i1], y[i:i1] # 训练并预测 w, b, _ logistic_regression(X_train, y_train, lrlr, max_itermax_iter) z_test X_test w b y_pred (sigmoid(z_test) 0.5).astype(int)[0] errors.append(y_pred ! y_test[0]) if (i1) % 50 0: print(f已完成{i1}/{n}次留一训练...) return np.array(errors) # 执行LOO谨慎运行约需2-3分钟 # errors_loo leave_one_out(X_dia, y_dia, lr0.2) # print(fLOO错误率{errors_loo.mean():.4f})性能权衡LOO方差小但偏差大训练集几乎全量模型过拟合风险高10折CV方差略大但偏差更小是工业界默认选择本资源中3.4.1.py已预跑LOO结果错误率0.283供你直接对比。4. 线性判别分析LDA手撕从投影方向到分类边界4.1 LDA核心思想再确认为什么它不是“另一个分类器”LDA线性判别分析常被误认为是分类算法实则是监督降维方法。西瓜书公式(3.39)给出最优投影方向$$ \boldsymbol{w}^* \mathbf{S}_w^{-1}(\boldsymbol{\mu}_0 - \boldsymbol{\mu}_1) $$其中Sw是类内散度矩阵μ₀、μ₁是两类均值。关键点LDA先将高维数据投影到一维或低维再在该子空间上用阈值分类。这与对率回归直接在原始空间找超平面有本质区别。4.2 LDA投影与分类在西瓜数据集上复现图3.5def lda_fit(X, y): 计算LDA投影方向w和阈值w0 # 分别提取正负样本 X0 X[y 0] X1 X[y 1] # 计算各类均值 mu0 np.mean(X0, axis0) mu1 np.mean(X1, axis0) # 计算类内散度矩阵Sw Σ0 Σ1 S0 np.cov(X0.T, biasTrue) # biasTrue 使用n而非n-1 S1 np.cov(X1.T, biasTrue) Sw S0 S1 # 计算投影方向 w Sw^{-1}(mu1 - mu0) try: w np.linalg.inv(Sw) (mu1 - mu0) except np.linalg.LinAlgError: # Sw奇异时添加微小扰动 w np.linalg.inv(Sw 1e-6 * np.eye(Sw.shape[0])) (mu1 - mu0) # 计算投影后的类中心 proj_mu0 w mu0 proj_mu1 w mu1 # 阈值取两类投影中心的中点 w0 (proj_mu0 proj_mu1) / 2 return w, w0 def lda_predict(X, w, w0): LDA预测投影后比较阈值 proj X w return (proj w0).astype(int) # 在西瓜数据集上运行LDA w_lda, w0_lda lda_fit(X, y) y_pred_lda lda_predict(X, w_lda, w0_lda) accuracy_lda np.mean(y_pred_lda y) print(fLDA在西瓜数据集准确率{accuracy_lda:.4f}) # 可视化LDA投影 def plot_lda_projection(X, y, w, w0): proj X w plt.figure(figsize(10, 4)) # 左图原始二维空间 投影线 plt.subplot(1, 2, 1) plt.scatter(X[y0,0], X[y0,1], cred, markerx, label坏瓜) plt.scatter(X[y1,0], X[y1,1], cgreen, markero, label好瓜) # 绘制投影方向线过原点方向为w xlim plt.xlim() ylim plt.ylim() x_line np.linspace(xlim[0], xlim[1], 100) y_line (w[1]/w[0]) * x_line if w[0] ! 0 else np.full_like(x_line, ylim[0](ylim[1]-ylim[0])/2) plt.plot(x_line, y_line, k--, labelLDA投影方向) plt.xlabel(密度) plt.ylabel(含糖率) plt.title(原始空间LDA投影方向) plt.legend() # 右图一维投影空间 plt.subplot(1, 2, 2) plt.hist(proj[y0], bins10, alpha0.6, label坏瓜投影, colorred) plt.hist(proj[y1], bins10, alpha0.6, label好瓜投影, colorgreen) plt.axvline(w0, colorblack, linestyle--, labelf阈值{w0:.3f}) plt.xlabel(投影值) plt.ylabel(频数) plt.title(投影空间LDA分类阈值) plt.legend() plt.show() plot_lda_projection(X, y, w_lda, w0_lda)关键逻辑np.cov(..., biasTrue)LDA理论推导使用总体协方差除以n非样本协方差除以n-1Sw奇异时加1e-6 * I是经典正则化技巧避免矩阵不可逆阈值w0取两类投影中心中点这是LDA最小化分类错误率的理论最优解。5. 避坑指南5个血泪经验总结的常见问题排查5.1 现象运行3.4.3.py时出现ValueError: shapes (17,2) and (3,) not aligned原因logistic_regression()函数中X w维度不匹配。西瓜数据集X是(17,2)但w被初始化为长度3的向量误加了偏置b进w。解决严格分离权重w长度n_features和偏置b标量所有计算中z X w b不要将b拼入w。5.2 现象Diabetes数据集10折CV错误率高达0.9远高于预期原因未对特征做标准化。Diabetes中“sex”特征是0/1二值而“bmi”是20-40的浮点梯度下降被bmi主导其他特征更新极慢。解决必须在load_diabetes_binary()中加入StandardScaler().fit_transform(X)且scaler需在每折CV内部重新拟合本资源已实现。5.3 现象LDA投影图中两类分布严重重叠阈值分类效果差原因西瓜数据集3.0α本身线性不可分书中明确指出LDA假设类条件概率服从高斯分布且协方差相同该假设在此数据上失效。解决这不是代码bug而是数据本质限制。此时应转向非线性方法如SVM核技巧或接受LDA在此数据上的理论局限——这正是西瓜书用此例的教学意图。5.4 现象sigmoid(z)返回nan后续计算全部中断原因z值过大如700导致exp(-z)下溢为01/(10)得infz过小如-700导致exp(-z)上溢为inf1/(1inf)得0再取log得nan。解决必须在sigmoid()中加入np.clip(z, -500, 500)这是数值稳定的黄金实践所有手写sigmoid函数的标配。5.5 现象3.5.py运行后决策边界是斜线而非西瓜书图3.5的垂直线原因图3.5中LDA投影方向恰好与x轴平行w[0,1]但你的w计算结果是[w1,w2]需将投影值映射回原始坐标系画线。解决决策边界方程为w1*x w2*y w0用y (w0 - w1*x)/w2绘制而非简单画yw0。本资源plot_lda_projection()已正确实现。6. 进阶技巧用梯度检查验证你的反向传播是否写对手写梯度下降最大的隐患是梯度计算错误——公式推导没错但代码实现漏了1/m、忘了转置、符号搞反。西瓜书第三章习题3.6要求“验证梯度”这恰恰是工业界模型开发的基石技能。我一般会在logistic_regression()训练前插入梯度检查模块def gradient_check(X, y, w, b, eps1e-5): 数值梯度 vs 解析梯度对比 # 解析梯度你写的dw, db z X w b p sigmoid(z) dw_analytic (1/len(X)) * X.T (p - y) db_analytic (1/len(X)) * np.sum(p - y) # 数值梯度对w每个分量扰动 dw_numeric np.zeros_like(w) for i in range(len(w)): w_plus w.copy() w_plus[i] eps z_plus X w_plus b loss_plus -np.mean(y * np.log(sigmoid(z_plus) 1e-15) (1-y) * np.log(1 - sigmoid(z_plus) 1e-15)) w_minus w.copy() w_minus[i] - eps z_minus X w_minus b loss_minus -np.mean(y * np.log(sigmoid(z_minus) 1e-15) (1-y) * np.log(1 - sigmoid(z_minus) 1e-15)) dw_numeric[i] (loss_plus - loss_minus) / (2 * eps) # 数值梯度对b扰动 b_plus b eps z_plus X w b_plus loss_plus -np.mean(y * np.log(sigmoid(z_plus) 1e-15) (1-y) * np.log(1 - sigmoid(z_plus) 1e-15)) b_minus b - eps z_minus X w b_minus loss_minus -np.mean(y * np.log(sigmoid(z_minus) 1e-15) (1-y) * np.log(1 - sigmoid(z_minus) 1e-15)) db_numeric (loss_plus - loss_minus) / (2 * eps) # 比较相对误差 rel_error_w np.max(np.abs(dw_analytic - dw_numeric) / np.maximum(np.abs(dw_analytic), np.abs(dw_numeric)) 1e-10) rel_error_b np.abs(db_analytic - db_numeric) / max(np.abs(db_analytic), np.abs(db_numeric), 1e-10) print(f梯度检查结果w相对误差{rel_error_w:.2e}, b相对误差{rel_error_b:.2e}) return rel_error_w 1e-4 and rel_error_b 1e-4 # 在训练前调用 w_init np.random.randn(X.shape[1]) * 0.01 b_init 0.0 is_correct gradient_check(X, y, w_init, b_init) print(f梯度检查通过{is_correct})为什么这招管用数值梯度用(f(xε)-f(x-ε))/2ε逼近导数不依赖公式推导是检验解析梯度的“后悔药”相对误差1e-4是业界通用阈值比绝对误差更鲁棒我从西电带毕业设计时就强制学生在交代码前跑这一段至今没放过一个梯度bug。这个技巧的价值在于它把“我相信我推导对了”变成“机器告诉我确实对了”。当你在3.4.2.py里改了损失函数、在logistics.py里加了L2正则、甚至自己写了个新优化器第一件事永远是跑梯度检查——不是为了炫技而是因为在机器学习里信任代码比信任自己更可靠。从那以后我每次写完梯度都强制走一遍gradient_check哪怕只是3行代码的改动。希望帮到你。本文还有配套的精品资源点击获取
返回列表