简介：围绕回归分析的人口数量预测形成的一篇统计学课程实验论文，面向统计、经济与数据分析方向的学习者，可用于课程实验报告参考、回归建模入门及人口与经济指标关系分析练习。文档以人口数量预测为主线，串联一元线性回归、多元线性回归与非线性模型：从指标选择、样本确定、散点图观察，到理论模型设定、参数估计与回归诊断；进而讨论数据中心化标准化、逐步回归、多重共线性检测和主成分分析，并延伸到曲线拟合与Logistic模型。压缩包仅含1个doc文件，大小约442KB，篇幅紧凑但目录完整，包含摘要、前言、各章分析、结论与参考文献，便于按章节查阅建模思路与论文结构。已有396人学习，适合需要参考完整统计实验论文框架、梳理回归分析流程或撰写同类课程的读者使用。

1. 拿到这份资源时我先翻到了结论页

一份统计系的课程实验论文，封面写着「基于回归分析的人口数量预测」，正文 17 页，前后用了 SPSS 22.0 和 MATLAB 两套工具，样本是 1980 到 2014 年共 35 条年度数据。多数人看到「课程论文」四个字就划走了，但这份东西其实是一条完整的建模流水线：一元线性回归、四元线性回归、逐步回归、多重共线性诊断、主成分分析、三次曲线、Logistic 曲线，七种模型在同一个数据集上跑了一遍，还给了每个模型对 2014 年的预测值和误差。对刚接触回归分析、想找一份能照着复现的样例的人来说，它省掉了「指标怎么选、变量怎么筛、模型怎么比」这三步最耗时间的摸索；对已经会跑 OLS 的人，它更像一份反面教材合集，里面好几个处理方式值得拿来当踩坑清单。人口数量预测本身不难，难的是判断哪个模型的结论能信，这份论文恰好把「不能信的地方」也一并留在了纸面上。

2. 数据结构与预处理：35 年样本、5 个指标怎么组织成能跑的数据

这份实验的数据结构非常朴素：一张宽表，35 行，6 列。真正动手复现时，第一件要做的事不是急着 fit 模型，而是把列名、单位、时间基准这三样东西固定下来，后面所有模型都依赖它。原文在指标选择上给的理由是「经济因素决定增殖条件和生存条件，粮食产量、出生率、死亡率同样影响人口增长」，落到操作层面就是四个自变量加一个因变量，符号记作 x1 到 x4 和 y。

2.1 指标选择与字段口径

论文里定的四个自变量是粮食产量、GDP、出生率、死亡率，因变量是总人口。这套指标组合在人口预测场景里不算最优——少了年龄结构、城镇化率、人均预期寿命这些更直接的量——但它的好处是数据都能从公开统计年鉴里找到连续的年度值，35 年不断档，这对课程实验来说比指标的理论完备性重要得多。

字段口径必须提前对齐，否则后面系数会莫名其妙地放大或缩小 100 倍：

表里最关键的一行是出生率和死亡率的单位。论文表头写的是「%」，但数值是 12.37、7.16 这种量级，稍微有点常识就知道这是千分率。这个小小的标注错误在后面会直接体现为回归系数看起来大得离谱，我在第 5 章会专门讲这个坑。

2.2 样本表的列结构与数据装载

把表格抄进 CSV 之后，第一段代码只做三件事：改列名、设索引、做类型和缺失检查。别小看这几行，从 Excel 导出的 CSV 经常把带千分位逗号的数值列读成字符串，直接进模型会报类型错误，排查起来很浪费时间。

import pandas as pd # 原始表头是中文带括号的，先做一次显式映射，避免列名里混进全角符号 COLS = { "年份": "year", "粮食产量（万吨）": "grain", "GDP（亿元）": "gdp", "出生率": "birth_rate", "死亡率": "death_rate", "总人口（万人）": "population", } df = pd.read_csv("population_1980_2014.csv", encoding="utf-8") df = df.rename(columns=COLS).set_index("year").sort_index() df = df.astype(float) # 防千分位逗号被读成 object print(df.shape) # 期望 (35, 5) print(df.isna().sum()) # 每列缺失数，全 0 才继续 print(df.describe().T[["min", "max"]])

set_index("year")这一步是为了让时间成为唯一主键，后面按年份切片、算时间变量 t 都不用再关心行号。astype(float)放在读入之后立刻执行，比等到报错再回头查省事。describe()打出来的 min/max 是第二道防线：如果 birth_rate 的最大值出现 123.7 这种数，说明单位已经错了，当场就能发现。

2.3 中心化与标准化：量纲对齐到底改了什么

原文把数据标准化单独列了一节，理由是「取消量纲不同、自身变异或者数值相差较大所引起的误差」。这句话在工程上要拆开看：标准化确实让粮食产量（万吨级）和出生率（千分率）在数值上可比，但它不会让模型的拟合优度、F 检验、t 检验发生任何变化——OLS 对自变量的可逆线性变换是不变的。

from sklearn.preprocessing import StandardScaler X = df[["grain", "gdp", "birth_rate", "death_rate"]] z = pd.DataFrame( StandardScaler().fit_transform(X), columns=X.columns, index=X.index, ) print(z.describe().loc[["mean", "std"]].round(3)) # mean≈0, std≈1

StandardScaler做的就是 (x - mean) / std，逐列独立计算。这里要注意两个容易被夸大的点：第一，标准化后重新跑回归，R²、调整 R²、F 值、各个 t 值和 p 值全部不变，变的只有非标准化系数 B 和它的标准误，标准化系数 Beta 本来就与量纲无关；第二，方差膨胀因子 VIF 在标准化前后完全一样，因为 VIF 只依赖自变量之间的相关系数矩阵，标准化不改变相关结构。

真正会变的是条件指数（condition index），它对变量的尺度敏感。原文表 2-5 里条件指数出现 33.798、36.436、86.828 三个大于 10 的值，其中一部分来自共线性，一部分来自 GDP 和其他变量量级差了四五个数量级。所以标准化的实际价值是：让诊断指标更干净，让系数之间能做大小比较，而不是「提高预测精度」。这个区别在写论文结论时很容易说错。

3. 一元线性回归：从散点图到系数表的一条完整链路

一元回归是整份实验里唯一一个能从头到尾讲清楚的模型：一个自变量、一张散点图、一张系数表、一个预测值，链路短，适合拿来验证数据管道通不通。原文选的第一个自变量是粮食产量，理由是散点图上样本点大致落在一条直线附近，随着粮食产量增加总人口增加。

3.1 散点图与模型设定

散点图这一步不能省。它是判断「线性形式是否成立」最便宜的手段，也是发现异常点最快的方式。原文图 1-1 显示 35 个点基本沿右上方向排布，没有明显的离群点，于是设定理论模型为 y = β0 + β1·x1 + ε。

import statsmodels.api as sm y = df["population"] X1 = sm.add_constant(df[["grain"]]) # 必须手动加截距列 m1 = sm.OLS(y, X1).fit() print(m1.summary()) print("截距:", round(m1.params["const"], 3)) print("斜率:", round(m1.params["grain"], 3)) print("R² :", round(m1.rsquared, 3))

sm.add_constant是 statsmodels 里最容易漏的一行。漏了它，模型会被强制过原点，拟合出来的斜率和 R² 都会偏离；而且 summary 里不会报错，属于典型的静默翻车。m1.params返回的是 Series，用列名取值比用位置索引稳，改列名时不会连带出错。

按原文数据跑出来，截距 53054.404（标准误 5500.013，t=9.646，p<0.001），斜率 1.468（标准误 0.117，Beta=0.909，t=12.513，p<0.001），R=0.909，R²=0.826，调整 R²=0.821，估计标准误 4985.997。这个 R² 意味着粮食产量单独一个变量就吃掉了总人口 82.6% 的变异，对一元模型来说已经很高了。

3.2 系数表逐列解读

SPSS 的系数表有五列值得逐一看：非标准化系数 B、标准错误、标准化系数 Beta、t 值、显著性。B 是代回方程用的，Beta 是用来比较不同自变量相对重要性的，t 和 p 是判断这个系数是否显著不为零的。

一元回归里有个巧合：Beta 值 0.909 和模型摘要里的 R 值完全相等。这不是笔误，一元线性回归中标准化斜率就等于皮尔逊相关系数，所以看到这两个数一样不用怀疑表格抄错了。斜率 1.468 的业务含义是「粮食产量每增加 1 万吨，总人口平均增加 1.468 万人」，这个解释在统计上成立，在因果上完全不成立——它只是两个同向增长的时间序列撞在了一起，并不能说明多打粮食就会多生人。写报告时这层边界要交代清楚，否则读者会把相关当因果。

截距 53054.404 在业务上同样讲不通，它对应的是「粮食产量为 0 时的总人口」，属于模型外推，实际意义为零。系数表里每个数都能算出来，但不是每个数都值得解释，这个判断力比会跑代码更重要。

3.3 预测与残差：一条直线的适用边界

把 2014 年的粮食产量 60702.6 万吨代进方程，得到预测值 142165.84 万人，而当年实际总人口是 136782 万人，残差 +5383.84 万人，相对误差约 3.9%。残差绝对值 5384 大约等于一个估计标准误（4986），从统计上看属于正常波动范围。

x_2014 = df.loc[2014, "grain"] y_hat = m1.params["const"] + m1.params["grain"] * x_2014 y_true = df.loc[2014, "population"] print(f"预测 {y_hat:.2f} / 实际 {y_true} / 误差 {y_hat - y_true:.2f}") print(f"相对误差 {(y_hat - y_true) / y_true:.4f}")

这段代码做的事情就是手工代公式，好处是能验证模型对象里的 predict 和手算结果一致，防止加常数项时搞错列顺序。df.loc[2014, ...]依赖前面把 year 设成了索引，这也是那一步的实际收益。

这条线的边界在于外推。粮食产量不会无限增长，总人口也不会随着粮食产量一直线性上升，所以用这条线预测 2025 年、2030 年，结果会迅速失真。一元回归在这份论文里的定位应该是「基线模型」，用来给后面几个更复杂的模型提供一个误差参照，而不是最终结论。

4. 多元回归与非线性的岔路：逐步回归、VIF、主成分、三次曲线

把所有自变量一次性塞进模型，是大多数人拿到数据后的第一反应。这份实验也是这么做的，然后依次遇到了三个问题：变量不显著、变量之间高度相关、模型形式本身不对。原文的解决路径是「四元回归 → 逐步回归 → 多重共线性诊断 → 主成分分析 → 三次曲线 → Logistic」，这条路径本身比任何一个单独的模型都更有参考价值，因为它展示了一个建模者遇到问题时的实际决策顺序。

4.1 四元模型建立与显著性判定

四元模型的公式是 y = β0 + β1x1 + β2x2 + β3x3 + β4x4 + ε，其中 x2 是 GDP，x3 是出生率，x4 是死亡率。跑出来的 R=0.972，R²=0.944，调整 R²=0.936，回归平方和 4448849243.984，残差平方和 264193847.616，F=126.295，p<0.001。整体模型显著，拟合度比一元模型高了一大截。

X4 = sm.add_constant(df[["grain", "gdp", "birth_rate", "death_rate"]]) m4 = sm.OLS(y, X4).fit() print(m4.summary2().tables[1].round(4)) # 只打系数表，tidy 格式 print("F:", round(m4.fvalue, 3), "p:", m4.f_pvalue)

summary2().tables[1]输出的是一个 DataFrame，比直接 print(summary()) 更适合贴进报告或者做二次筛选。这一步要盯的是每个变量自己的 p 值，而不是整体 F 检验：整体显著只说明「至少有一个变量有用」，不说明每个都有用。

系数表里 GDP 的 p=0.531，死亡率的 p=0.275，两个都没通过 0.05 的显著性检验；出生率系数 -1460.798，p<0.001，粮食产量系数 1.045，p<0.001。GDP 的系数是负的（-0.005），这个符号和直觉完全相反——经济越发达，人口应该越多才对。这种符号反转是多重共线性的典型症状，不是数据错了，是变量之间打架。四元模型对 2014 年的预测值是 138182.24，误差 1400.24，比一元模型的 5384 好了不少，但这一步的好转有运气成分，不能当成模型更优的证据。

4.2 逐步回归的进出机制与代码复现

原文对逐步回归的定义写得很准：「有进有出」，每引入一个变量就对已引入的变量逐个检验，变不显著就剔除，保证最后的子集是最优回归子集。在 SPSS 里的操作路径是「分析 → 回归 → 线性 → 方法选步进」，准则用 F 概率，进入 0.05、剔除 0.10。

实际跑出来三步：第一步只进出生率，截距 165146.388，系数 -2687.949，t=-12.794；第二步进粮食产量，截距 108407.805，出生率系数变成 -1534.270，粮食产量系数 0.816；第三步进死亡率，截距 134678.878，出生率 -1366.202，粮食产量 1.012，死亡率 -5681.557，p=0.035。GDP 全程没被选进来。

def forward_stepwise(X, y, threshold_in=0.05): """只做前向选择的最小实现；SPSS 的步进法还包含后向剔除""" included = [] while True: excluded = [c for c in X.columns if c not in included] best_p, best_feat = threshold_in, None for feat in excluded: model = sm.OLS(y, sm.add_constant(X[included + [feat]])).fit() p = model.pvalues[feat] if p < best_p: best_p, best_feat = p, feat if best_feat is None: break included.append(best_feat) return included print(forward_stepwise(df[["grain", "gdp", "birth_rate", "death_rate"]], y))

threshold_in对应 SPSS 里的「进入 0.05」，调大它会让更多变量被放进来，调小则更严格。这个函数没有实现原文说的「后向剔除」，严格意义上的逐步回归还需要在每轮加入后对已入变量重新检验，把 p 值超过 0.10 的踢出去。写课程实验够用，写正式分析建议直接用mlxtend里的stepwise_selection或者 R 的step()。

逐步回归的最终方程是 y = 134678.878 + 1.012x1 - 1366.202x3 - 5681.557x4，2014 年预测值 138506.22，误差 172.22。这个误差是全文所有模型里最小的。但原文接下来的判断出了问题：它说「不用逐步回归的预测值更接近真实值」，这个结论和数字是矛盾的，我在第 5 章会细说为什么不能这么比。

4.3 多重共线性诊断：特征值、条件指数与 VIF

原文用了三个指标判断共线性：特征值、条件指数、方差膨胀因子。表 2-5 里特征值有 3 个接近 0（0.028、0.003、0.003），条件指数有 3 个大于 10（33.798、36.436、86.828），死亡率在相关系数矩阵中的数值接近 1。三个信号同时出现，结论就是存在中等偏严重的共线性。

from statsmodels.stats.outliers_influence import variance_inflation_factor Xc = sm.add_constant(df[["grain", "gdp", "birth_rate", "death_rate"]]) vif = pd.DataFrame({ "feature": Xc.columns, "VIF": [variance_inflation_factor(Xc.values, i) for i in range(Xc.shape[1])], }) print(vif.sort_values("VIF", ascending=False))

variance_inflation_factor的第二个参数是列索引，必须从 0 开始遍历整张矩阵，包括常数列。常数列的 VIF 是 1，不用管它，看其余四行即可。经验阈值是 VIF > 10 认为共线性严重，5 到 10 之间要警惕，小于 5 基本可以放心。

剔除 GDP 之后重跑，出生率 VIF=2.868、粮食产量 VIF=4.146、死亡率 VIF=1.808，全部小于 10，说明剩下的三个变量之间的相关程度可以接受。注意 VIF 的计算公式是 1/(1-R²_j)，其中 R²_j 是把第 j 个自变量当因变量、其余自变量当自变量回归得到的决定系数；原文提到的容许度就是 VIF 的倒数，两者说的是同一件事。这一步的正确做法是「剔除一个、重算一次、再看」，而不是一次剔掉所有 VIF 高的变量，因为变量之间会互相影响。

4.4 主成分降维的收益与代价

主成分分析在这份论文里的作用是降维 + 消共线性。原始四个变量做 PCA，前两个主成分的特征值分别是 3.040 和 0.727，解释方差比例 76.004% 和 18.176%，累计 94.180%，所以取前两个主成分建回归模型。

from sklearn.decomposition import PCA pca = PCA(n_components=2) scores = pca.fit_transform(StandardScaler().fit_transform(X)) print(pca.explained_variance_ratio_, pca.explained_variance_ratio_.cumsum()) print(pd.DataFrame(pca.components_, columns=X.columns))

explained_variance_ratio_给出每个主成分解释的方差占比，components_给出载荷矩阵，也就是每个主成分是原始变量的什么线性组合。载荷矩阵是把主成分翻译回业务语言的唯一入口，不看它就没法解释模型。

原文最终给出的方程是 y = 53075.711 + 1.468x1 + 2.494×10⁻⁵x2，2014 年预测 142180.53，误差 5398.53。这个结果几乎退化成一元回归了——粮食产量的系数 1.468 和一元模型一模一样，GDP 的系数小到可以忽略。也就是说，降维之后信息重新集中到了粮食产量这一个变量上，PCA 在这里没有带来新信息，只带来了「系数看不懂」的代价。

主成分回归的通用问题是可解释性：主成分是原始变量的线性组合，往往没有明确的业务含义，系数没法回答「出生率每上升一个千分点人口怎么变」这种问题。原文自己也承认了这一点，说主成分的解释「多少带点模糊性，是变量降维过程中不得不付出的代价」。用不用 PCA，取决于你要的是预测精度还是变量解释，两者在这类场景里很难兼得。

4.5 三次曲线与 Logistic：非线性分支的拟合代码

线性模型全跑完之后，原文换了个思路：不再找自变量，直接把时间 t 当自变量，拟合总人口随时间变化的趋势。以 1980 年为基准，t 取 1 到 35。散点图显示总人口大致符合三次函数形态，于是设定 y = β0 + β1t + β2t² + β3t³ + ε。

import numpy as np t = np.arange(1, 36) # 1980 对应 t=1，2014 对应 t=35 coef = np.polyfit(t, df["population"].values, deg=3) print(dict(zip(["t^3", "t^2", "t", "const"], coef.round(4)))) y_hat_2014 = np.polyval(coef, 35) print(f"2014 预测 {y_hat_2014:.2f}")

np.polyfit返回的系数是按幂次从高到低排列的，和 SPSS 表格里的顺序刚好相反，代回方程时别接错。deg=3对应三次多项式，调成 2 就是抛物线，调成 4 会开始过拟合。原文得到的方程是 y = 96103.588 + 1849.094t - 17.683t² - 0.070t³，R=1.000，R²=0.999，标准误 365.498，F=11749.732，p<0.001，2014 年预测 136174.95，误差 -607.05。

这个模型里有个细节容易被忽略：三次项的 p=0.371，没通过显著性检验。也就是说降一次、二次之后再加三次项，对拟合的贡献在统计上不显著。但 R² 已经到 0.999，看起来漂亮得不行。这种「整体极显著、个别项不显著」的组合，在外推时会暴露问题——三次项系数是负的，意味着曲线在足够远的地方会掉头向下。

Logistic 模型的定位是「考虑人口增长的有限性」，公式是 p(t) = pm / (1 + b·e^(-k(t-t0)))，其中 pm 是人口上限。原文用 MATLAB 的 fittype 做拟合：

x = [2014 2013 2012 2011 2010 2009 2008 2007 2006 2005 2004 2003 2002 2001 ... 2000 1999 1998 1997 1996 1995 1994 1993 1992 1991 1990 1989 1988 1987 ... 1986 1985 1984 1983 1982 1981 1980]; y = [136782 136072 135404 134735 134091 133450 132802 132129 131448 130756 ... 129988 129227 128453 127627 126743 125786 124761 123626 122389 121121 ... 119850 118517 117171 115823 114333 112704 111026 109300 107507 105851 ... 104357 103008 101654 100072 98705]; x = x'; y = y'; st_ = [2000 127627 1]; % 初值：上限、形状参数、斜率 ft_ = fittype('a/(1+b*exp(-k*(x-1980)))', ... 'dependent', {'y'}, 'independent', {'x'}, ... 'coefficients', {'a', 'b', 'k'}); cf_ = fit(x, y, ft_, 'Startpoint', st_)

st_里的三个数是初值，分别对应 a（上限）、b（形状）、k（斜率）。fittype第一个参数是模型表达式字符串，dependent和independent声明变量名，coefficients声明待估参数，这几个名字必须和表达式里出现的符号完全一致，大小写也算。

运行结果 a=1.217e5（置信区间 1.144e5 到 1.29e5），b=1.176e5（置信区间 -4.741e9 到 4.742e9），k=13.21（置信区间 -4.032e4 到 4.035e4）。b 和 k 的置信区间横跨了正负几个数量级的范围，这说明参数根本不可识别——模型没有收敛到一个有意义的解。把 1980 年代进去算一下就知道：y = 121700 / (1 + 117600) ≈ 1.03，和实际值 98705 差了五个数量级。原文的结论是「1995 年人口实际值就超过了预测上限，所以该模型不适合预测」，这个观察是对的，但原因不是模型不适合，是这次拟合本身失败了。

用 Python 给一组量级合理的初值，结果会正常得多：

from scipy.optimize import curve_fit years = df.index.values.astype(float) pobs = df["population"].values def logistic(x, a, b, k): return a / (1 + b * np.exp(-k * (x - 1980))) # a 取比当前人口大 20% 的上限，b 由 (a/y0 - 1) 反推，k 取 0.02 量级 a0 = pobs[0] * 1.2 b0 = a0 / pobs[0] - 1 p0 = [a0, b0, 0.02] popt, pcov = curve_fit(logistic, years, pobs, p0=p0, maxfev=20000) perr = np.sqrt(np.diag(pcov)) print("参数:", popt.round(4)) print("标准误:", perr.round(4))

p0的三个初值是有讲究的：a 必须大于样本中的最大值，b 由 a/y(1980)-1 反推，k 通常在 0.01 到 0.05 之间。maxfev是最大迭代次数，数据尺度大的时候默认值经常不够，报「 Optimal parameters not found」时先加这个。pcov是协方差矩阵，对角线开根号就是每个参数的标准误，判断拟合是否可信看这个比看曲线图快。

5. 避坑与排查：六个我在这份论文里踩过的坑

这份实验的价值有一半在模型，另一半在它暴露出来的问题。下面六条都是复现时会真实撞上的，按「现象 → 原因 → 解决」记下来，下次做同类分析可以直接对照。

5.1 数据口径层面的三个坑

出生率单位标错。现象：系数表里出生率的系数是 -1460.798，死亡率是 -4060.272，看起来像是「出生率每上升 1 个单位人口减少一千多万」，量级完全不合理。原因：表头把千分率标成了百分数，如果照着「%」去理解，就得把系数除以 100 才对应一个百分点的变化，但代码里用的还是原始千分率数值，两套口径混在一起。解决：在字段注释里写死单位，读入后立刻做一次assert df["birth_rate"].max() < 100，因为千分率不可能到三位数；所有报告里的系数解释都带上单位，写成「每千分点」而不是「每百分点」。

时间基准有三套。现象：三次曲线用的是 t=1..35，Logistic 用的是 (x-1980) 从 0 开始，一元和多元用的是原始年份 1980..2014，三个模型之间系数没法互相对照，外推到 2015 年时不知道该代 36 还是 35。原因：每个模型各自造了一套时间变量，没有统一入口。解决：表里只保留 year 一列，需要 t 的时候统一用t = year - 1980 + 1，写成函数而不是到处手写，Logistic 里也改用同一个 t 再做中心化。

GDP 是当年价。现象：GDP 系数为负（-0.005），和「经济越好人口越多」的直觉相反。原因：两个因素叠加，一是 GDP 和粮食产量、出生率之间存在强共线性，系数符号被带偏；二是当年价没有做价格平减，不同年份的数值不可比。解决：换成可比价序列，或者直接用 GDP 增长率代替绝对值；如果数据拿不到，就在报告里明确写「GDP 系数不可解释」，不要硬编业务含义。

5.2 模型与工具层面的三个坑

拿单点误差挑模型。现象：原文先用四元模型算出 2014 年误差 1400，再用逐步回归算出误差 172，最后写「不用逐步回归的预测值更接近真实值」，前后矛盾。原因：2014 这个点既是评估集又参与了模型选择，哪条路走下来误差小带有很大随机性，用单个点判断模型优劣本身就不成立。解决：把最后 5 年（2010–2014）留出来做滚动回测，每次用前面的数据拟合、预测下一年，把 5 个误差取平均；同时看 AIC/BIC，这两个指标对变量个数有惩罚，比单点误差稳。血泪经验是，任何在测试集上比较过的模型，都不能再用同一批测试集报最终指标。

Logistic 初值给错导致拟合失败。现象：MATLAB 的 fit 返回 b=1.176e5、k=13.21，置信区间横跨 ±4.7e9，曲线几乎是水平线。原因：初值[2000 127627 1]和模型形式不匹配，a 取 2000 但人口上限应该是十几万，k 取 1 但实际斜率是 0.02 量级，优化器一出发就掉进了平坦区域。解决：按数据反推初值，a 取样本最大值乘 1.2，b 用 a/y0-1 算，k 从 0.02 试起；拟合完先看参数标准误，标准误比参数本身还大就说明没收敛，换初值重来。

三次项不显著却保留。现象：三次曲线 R²=0.999，但 t³ 的系数 -0.070，p=0.371，没通过检验。原因：判断标准用的是 R² 而不是系数显著性，R² 会被高次项天然抬高。解决：报告里同时给 R² 和各项 p 值，别只报一个；外推前先算一下导数在预测区间内会不会变号，如果会，说明曲线在远处要掉头，这个预测就不能用。黑匣子式的「R² 高就是好模型」在多项式拟合里最容易翻车。

5.3 复现前的自检清单

把上面这些坑固化成几条断言，每次换数据集先跑一遍，能省掉大半回查时间：

assert df.shape == (35, 5), "行列数不对，先确认年份是否连续" assert df.index.is_monotonic_increasing, "年份必须升序，时间变量才成立" assert df.notna().all().all(), "有缺失值，先决定插补还是删行" assert (df[["birth_rate", "death_rate"]] < 100).all().all(), "出生/死亡率疑似单位错误" assert df["population"].is_monotonic_increasing, "总人口非单调，检查数据源或口径" assert (df[["grain", "gdp"]] > 0).all().all(), "产量/GDP 不该出现非正值"

六条断言对应六类最常出的问题：形状、排序、缺失、单位、单调性、非正值。它们不检查模型对不对，只检查数据能不能进模型，但恰恰是这一步出错最难回头查。后悔药是提前写好这些断言，而不是等模型跑出离谱结果再回头翻原始表。

6. 二次开发：把论文改成一个能反复跑的预测脚本

这份论文目前的状态是「一次性分析」：换一批数据，所有步骤都得手动重来一遍。要把它变成能复用的东西，改造方向有三个，按性价比排序。

滚动回测的实现不复杂，核心是「每次只用预测年份之前的数据拟合」：

def rolling_backtest(df, model_fn, start_year=2010): """model_fn(train_df) -> 预测函数；返回逐年误差列表""" errors = [] for year in range(start_year, df.index.max() + 1): train = df.loc[: year - 1] # 严格只用当年之前的数据 pred = model_fn(train)(year) # 由训练集产出的预测器预测该年 errors.append(pred - df.loc[year, "population"]) return errors def cubic_model(train): t = train.index.values - 1980 + 1 coef = np.polyfit(t, train["population"].values, 3) return lambda year: np.polyval(coef, year - 1980 + 1) print([round(e, 1) for e in rolling_backtest(df, cubic_model)])

start_year决定回测窗口，样本 35 条的话从 2010 开始能留下 5 个测试点，够看趋势了。model_fn接收训练集、返回一个「年份 → 预测值」的闭包，这个接口设计的好处是加新模型只需要再写一个函数，主流程一行不用改。注意三次多项式每次只用 31 个点拟合，系数会和全样本拟合出来的不一样，这是正常的，回测要的就是「每年重新估参」的过程。

另一个值得做的是换掉三次曲线这个主模型。三次多项式没有增长上限，外推到足够远的年份会掉头向下，物理上不合理。可以试灰色 GM(1,1)，它对短序列、单调序列很友好，代码不到二十行；也可以试带外生变量的 ARIMAX，把出生率、死亡率作为协变量放进去，比单纯用时间 t 更有解释力。这几种模型的评估都走同一个rolling_backtest，横向比较才有意义。

还有一个容易被忽略的技巧：把「模型输出」和「报告输出」分开。前者返回一个包含预测值、区间、参数的字典，后者负责渲染表格和文字。分开之后，换数据集只需要换输入文件，报告格式不用动；做敏感性分析时也能批量跑参数组合，不用每次手动改图。

从那以后我每次做时间序列类的模型对比，都强制留最后 5 年做滚动回测，哪怕样本只有三十几条也照留不误——用单点误差挑模型这件事，翻车概率比想象中高得多，而且翻车之后很难自证。这份论文里最值钱的部分不是哪个模型误差最小，而是它把「误差最小的那个可能是运气」这件事原原本本留在了纸面上。希望帮到你。

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