企业官网建设流程全解析

简介：这份面向统计学与回归分析课程学习者的实验论文，围绕人口数量预测这一具体问题，系统展示从指标选取、样本确定到模型建立与检验的完整分析路径。压缩包内为1个doc文件，大小约442KB，内容按一元线性回归、多元线性回归和非线性模型组织，便于按章节查阅与复现。已有397人学习下载，适合需要完成课程论文、实验报告或掌握回归建模流程的读者参考。文中涵盖总人口与粮食产量的散点图、理论模型设定、回归诊断，以及数据中心化标准化、逐步回归、多重共线性检测和主成分分析等多元建模环节；后半部分还涉及曲线拟合与Logistic模型，可帮助读者理解不同回归方法在人口预测中的适用条件与结果解读方式。整体结构清晰，兼具课程实验的规范性与方法演示价值。

1. 从一份课程实验说起：回归分析做人口数量预测到底能落地到什么程度

去年帮人看一份统计系的课程实验，题目是「基于回归分析的人口数量预测」，拿 1980 到 2014 年共 35 条年度样本，用粮食产量、GDP、出生率、死亡率去拟合总人口，一路把一元线性、多元线性、逐步回归、主成分分析、三次曲线、Logistic 全跑了一遍。很多人第一反应是：这不就是套公式交作业？真跑一遍就知道，同一个数据集上四个模型的 2014 年预测值能从 136175 万人拉到 142180 万人，差出六千万人，谁对谁错根本不是 R² 说了算。这份材料值钱的地方不在公式，而在它把一个「指标怎么选 → 模型怎么建 → 诊断怎么做 → 结论怎么下」的完整链路摊开了，适合刚上手回归的统计、数据分析方向读者照着复现，也适合已经会调库、但没认真做过模型取舍的人拿去对照。

2. 项目全貌与技术栈：一张 35 行数据表撑起的六套模型

这份课程实验的体量比看上去小，核心数据只有一张表、35 行、6 列，但围绕它展开的建模路径有六条，每条都对应一套独立的工具操作和诊断输出。理解这个项目，先理解它的「数据 — 模型 — 结论」三段结构：数据段固定不变，模型段是六个并列分支，结论段做横向误差对比。这种结构在教学材料里很常见，但恰恰是真实项目里最该学的部分——先穷举模型，再收敛到一个能解释、能外推的。

2.1 目录结构与文件职责

原始材料是 Word 论文加 SPSS 输出截图，直接拿来做二次开发不好用。我按可复现的思路把它重排成下面这个骨架，所有路径、命名都是常规做法，你可以照着自己的习惯改：

population-regression/ ├── data/ │ ├── raw/population_1980_2014.csv # 原始样本，35 行 6 列，不做任何加工 │ └── processed/population_scaled.csv # 标准化后数据，供多元回归和 PCA 用 ├── notebooks/ │ ├── 01_uni_linear.ipynb # 一元线性：散点图 + 模型摘要 + 系数 │ ├── 02_multi_linear.ipynb # 多元线性 + 逐步回归 │ ├── 03_collinearity_pca.ipynb # VIF 诊断 + 主成分分析 │ └── 04_curve_logistic.ipynb # 三次曲线 + Logistic 拟合 ├── spss/ # SPSS 22.0 的 .sav 与导出结果，留档用 ├── matlab/logistic_fit.m # Logistic 曲线拟合脚本 ├── figures/ # 散点图、拟合曲线图，论文插图来源 ├── report/ # 论文正文与图表编号对照 └── requirements.txt

这个分层的关键在于data/raw和data/processed必须分开。原始表里的粮食产量单位是万吨、GDP 单位是亿元、出生率和死亡率是百分比、总人口是万人，量纲差了三四个数量级，任何一次标准化都会覆写数值，如果直接改原表，后面想回头核对原始散点就没依据了。spss/目录是给「论文要交 SPSS 输出截图」这个场景留的，做工程化改造时可以整个删掉，但建议保留一个版本，因为 SPSS 的模型摘要表格式和 Python 输出不完全一致，写报告时对不上会很麻烦。

2.2 字段清单与建模主线

原始样本的字段定义如下，这套命名在后面的所有代码块里保持一致：

六个分支的模型结构、拟合优度和 2014 年预测结果对比如下，这张表是整个项目最有信息量的东西，后面几章都是围绕它展开的：

注意最后一行 Logistic，它在 1995 年前后就被实际人口突破了预测上限，这个模型在这份数据上直接失效，原因第 5 章会说。表的第二列也值得盯一眼：越往后模型的自变量越少，逐步回归从 4 个压到 3 个反而误差更小，这说明原始的多元回归里确实有变量在添乱。

2.3 工具链分工：SPSS、MATLAB 与 Python 复算

原始材料的工具链是三段式的：SPSS 22.0 负责一元线性、多元线性、逐步回归、共线性诊断和主成分分析，MATLAB 负责 Logistic 曲线拟合，三次曲线也是 SPSS 的曲线估计模块出的。这套组合在统计系课程里非常典型，SPSS 点几下就出模型摘要和系数表，截图直接进论文，MATLAB 的fit函数处理非线性拟合比 SPSS 顺手。

问题在于可复现性。SPSS 的.spv输出是二进制，MATLAB 脚本里的初值靠手填，换台机器、换个数据区间，结果就不一样了。我一般会做一层 Python 复算，用statsmodels把每个模型的系数和显著性打印出来，和 SPSS 截图逐个对，对不上的地方就是理解有偏差的地方。下面这条命令装齐复算需要的依赖：

pip install pandas numpy statsmodels scikit-learn scipy matplotlib

statsmodels负责带完整统计量的 OLS（包括 t 值、p 值、R²、AIC），scikit-learn负责标准化和 PCA，scipy的curve_fit替掉 MATLAB 的fit。这么分是因为statsmodels的summary()输出格式最接近 SPSS 的模型摘要表，对账成本最低；而sklearn的StandardScaler和PCA在数值实现上更稳，处理量纲差异大的多变量数据不容易出幺蛾子。

数据读入这一步没什么花活，但列名和类型要一次定死，不然后面每个 notebook 都要重复清洗：

import pandas as pd # 六个字段全部按原始单位读入，不做任何转换，保持和统计年鉴一致 df = pd.read_csv("data/raw/population_1980_2014.csv") # 显式声明数值类型，避免出生率、死亡率被识别成字符串 num_cols = ["year", "grain_output", "gdp", "birth_rate", "death_rate", "total_pop"] df[num_cols] = df[num_cols].apply(pd.to_numeric, errors="coerce") # 时间自变量 t：以 1980 年为基准年记 1，2014 年记 35，供三次曲线和 Logistic 使用 df["t"] = df["year"] - 1979 print(df.shape) # 期望输出 (35, 7) print(df.isna().sum()) # 期望全部为 0

这段代码有两个不起眼但重要的点。errors="coerce"会把无法转换的值变成 NaN 而不是直接抛异常，配合后面的isna().sum()就能一次性看出原始表里有没有空行、有没有千分位逗号混进来。t的基准年定在 1980 记 1 而不是记 0，是因为 Logistic 模型的标准形式里带(x - 1980)这一项，基准年和原论文对齐后，两个模型的系数可以直接横向比。

提示：df.shape打印出来如果是 (35, 6) 而不是 (35, 7)，说明t列没加上，检查一下year列的类型是不是被读成了字符串。

3. 一元线性回归跑通：从散点图到 2014 年预测值

一元线性是整个项目的入口，它的作用不是给出准确预测，而是先回答一个更基础的问题：粮食产量和总人口之间到底有没有线性关系，关系有多强。这一步做不扎实，后面多元回归加变量就是盲加。

3.1 指标选择与样本口径

原始材料在指标选择上给了一段经济学的定性分析，大意是经济因素决定了人口的增殖条件和生存条件，粮食产量、GDP、出生率、死亡率都是候选变量。这段论述本身没问题，但落到建模上，选变量的真正依据是「可得性 + 相关性 + 不重复」。

可得性最好理解：这五个指标都是年度统计口径，35 年连续无缺失。相关性要先看散点图，不能凭感觉。不重复指的是变量之间不能互相替代——粮食产量和 GDP 在 1980 到 2014 这个区间里高度同向增长，两个都放进去就会触发多重共线性，这正是第 4 章要解决的问题。

样本区间定在 1980 到 2014 年，35 个观测值。这个样本量对一元线性回归够用（一般要求 n ≥ 30），对多元回归就偏紧了——4 个自变量加常数项共 5 个待估参数，35 个样本的自由度是 30，还在可接受范围，但如果再加两个变量就会明显过拟合。这个约束在后面的逐步回归里体现得很直接：模型宁愿剔除变量也不愿意保留。

注意：样本区间跨越了人口增长从快速到放缓的完整过程，这段数据的趋势不是一条直线，这正是后面三次曲线能赢过线性模型的根本原因。选样本区间时如果只取其中十年，结论会完全不同。

3.2 散点图与理论模型设定

散点图是这一步唯一不能跳的环节。原论文的图 1-1 显示，粮食产量横轴、总人口纵轴，样本点大致落在一条从左下到右上的带状区域里，没有明显的弯曲。基于这个观察，理论模型设成一元线性形式：

$$ y = \beta_0 + \beta_1 x_1 + \varepsilon $$

其中 y 是总人口，x1 是粮食产量，β0 是截距，β1 是回归系数，ε 是随机误差项。用 Python 把散点画出来并拟合：

import matplotlib.pyplot as plt import statsmodels.api as sm # 自变量需要手动加常数项，statsmodels 默认不包含截距 X = sm.add_constant(df["grain_output"]) y = df["total_pop"] model_uni = sm.OLS(y, X).fit() print(model_uni.summary()) # 散点 + 拟合直线，用于核对原论文图 1-1 的形态 fig, ax = plt.subplots(figsize=(7, 5)) ax.scatter(df["grain_output"], df["total_pop"], s=18, label="sample") ax.plot(df["grain_output"], model_uni.predict(X), color="crimson", label="OLS fit") ax.set_xlabel("grain output (10k tons)") ax.set_ylabel("total population (10k)") ax.legend() fig.savefig("figures/fig1_1_scatter.png", dpi=150)

sm.add_constant这行是最容易漏的，漏了之后模型强制过原点，截距变成 0，R² 会莫名其妙掉一大截，而且系数解释完全不同。model_uni.summary()的输出里有三块要重点看：R² 和调整后 R²、F 统计量的 p 值、每个系数的 t 值和 p 值。

拟合结果和原论文一致：R = 0.909，R² = 0.826，调整后 R² = 0.821，标准误 4985.997。截距 53054.404（t = 9.646，p < 0.001），粮食产量系数 1.468（t = 12.513，p < 0.001）。写成方程就是：

$$ \hat{y} = 53054.404 + 1.468 x_1 $$

系数 1.468 的含义是：粮食产量每增加 1 万吨，总人口平均增加 1.468 万人。这个解释在统计上成立，但在业务上很别扭——粮食增产不会直接导致人口增长，两者都是被时间趋势带着走的。这就是典型的伪回归风险，一元线性在这里只能当「相关性强弱」的探针用，不能当因果结论。

3.3 回归诊断与预测偏差

把 2014 年的粮食产量 60702.6 万吨代进去，得到预测值 142165.84 万人，而实际总人口是 136782 万人，绝对误差 5383.84 万人，相对误差约 3.9%。对一个只用一个自变量的一元模型来说，这个误差不算离谱，但方向很值得注意：预测值系统性偏高。

偏高的原因从残差图能看出来。样本前期（1980 年代到 1990 年代）实际人口低于拟合线，后期（2000 年以后）实际人口也低于拟合线，中间段高于拟合线，残差呈现明显的倒 U 形。这说明真实的增长曲线是「先快后慢再趋平」，而直线只能取一个平均斜率，早期低估、晚期高估，外推一年就多算一大截。

诊断这一步我一般固定看四个量，缺一不可：

import numpy as np resid = model_uni.resid fitted = model_uni.fittedvalues # 1. 残差正态性：偏度和峰度接近 0 才谈得上 t 检验有效 print("skew:", round(float(pd.Series(resid).skew()), 3)) print("kurt:", round(float(pd.Series(resid).kurt()), 3)) # 2. 残差与拟合值的关系：理想情况是随机散布，出现喇叭形说明异方差 print("corr(resid, fitted):", round(float(np.corrcoef(resid, fitted)[0, 1]), 3)) # 3. Durbin-Watson：时间序列数据必看，接近 2 才说明残差无自相关 from statsmodels.stats.stattools import durbin_watson print("DW:", round(durbin_watson(resid), 3)) # 4. 单点外推：把 2014 年粮食产量代入，看误差量级 x_2014 = sm.add_constant(pd.DataFrame({"grain_output": [60702.6]}), has_constant="add") print("pred 2014:", round(float(model_uni.predict(x_2014)[0]), 2))

四个量的判读标准：偏度绝对值小于 1、峰度绝对值小于 2 算近似正态；corr(resid, fitted)绝对值大于 0.3 就要怀疑异方差；Durbin-Watson 落在 1.5 到 2.5 之间算可接受，明显小于 1.5 说明残差正自相关，标准误会被低估，p 值不可信。这份数据是逐年递增的时间序列，DW 值大概率偏低，这也是为什么一元线性模型的显著性看起来漂亮，预测却不准。

提示：sm.add_constant对单行 DataFrame 要显式传has_constant="add"，否则当数据里没有常数列时会静默跳过，预测结果直接偏移一个截距。

4. 多元线性回归与降维：标准化、逐步回归、VIF 与 PCA 的取舍

一元模型只用了粮食产量一个变量，误差 5383 万人。把 GDP、出生率、死亡率都放进去，R² 从 0.826 跳到 0.944，但 2014 年预测误差只降到 1400 万人，而且四个变量里有两个不显著。这一章的核心不是「怎么把变量加进去」，而是「加进去之后怎么把没用的踢出来」。

4.1 数据中心化与标准化

多元回归的第一个动作是中心化和标准化。原因很直白：粮食产量在 3 万到 6 万之间，GDP 在 4000 到 60 万之间，出生率死亡率在 6 到 22 之间，量纲差了三四个数量级，直接进模型会让数值求解的矩阵条件数变得很差，系数估计不稳定。中心化减均值，标准化再除标准差：

$$ z_{ij} = \frac{x_{ij} - \bar{x}_j}{s_j} $$

手写一版便于理解，再用 sklearn 对一遍：

# 手写版：逐列减均值除标准差，ddof=1 与 SPSS 的样本标准差口径一致 cols = ["grain_output", "gdp", "birth_rate", "death_rate"] Z = df[cols].copy() for c in cols: Z[c] = (df[c] - df[c].mean()) / df[c].std(ddof=1) # sklearn 版：ddof 默认是 0（总体标准差），和 SPSS 会有一点差异 from sklearn.preprocessing import StandardScaler scaler = StandardScaler() Z_sk = pd.DataFrame(scaler.fit_transform(df[cols]), columns=cols) print("max abs diff:", float((Z - Z_sk).abs().max().max()))

两版结果的差异来自ddof：pandas.std()默认 ddof=1（样本标准差），StandardScaler默认 ddof=0（总体标准差）。35 个样本下两者差不到 1.5%，对标准化系数的影响微乎其微，但如果你要拿标准化系数和 SPSS 的 Beta 列逐位对账，就得统一口径，否则最后一位永远对不上，白白浪费时间排查。

标准化之后系数变成「标准化回归系数」，可以直接比大小判断哪个变量对总人口的影响更强。这一步很关键，因为原始量纲下的系数不能横向比较——出生率的系数是 -1460.798，看着比粮食产量的 1.045 大三个数量级，但出生率本身只波动十几个点，粮食产量波动几万吨，量纲一除，实际影响力反而是出生率更大。

4.2 多元回归模型建立与逐步回归

四个变量全进的模型形式是：

$$ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \varepsilon $$

跑出来的结果：R = 0.972，R² = 0.944，调整后 R² = 0.936，F = 126.295，整体 p < 0.001，模型整体显著。系数表里粮食产量 p < 0.001、出生率 p < 0.001，但 GDP 的 p = 0.531、死亡率 p = 0.275，两个都没过 0.05。也就是说，模型整体显著，但有两个变量是「搭便车」进来的，它们的贡献被其他变量吸收了。

逐步回归要解决的就是这个。它的逻辑是「有进有出」：每次引入一个变量，引入后对所有已引入变量重新做显著性检验，哪个不显著就剔掉，直到没有变量能进也没有变量该出。手写一版双向逐步回归，用 p 值做进出标准：

def stepwise_ols(X, y, p_enter=0.05, p_remove=0.10): """双向逐步回归：p 值小于 p_enter 进，大于 p_remove 出""" included = [] candidates = list(X.columns) while True: changed = False # 前向：从候选里挑 p 值最小且显著的加进来 excluded = [c for c in candidates if c not in included] best_p, best_var = 1.0, None for c in excluded: trial = included + [c] m = sm.OLS(y, sm.add_constant(X[trial])).fit() p = m.pvalues[c] if p < best_p: best_p, best_var = p, c if best_var is not None and best_p < p_enter: included.append(best_var) changed = True # 后向：已引入变量里 p 值超过阈值的剔掉 if included: m = sm.OLS(y, sm.add_constant(X[included])).fit() worst_p = m.pvalues.drop("const").max() if worst_p > p_remove: included.remove(m.pvalues.drop("const").idxmax()) changed = True if not changed: break return included, sm.OLS(y, sm.add_constant(X[included])).fit() X_all = df[cols] selected, model_step = stepwise_ols(X_all, df["total_pop"]) print("selected:", selected) print(model_step.summary())

这段代码有两个参数要盯住：p_enter=0.05是入选门槛，p_remove=0.10是剔除门槛。剔除阈值必须比入选阈值宽，否则会出现「刚进来就被踢、踢完又进来」的死循环，这是逐步回归实现里最常见的翻车点。m.pvalues.drop("const")排除常数项，因为截距不参与变量筛选。

运行结果是剔掉 GDP，保留粮食产量、出生率、死亡率三个变量，R² 反而升到 0.972，调整后 R² 升到 0.952。方程写成：

$$ \hat{y} = 134678.878 + 1.012 x_1 - 1366.202 x_3 - 5681.557 x_4 $$

2014 年预测值 138506.22 万人，绝对误差 172.22 万人，比四变量模型的 1400.24 万人小了一个数量级。变量少了、误差反而小了，这不是玄学，是剔掉了共线性带来的方差膨胀。

不过原论文在这里下了一个我认为站不住的结论，它说「不用逐步回归的预测值更接近真实值」。按它自己给的数据，四变量多元回归误差 1400.24 万人，逐步回归误差 172.22 万人，明明是逐步回归更接近。写材料时这种地方最容易糊过去，但如果你照着复现，一定要以自己跑出来的数为准。

4.3 多重共线性检测与 VIF

GDP 被剔除不是偶然，它和粮食产量、总人口在时间维度上同向增长，共线性很严重。检测共线性最常用的是方差膨胀因子：

$$ \mathrm{VIF}_j = \frac{1}{1 - R_j^2} $$

其中 R²_j 是把第 j 个自变量当因变量、其余自变量当自变量回归得到的决定系数。经验判断：VIF 大于 10 说明共线性严重，5 到 10 之间要警惕，小于 5 基本可以放心。SPSS 的共线性诊断表里还给特征值和条件指数，特征值有 3 个接近 0、条件指数有 3 个大于 10，这两个信号同时出现，基本可以确诊。

Python 里手算 VIF 只要几行：

from statsmodels.stats.outliers_influence import variance_inflation_factor def vif_table(X): Xc = sm.add_constant(X) out = [] for i, col in enumerate(Xc.columns): if col == "const": continue out.append({ "var": col, "VIF": round(variance_inflation_factor(Xc.values, i), 3), }) return pd.DataFrame(out).sort_values("VIF", ascending=False) print(vif_table(df[cols])) # 四变量全进，看谁是重灾区 print(vif_table(df[["grain_output", "birth_rate", "death_rate"]])) # 剔除 GDP 后复查

variance_inflation_factor的第二个参数是列索引，不是列名，且必须包含常数列，否则计算出的 VIF 会偏大。四变量全进时，GDP 的 VIF 一般会明显偏高；剔掉 GDP 后，剩下三个变量的 VIF 全部落在 10 以下（原论文的表 2-6 给的具体值是 2.868、4.146、1.808），说明不需要再剔变量了。

这里有个常被忽略的细节：VIF 只能发现「线性」共线性，对非线性关系和交互项无能为力。如果两个变量是二次关系，VIF 可能显示正常，但模型依然不稳。所以 VIF 通过不代表模型一定健康，还要回头看残差。

4.4 主成分分析降维建模

主成分分析是另一条处理共线性的路：不剔变量，而是把相关的原始变量正交变换成一组互不相关的主成分，用前几个主成分做回归。原论文用 SPSS 跑出来前两个主成分的累计贡献率 94.18%，第一个主成分单独占 76.004%，第二个占 18.176%。

Python 复现：

from sklearn.decomposition import PCA X_std = StandardScaler().fit_transform(df[cols]) # PCA 前必须标准化 pca = PCA(n_components=4) scores = pca.fit_transform(X_std) # 累计方差贡献率，对应 SPSS 的「解释的总方差」表 cum = np.cumsum(pca.explained_variance_ratio_) print("explained ratio:", np.round(pca.explained_variance_ratio_, 4)) print("cumulative:", np.round(cum, 4)) # 取前两个主成分做回归 X_pca = sm.add_constant(pd.DataFrame(scores[:, :2], columns=["PC1", "PC2"])) model_pca = sm.OLS(df["total_pop"], X_pca).fit() print(model_pca.summary())

n_components=4先算全量再按累计贡献率截断，比直接写n_components=2好——后者拿不到完整的贡献率表，没法判断到底该取几个。fit_transform输入的必须是标准化后的矩阵，如果拿原始量纲的数据直接送进去，方差最大的 GDP 会独占第一主成分，PCA 就退化成了「挑出 GDP」。

按原论文的路径，用前两个主成分回归，得到方程 ŷ = 53075.711 + 1.468x1 + 2.494e-5x2，2014 年预测值 142180.53 万人，误差 5398.53 万人，比逐步回归差了一个数量级。原因有两个：一是主成分是原始变量的线性组合，回归完还要「翻译」回原始变量，这个翻译过程会丢信息；二是被丢掉的那 5.82% 方差里可能恰好包含了对人口预测最关键的那部分变异。主成分分析降维是为「变量太多、样本太少」准备的，这份数据只有 4 个变量、35 个样本，用它属于杀鸡用牛刀，效果反而不如老老实实做逐步回归。

注意：主成分回归的可解释性天然弱于原始变量回归。如果下游要给业务方解释「出生率每降一个千分点，人口怎么变」，用主成分模型基本没法回答，这时候逐步回归是更好的选择。

5. 非线性模型与踩坑排查：三次曲线为什么赢过 Logistic

前四章的模型都是线性的，自变量要么是粮食产量，要么是几个经济指标。这一章换思路，把时间 t 当唯一自变量，直接拟合人口随时间的增长曲线。这个转向很关键：人口的长期趋势本质是时间的函数，经济指标只是伴生现象。

5.1 时间变量曲线拟合与三次函数

以 1980 年为基准年记 t = 1，2014 年记 t = 35，画出总人口对 t 的散点图，形态是「先快速上升、中段放缓、末端趋平」，一条直线显然不够。原论文试了线性、复合、三次等多种曲线形式，三次曲线的 R = 1.000，R² = 0.999，调整后 R² 也接近 1，明显优于其他形式。

三次模型形式：

$$ y = \beta_0 + \beta_1 t + \beta_2 t^2 + \beta_3 t^3 + \varepsilon $$

Python 里不用专门的曲线拟合库，直接构造 t 的幂次列做 OLS 就行：

# 构造三次多项式特征，t 已在前面的数据准备中算好 df["t2"] = df["t"] ** 2 df["t3"] = df["t"] ** 3 X_cubic = sm.add_constant(df[["t", "t2", "t3"]]) model_cubic = sm.OLS(df["total_pop"], X_cubic).fit() print(model_cubic.summary()) # 预测 2014 年（t = 35） t_new = pd.DataFrame({"const": [1.0], "t": [35], "t2": [35**2], "t3": [35**3]}) print("pred 2014:", round(float(model_cubic.predict(t_new)[0]), 2))

t2、t3两列必须显式构造，statsmodels的公式接口里虽然可以写y ~ t + I(t**2) + I(t**3)，但显式建列在需要往模型里加新样本预测时更方便——直接构造一行同样结构的数据就能喂进去。

拟合结果：截距 96103.588，t 的系数 1849.094（t 值 28.230，p < 0.001），t² 的系数 -17.683（t 值 -4.214，p < 0.001），t³ 的系数 -0.070（t 值 -0.908，p = 0.371）。方程是：

$$ \hat{y} = 96103.588 + 1849.094 t - 17.683 t^2 - 0.070 t^3 $$

2014 年预测值 136174.95 万人，绝对误差 607.05 万人，相对误差约 0.44%，在所有模型里拟合优度最高。但这里有个坑必须指出来：t³ 项的 p 值是 0.371，没通过显著性检验。原论文写「每个不同次幂的 t 也经过了显著性检验」，这和它自己给的 0.371 对不上。实践中的处理是：留在模型里（因为去掉三次项后曲线形态会明显变差），但在报告里如实标注该项不显著，别硬说全显著。

还有一处更值得討論：三次曲线误差 607.05 万人，逐步回归误差 172.22 万人，单看这一个点，逐步回归其实更准。原论文说三次曲线「远小于」其他模型，这个判断至少对逐步回归不成立。这说明单点误差不能作为模型选择的唯一标准——三次曲线的价值在于它拟合出了「增速递减到趋平」的趋势形态，往前外推几年的走势是合理的；逐步回归虽然在这一个点上更准，但它靠的是经济指标，而这些指标本身也在变，外推风险更大。

5.2 Logistic 模型与初值设定

Logistic 模型假设人口增长有上限，增长率随人口规模上升而下降：

$$ p_t = \frac{p_m}{1 + e^{a + bt}} $$

其中 p_m 是人口上限（环境容纳量），a、b 是待定参数。原论文用 MATLAB 的fit函数拟合，脚本大致是：

x = [2014 2013 ... 1980]'; % 年份倒序，转列向量 y = [136782 136072 ... 98705]'; % 对应总人口，单位万人 st_ = [2000 127627 1]; % 三个参数的初值 [a b k] ft_ = fittype('a/(1+b*exp(-k*(x-1980)))', ... 'dependent', {'y'}, 'independent', {'x'}, ... 'coefficients', {'a', 'b', 'k'}); cf_ = fit(x, y, ft_, 'Startpoint', st_);

fittype里把基准年写死成(x-1980)，和前面 t 的基准年定义一致。Startpoint给的是[a b k]三个初值，原脚本里 a 填 2000、b 填 127627、k 填 1，注释写「任意取数」。这个「任意取数」是后面出问题的根源：非线性最小二乘对初值极度敏感，初值离真解太远时，算法可能收敛到局部极小，也可能发散。

Python 用scipy.optimize.curve_fit复现：

from scipy.optimize import curve_fit def logistic(t, pm, a, b): # pm 为人口上限，a、b 为形状参数，t 已按 1980 基准年偏移 return pm / (1 + np.exp(a + b * t)) # 初值：pm 取样本最大值上浮 10%，a、b 给小量，不要动辄上万 p0 = [float(df["total_pop"].max()) * 1.1, -5.0, 0.05] popt, pcov = curve_fit(logistic, df["t"], df["total_pop"], p0=p0, maxfev=20000) pm, a, b = popt print(f"pm={pm:.1f}, a={a:.4f}, b={b:.4f}") print("pred 2014:", round(float(logistic(35, *popt)), 2))

maxfev=20000是最大迭代次数，默认值在参数尺度悬殊时经常不够，会抛RuntimeError: Optimal parameters not found。p0的选取也有讲究：把 pm 的初值设在样本最大值略微上浮的位置，a、b 给量级合理的小数，比原脚本那种「2000、127627、1」的乱填成功率高得多。

按原论文的结果，MATLAB 拟合出 a = 1.217e5、b = 1.176e5、k = 13.21，置信区间宽到没有意义（b 的区间是 -4.741e9 到 4.742e9）。更致命的是上限估计值 1.217e5 万人，而 1995 年的实际人口就已经超过这个数，模型从一开始就被数据否定了。原论文把原因归结为「数据过少、抽样误差大」，这个解释只对了一半，更本质的原因是：Logistic 假设增长有上限且 S 形对称，而这份样本区间恰好只覆盖了 S 形的上升段和拐点附近的一小段，模型参数根本没法识别。

5.3 踩坑排查：五条血泪经验

这份材料从复现到对账，踩过的坑基本都在下面五条里，每条都是「现象 → 原因 → 解决」的完整链路。

现象一：标准化之后的标准化系数和 SPSS 对不上最后一位。原因：pandas.std()默认 ddof=1，StandardScaler默认 ddof=0，35 个样本下两者差约 1.4%，系数传到第四位小数就分叉了。解决：对账阶段统一用ddof=1，手写标准化或者给StandardScaler外面包一层pandas的std，别混用两套口径。

现象二：逐步回归跑出来变量反复进出，程序卡死不退出。原因：入选阈值和剔除阈值设成了同一个值（比如都是 0.05），一个变量刚进来 p 值略超 0.05 就被踢，踢完再算又低于 0.05 又被加回来，形成死循环。解决：剔除阈值必须严格大于入选阈值，常规设法是入选 0.05、剔除 0.10，并且每轮循环记录一次已入选集合，发现集合重复就强制 break。

现象三：VIF 算出来全是 1.0，明显不对。原因：调用variance_inflation_factor时传入的矩阵没有加常数列，或者传了常数列但索引i用错了列的位置。没有常数项时，各个自变量的 VIF 会被严重低估。解决：先sm.add_constant再用Xc.values和列索引i，并且在函数里显式跳过const列，不要用列名的位置去猜索引。

现象四：三次曲线 R² = 0.999，但拿它预测 2020 年人口，结果反而往下掉。原因：三次项系数是负的（-0.070），当 t 继续增大，负的三次项会逐渐压过一次项和二次项，曲线在 t = 35 之后不久就掉头向下。这不是数据错了，是三次多项式的外推行为本身不可控。解决：三次曲线只用于拟合样本区间内的趋势，外推必须限制在样本外一到两年，再远就要换 Logistic 或灰色模型，并且在报告里明确标注外推区间。

现象五：Logistic 拟合时报Optimal parameters not found，或者算出的上限比历史实际值还小。原因有两个，一是初值给得太离谱，fit/curve_fit在参数尺度相差几个数量级时很容易卡在局部极小；二是样本区间太短，只覆盖了 S 形曲线的一小段，上限参数 p_m 根本没被数据约束住。解决：初值按「p_m 取样本最大值上浮 10% 到 20%、a 取负数、b 取小正数」设定，maxfev调大；如果换了几组初值上限还是被实际值突破，就承认这个数据不支持 Logistic，换模型，别硬拟合。

提示：这五条里最容易反复踩的是第二条和第四条。逐步回归的死循环不会报错，只会一直跑；三次曲线的外推掉头也不会报警，只有你把预测值画出来才会看见。两个问题都建议加一行断言或者画一张外推图兜底。

6. 二次开发：把这套课程实验改成能复用的预测脚本

原始材料是论文形态，每个模型占一节，代码零散在 SPSS 截图和 MATLAB 脚本里。真要拿去用，第一步是把六个模型抽成统一接口，让「换数据、换模型、出对比表」变成三行调用。

我的做法是定义一个fit_and_forecast约定，所有模型实现同一个函数签名，输入是训练集和预测年份，输出是预测值、模型对象和诊断指标。这样横向对比就不用一个个 notebook 手动抄数了：

from dataclasses import dataclass @dataclass class ForecastResult: name: str pred: float r2: float params: dict note: str = "" def fit_uni_linear(train, target_year_grain): X = sm.add_constant(train["grain_output"]) m = sm.OLS(train["total_pop"], X).fit() pred = float(m.params["const"] + m.params["grain_output"] * target_year_grain) return ForecastResult("uni_linear", pred, m.rsquared, dict(m.params)) def fit_stepwise(train, target_row): cols = ["grain_output", "birth_rate", "death_rate"] m = sm.OLS(train["total_pop"], sm.add_constant(train[cols])).fit() pred = float(m.predict(sm.add_constant(target_row[cols], has_constant="add"))[0]) return ForecastResult("stepwise", pred, m.rsquared, dict(m.params))

dataclass的好处是把「模型名、预测值、R²、参数、备注」绑成一个对象，后面汇总成 DataFrame 只要一行pd.DataFrame([r.__dict__ for r in results])。target_row传一个单行 DataFrame，字段名和训练集保持一致，这样预测阶段不需要再关心每个模型用了哪几个变量。

第二步是把「模型选择」这件事从人工判断变成可配置的规则。现在这套实验里，一元线性只是探针、多元线性是基线、逐步回归是主力、三次曲线是趋势参照、Logistic 是反例，这个定位应该写进代码而不是写进论文。我一般会加一个MODEL_REGISTRY字典，把每个模型的适用条件和已知缺陷挂上去，跑完自动打印一句结论，避免下次换数据时又从头判断一遍。

第三步，也是最值钱的一步：把误差评估从「单点绝对误差」换成「滚动外推误差」。这份数据只有 2014 年一个真值可以对照，六个模型比来比去其实样本量是 1，统计上没有意义。改成用 2005 到 2013 年做训练、逐年外推预测 2006 到 2014 年，得到 9 个误差点，再算 MAE 和 RMSE，模型的排序才站得住。我试过一次，排序确实变了——逐步回归在单点上最准，但滚动外推的稳定性不如三次曲线，原因还是那句话：经济指标自身也在漂移，靠它们预测人口等于把两个不确定性叠在一起。

最后一个技巧是给外推加护栏。三次曲线在 t 超过 35 之后会掉头，Logistic 的上限可能被历史值突破，这两类问题在代码层面都能提前拦住：拟合完先算一阶导数，在预测区间内检查符号是否变化，变了就抛警告；Logistic 拟合完检查pm是否大于训练集最大值，小于就直接标记为不可用。这类检查写起来不超过二十行，但能省掉后面无数次「这个数怎么这么怪」的排查。

从那以后我每次做回归类项目，都会在拟合完的第一时间把「模型预测值、真值、残差、外推区间」四列打在同一张表里，不看 R²，先看误差的方向和量级对不对。这套习惯就是从这份实验里被三次曲线的 0.999 和 Logistic 的上限突破教出来的。希望帮到你。

本文还有配套的精品资源，点击获取