简介围绕中国人口增长预测这一数学建模问题这份文档提供了一套较完整的分析与建模方案适合数学建模参赛者、人口与统计方向学习者以及使用MATLAB做预测实验的读者参考。内容以Logistic阻滞增长模型和按年龄分布的Leslie模型为主线先比较不同起止年份数据对拟合效果的影响再用Leslie矩阵预测总人口、劳动年龄人口与老年人口并延伸到人口老龄化、抚养比、育龄妇女人数等指标最后对两个模型做残差分析、优缺点评价与推广讨论。压缩包仅含1个doc文件大小约988KB篇幅集中、结构清晰便于直接查阅建模假设、符号说明、公式推导与预测结论。已有5156人学习可作为课程论文、竞赛复盘或人口预测入门案例帮助读者理解如何把MATLAB计算与人口结构分析结合起来。1. 一条 0.9987 的曲线预测不了老龄化Leslie 模型到底补上了哪一维先给一个反直觉的结论用 1980—2005 年总人口数据拟合出来的 Logistic 曲线可决系数 0.9987残差几乎贴着零轴走但它对“老龄化”这三个字一个字都说不出来。原因很简单Logistic 把一国人压成了一个标量 N(t)所有年龄被平均掉了。Leslie 人口增长模型真正解决的问题不是把总数预测得更准而是把“年龄”抬升为状态向量的一维让生育率、存活率、性别比在同一个矩阵乘法里闭合。搞懂它你才能解释为什么同一批数据能算出 2023 年达峰、为什么总和生育率要从 1.42 硬拉到 1.8、为什么峰值年份会比官方预测提前十年。这几个问题恰恰是这类模型最容易被当成黑匣子的地方。2. Logistic 阻滞增长率为什么会自己踩刹车2.1 从指数增长到 r(N)r(1−N/K)一条线性假设如何决定 S 型指数增长模型只有一个假设增长率恒定dN/dt rN解出来是 N N₀e^{rt}。它在短期能蒙对长期一定翻车因为 r 恒等意味着资源无限。阻滞增长模型改的就是这一条把 r 从常数变成 N 的函数且必须是减函数——人越多人均资源越少增长率越低。最简单的假定是线性r(N) r − sN。再用一个边界条件把它钉死当 N 达到资源与环境能容纳的上限 K 时增长率归零即 r(K) 0代入得 s r/K于是r(N) r(1 − N/K)把它代回 dN/dt r(N)·N得到 Logistic 方程 dN/dt rN(1 − N/K)分离变量求解后是N(t) K / (1 (K/N₀ − 1)·e^{−rt})这条曲线有三个可验证的特征N 单调趋近 K拐点在 N K/2此处增长率最大t → ∞ 时 N → K。模型里所有的“人口学”都被压缩进了两个参数——K 和 r。这也是它的第一个软肋K 不是一个能独立观测到的物理量它只能从数据里反演出来数据段一变K 就变。注意这里的“线性阻滞”是一阶近似不是自然定律。选线性的理由是可解析、参数少、可解释而不是它更真。这一点决定了后面三段拟合会给出三个不同的 K。整份推导里最值得琢磨的设计取舍在这儿模型假设人口增长只受“总量拥挤”影响完全不区分 0 岁和 80 岁。这在总量预测上够用但一旦你要回答“劳动年龄人口什么时候见顶”“抚养比会不会翻倍”这个模型在结构上就答不了不是精度不够是维度不够。2.2 非线性最小二乘nlinfit 怎么把 K 和 r 拧出来Logistic 的解对参数是非线性的没法用普通最小二乘的闭式解只能用迭代。思路是把观测值 x_i 和模型值 f(t_i; K, r) 的差平方和作为目标函数用 Gauss-Newton 或 Levenberg-Marquardt 一步步往下走。这类算法对初值极其敏感初值给歪了要么不收敛要么收敛到一个毫无意义的局部极小。原文给初值的办法很土但有效拿两个年份的数据点联立方程先解出一组 (K, r) 当起点再交给 nlinfit 精修。t 0:51; % 以 1954 年为初始年 t0 x [60.2 61.5 62.8 64.6 66.0 67.2 66.2 65.9 67.3 69.1 ... 70.4 72.5 74.5 76.3 78.5 80.7 83.0 85.2 87.1 89.2 ... 90.9 92.4 93.7 95.0 96.259 97.5 98.705 100.1 101.654 ... 103.008 104.357 105.851 107.5 109.3 111.026 112.704 ... 114.333 115.823 117.171 118.517 119.85 121.121 122.389 ... 123.626 124.761 125.786 126.743 127.627 128.453 129.227 ... 129.988 130.756]; % 单位千万 % 用两个年份的点联立先解一组粗糙初值 [c, d] solve(c/(1(c/60.2-1)*exp(-5*d))67.2, ... c/(1(c/60.2-1)*exp(-20*d))90.9, c, d); b0 [241.9598, 0.02985]; % 迭代起点避免从零猜 fun inline(b(1)./(1(b(1)/60.2-1).*exp(-b(2).*t)),b,t); [b1, r1] nlinfit(t, x, fun, b0); % r1 是残差向量 y 180.9871 ./ (1 (180.9871/60.2 - 1) .* exp(-0.0336 .* t)); plot(t, x, *, t, y, -or); R 1 - sum(r1.^2) / sum((x - mean(x)).^2); % 可决系数 W sum(abs(r1)); % 残差绝对值之和几个参数要看清。b(1)就是 K拟合出 180.9871千万即 18.1 亿b(2)是 r0.0336含义是人均增长率在人口接近零时的理论上限。r1是 52 个残差R是可决系数W是残差绝对值之和。nlinfit返回的b1和手写在y里的两个数应当一致如果不一致说明迭代没收敛到你以为的那个点——这是排查拟合问题时第一个该看的信号。这里有个常被忽略的细节可决系数 R 1 − SSE/SST 的写法必须是两个求和之比不能用逐元素相除再求和。时间序列拟合里 R² 天然偏高因为趋势项主导了总离差0.99 以上并不稀奇它只能说明曲线抓住了趋势说明不了残差是否服从正态、参数是否可解释。2.3 三段起始年份对照0.9987 是拟合出来的还是选出来的原文做了一件很关键的事同一套 Logistic 形式分别用 1954 年、1963 年、1980 年作为初始年份各拟合一次比较结果。起始年份拟合窗口K千万r残差绝对值之和可决系数19541954—2005180.98710.033657.9992未给出19631963—2005151.45130.048427.8046未给出19801980—2005由 nlinfit 给出由 nlinfit 给出10.16990.9987三个 K 从 181 千万掉到 151 千万差了近 30 千万相当于三个亿。为什么因为最小二乘法有一个隐含前提残差独立同分布、近似正态。1954—1962 年这一段既包含第一次出生高峰形成的“由缓到快”增长又夹着一段数据本身的异常波动残差在这段出现系统性偏离正态假设不成立。早期数据在时间轴上杠杆最长它把 K 往上拽拽出一个偏大的饱和水平。1963 年之后增长基本按自然规律走1980 年之后生育政策稳定没有大的外部干扰残差大致关于零轴对称绝对值之和也最小。这就是“0.9987”的真正来源——不是模型更聪明而是数据段更符合模型的前提假设。参数可信度是被数据口径决定的不是被算法决定的。三条曲线对同一批年份给出的预测差得很明显单位千万年份1954 起始曲线1963 起始曲线1980 起始曲线2000126.7649126.3338126.47302021149.4251141.2856142.34892048166.7683148.5580150.27751954 起始那条曲线一路飙到 166.7 千万明显偏乐观1963 起始偏保守1980 起始居中和专家预测最接近。原文据此采信 1980 段的结果2010 年 13.55357 亿、2020 年 14.18440 亿、2033 年 14.70172 亿。这套“分段拟合 交叉比对”的做法本质上是在用数据口径的稳定性替代统计检验在没有更长可靠序列的情况下这是少有的可行解。3. Leslie 矩阵把出生和死亡压缩成一次矩阵乘法3.1 年龄离散化与状态向量为什么只数女性Logistic 的一个标量装不下年龄结构Leslie 换了个载体一个向量。把人口按年龄等间隔切成 m1 组原文按 1 岁一组0 到 90 岁90 岁以上合并时间也离散成同样步长于是“一年”就是一次向量更新。设 n_i(t) 是第 t 年第 i 年龄组的人口数组成列向量 n(t) (n₀, n₁, …, n_m)ᵀ。原文只追踪女性理由是出生数由育龄女性数量乘生育率决定男性数量在模型假设下不构成对出生数的直接约束性别比不消失而是被折算进生育率里。这个处理让矩阵规模直接减半代价是男性人口只能事后按比例反推一旦出生性别比随时间漂移反推就会有系统偏差。模型还假设生育率和死亡率不随时间变化。这条假设是整个模型最脆的地方——它等价于说“今年 25 岁的人怎么生、怎么死二十年后 25 岁的人也一模一样”。在十年尺度上勉强能接受拉到五十年尺度就是明显的硬伤。原文在模型评价里也承认了这一点改进方向是把死亡率做成时变的后面会讲。3.2 第一行生育率、次对角线存活率矩阵结构就是人口学假设Leslie 矩阵的结构本身就在陈述人口学。设 b_i 是第 i 年龄组每个女性平均生育女儿的人数s_i 是第 i 年龄组的存活率 1 − 死亡率则⎡ b₀ b₁ b₂ … b_{m-1} b_m ⎤ ⎢ s₀ 0 0 … 0 0 ⎥ L ⎢ 0 s₁ 0 … 0 0 ⎥ ⎢ 0 0 s₂ … 0 0 ⎥ ⎢ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⎥ ⎣ 0 0 0 … s_{m-1} 0 ⎦第一行全是生育率含义是“t 年各年龄组女性生出的女婴全部进入 t1 年的 0 岁组”——所以 0 岁组的人口是各行生育率的加权和权重就是各年龄组女性人数。次对角线全是存活率含义是“t 年 i 岁的人活到 t1 年就是 i1 岁”所以第 i1 行第 i 列放 s_i。其余位置为 0因为一个人不可能一年长两岁也不可能往回长。矩阵成为 Leslie 矩阵需要两个条件b_i ≥ 0且至少一个严格大于 0s_i ∈ (0, 1]。前一条保证有出生后一条保证死亡率不为 1否则整列清零那一代人直接断掉。注意第一行最后一个元素 b_m 对应 90 岁组现实中基本为 0但它留在那里表示“末组仍在生育”的可能性删掉也不影响结果。这里有一个容易看漏的设计矩阵里没有“死亡”这个动作死亡是以“不进入下一行”的方式隐式表达的。s_i 0.98 意味着 2% 的人凭空消失不会出现在任何一行里。这使得 Leslie 模型天然是一个只增不减的入口结构——总人口的变化完全由出生数和各年龄存活率决定。3.3 n(t1)L·n(t) 的迭代长期增长率藏在主特征值里有了矩阵递推只剩一行n(t1) L·n(t)展开就是 n(t) Lᵗ·n(0)。原文实现的核心循环非常短function [N, L] leslie_forecast(f, s, n0, T) % f : 1×(m1) 各年龄组女性平均生育女儿数第一行 % s : 1×m 各年龄组女性存活率次对角线 % n0 : (m1)×1 初始年份分年龄女性人口 % T : 预测年数 m numel(n0) - 1; L zeros(m 1, m 1); L(1, :) f; % 生育0 岁组由各行生育率加权得到 L(2:end, 1:m) diag(s); % 存活i 岁到 i1 岁 n n0; N zeros(1, T 1); N(1) sum(n); for t 1:T n L * n; % 一次乘法 一年的出生 存活 死亡 N(t1) sum(n); end endL(1,:) f决定了出生L(2:end,1:m) diag(s)把存活率放到次对角线上其余为 0。循环里没有任何分支判断因为所有的人口学规则都已经编码进矩阵结构了。sum(n)得到的是女性总数要还原总人口还得乘性别比系数——原文的预测值 2010 年 14.2609 亿、2020 年 14.9513 亿、2023 年峰值 14.985 亿都是这么还原出来的。真正有意思的是 L 的谱性质。L 是一个非负矩阵且结构上不可约每一代人最终都能通过生育链条连回 0 岁组满足 Perron-Frobenius 定理的适用条件因此存在唯一一个正的主特征值 λ₁对应一个正特征向量 v。长期来看 n(t) ≈ c·λ₁ᵗ·v也就是说λ₁ 1总人口长期指数增长增长率为 λ₁ − 1λ₁ 1总人口长期指数衰减v 归一化之后就是稳定年龄分布不管初始人口长什么样迭代足够多步之后形状都会收敛到 v。这解释了 Leslie 模型一个反直觉的表现即使总人口还在涨老年人口占比也可能在快速上升。因为老龄化不是“总数问题”是向量形状从 n(0) 向 v 收敛的过程加上生育率低于更替水平后 λ₁ 被压到 1 以下v 本身偏老。原文预测 2020 年 60 岁以上 2.907 亿、占比 19.443%65 岁以上 2.0628 亿、占比 13.797%到 40 年代中后期 60 岁以上达 4.45 亿、占比 33.277%本质上是这套向量迭代的自然结果而不是额外加的规则。4. 从抽样表到矩阵元素b 和 s 是怎么一项项算出来的4.1 市、镇、乡三层抽样如何还原分年龄女性人口Leslie 矩阵的两个关键输入——初始女性人口向量 n₀ 和生育率行 b——都不是现成的得从抽样调查数据里一步步折算。原文用的是 2001 年的抽样表按市、镇、乡三层分别给出男女样本量层级男女合计市147907147465295372镇8027977976158255乡394690372242766932合计6228765976831220559按这个样本合计约 122 万人算市、镇、乡的样本占比分别约为 24.2%、13.0%、62.8%。用这个比例乘上 2001 年全国总人口 127.627 千万就得到三层的总人口 P_市、P_镇、P_乡。再乘上各层分年龄女性的占比 q_{r,i}就还原出全国第 i 岁女性人数n_i Σ_{r∈{市,镇,乡}} P_r · q_{r,i}这一步的关键在于“分层”不是可选项。市、镇、乡的年龄结构、生育率、死亡率差异都很大如果直接把三层样本混在一起算平均占比会得到一个谁都不像的中间结构尤其是乡村生育率被城镇拉低预测出的出生数会系统性偏少。原文分三层各建一套参数再相加计算量大但结构误差小得多。4.2 女婴出生率的折算性别比 103~107 是怎么进矩阵的生育率 b_i 的定义不是“生多少个孩子”而是“生多少个女儿”这一步最容易算错。原文的流程是先算第 i 年龄组育龄妇女按年龄别生育率生育的孩子总数含男女再乘女婴占比最后除以该年龄组女性总数。女婴占比由出生人口性别比换算。性别比的定义是“以女婴为 100 时对应的男婴数”正常范围 103~107。若某层性别比为 SRB则女婴占比为 100/(100 SRB)。取 SRB 106女婴占比约 48.5%。% 第 i 年龄组从“生孩子”到“生女儿”再到“每个女性生几个女儿” births_i n_f_i .* fert_i; % 含男女的出生数 share_girl 100 ./ (100 srb_i); % 女婴占比srb_i 为实测性别比 b_i births_i .* share_girl ./ n_f_i; % 每个女性平均生育女儿数 % 总和生育率各年龄组生育率之和含男女口径 TFR_raw sum(fert_i);这里有个细节值得说清楚。fert_i是年龄别生育率口径包含男女b_i是“女儿数”口径只含女。两者差一个share_girl因子差得不多约 3%但会一路累积到五十年后的总人口上。原文在数据描述里专门强调了“此处应该是包含男孩和女孩”说明这一步在当时是容易混淆的地方。另一点srb_i用的是分层的实测值而不是统一取 106因为市、镇、乡的出生性别比差异明显统一取值会抹掉这个结构。4.3 总和生育率从 1.42 修正到 1.8一次显式的数据校准整套建模里最有争议、也最能体现“工程判断”的一步在这里。把上面算出的各年龄组生育率直接求和得到的总和生育率小于 1.8原文说“误差很大”进一步算出的人口更替率只有 1.42五年均值而当时中外专家给出的更替水平是 1.8。差了 0.38这不是随机误差是系统偏差。抽样调查里高年龄、多孩的生育行为存在明显漏报越是不合规的生育越不容易被统计到于是生育率被系统性压低。如果直接用 1.42 建矩阵λ₁ 会被压得很低总人口会在远早于现实的年份开始下降。原文的处理是显式修正把生育率整体缩放到总和生育率等于 1.8 的水平用这个校准后的 b 重建 Leslie 矩阵。修正后预测出的峰值年份比专家报告提前了约 10 年原文把原因归结为两点一是出生性别比中女性占比偏小二是原始数据本身不够真实。这一步的设计动机值得单独拎出来讲。用外部锚点覆盖数据算出来的值在统计上是不“干净”的但它换来了两样东西一是可解释性1.8 是有明确人口学含义的更替水平读者能判断这个数靠不靠谱二是可比性不同年份、不同来源的模型可以用同一把尺子量。代价是预测结果的绝对水平被锚点决定了如果 1.8 本身选错整个五十年预测都会跟着偏。这也解释了原文为什么要同时保留“原始数据矩阵”和“1.8 矩阵”两套结果——不是二选一是给读者一个区间。5. 避坑与排查人口预测最容易翻车的五处5.1 起始年份选错曲线会“假装”拟合得很好现象同一套 Logistic 形式换一个起始年份K 从 181 千万跳到 151 千万2048 年的预测值差了 18 千万但两条曲线的可决系数都在 0.99 以上。原因最小二乘的统计基础是残差独立同分布、近似正态。序列前段如果有出生高峰或数据异常残差在那里出现系统性偏离而早期数据在时间轴上杠杆最长会把 K 往上拽。高 R² 只说明曲线抓住了趋势不说明残差结构合格。解决按制度口径或数据口径分段拟合画残差散点图看残差是否关于零轴对称、绝对值之和是否收敛。原文三段对照里1980 段的残差绝对值之和是 10.16991954 段是 57.9992差了五倍多这个差距比 R² 更能说明问题。5.2 生育率没折算性别比总和生育率系统性偏低现象算出来的总和生育率总比外部数据低一截或者女婴数量比实际多出 3% 左右。原因年龄别生育率的口径含男女直接当“女儿数”用会高估女婴出生反过来如果只算女婴却忘了除以该年龄组女性总数b_i 的量纲就错了。解决严格分两步走——先算含男女的出生数再乘 100/(100SRB) 折算成女婴最后除以该组女性人数。性别比用分层实测值不要统一取 106。5.3 90 岁截断与末组处理存活率对角线的边界现象长期预测里总人口偏低且时间越长偏得越多但短期几年看不出来。原因矩阵最后一行只有一个 s_{m-1}末组90 岁以上没有对应的存活率迭代一次之后这批人全部消失。虽然 90 岁以上人口占比很小但它是唯一一个“只出不进”的组每年清零会持续累积。解决把末组当成开区间处理给最后一列补一个接近 0 但不为 0 的存活率或者把末组人口按固定比例并入 m−1 组。检查方法是单独追踪末组人数的时间序列如果它每年归零就是这里漏了。5.4 峰值年份提前 10 年数据口径不一致的锅现象Leslie 模型预测 2023 年达到峰值 14.985 亿而官方报告说峰值在 2033 年前后差了整整十年。原因不是矩阵算错了是输入参数偏了。原始数据算出的更替率只有 1.42远低于专家给出的 1.8同时出生性别比中女性占比偏小两者叠加把 λ₁ 压得更低达峰自然提前。解决把生育率缩放校准到 1.8 重建矩阵并跑一组 k 1.6 / 1.8 / 2.0 / 2.2 的包络看峰值年份和峰值高度的分布区间。单点预测没有意义区间才有。5.5 用可决系数判断 Leslie 模型好坏指标用错了地方现象想给 Leslie 模型也算一个 R²算出来要么很低要么根本对不上。原因R² 衡量的是模型对历史序列的拟合优度前提是参数通过最小二乘从历史数据反演。Leslie 是递推结构模型b 和 s 是从生育率、死亡率直接算出来的不是拟合出来的没有“拟合值 vs 观测值”的对应关系。解决Leslie 的验证走另外三条路——回测用早期数据建模型预测已知年份比误差、结构检验各年龄组人数是否单调、是否出现负值、灵敏度分析扫 k 看结果怎么变。拿 R² 去评价它属于用错尺子。6. 灵敏度与取舍k1.8 不是答案是一个平衡点6.1 四组总和生育率下的分叉负增长与膨胀之间把总和生育率 k 设成 1.6、1.8、2.0、2.2 四档其余参数固定重跑 Leslie 模型会看到两条曲线朝相反方向分叉。k 1.6 时人口增长最慢达峰后下降速度最快出现明显的负增长同时老龄化指数上升最陡年龄结构迅速老化。k 2.2 时人口增长快、峰值高、达峰后下降平缓老龄化被稀释但总人口会持续膨胀对资源环境形成压力。k 1.8 附近总人口曲线相对平稳老龄化速度处于中间位置。k总和生育率总人口走势老龄化速度结构稳定性1.6早达峰、快速负增长最快差1.8平稳达峰、缓降中等较好2.0峰值偏高、下降慢较慢一般2.2持续膨胀最慢差资源压力这张表说明的是一件很朴素的事控制总人口和控制老龄化是两个方向相反的力k 就是那个杠杆。只看总数会低估老龄化的代价只看老龄化会放任总量膨胀。原文说“只有掌握好一个平衡点”这句话在数学上就是 λ₁ 接近 1 且 v 的形状不过度偏老——两个条件同时满足的 k 区间很窄。6.2 把常数矩阵改成 L(k)生育模式 h_i 才是真正的控制变量固定矩阵的 Leslie 模型有一个理论缺陷b_i 不随时间变意味着生育模式五十年不变。现实中结婚年龄、生育间隔、孩次结构都在变。改进方向是把生育率行拆成“水平 × 形状”b_i k · h_i其中 h_i 是归一化的生育模式Σh_i 1k 是平均生育胎数。% 不直接改 b_i而是改生育模式 h_i 和水平 k h fert_shape / sum(fert_shape); % 归一化Σh 1 L zeros(m 1, m 1); L(1, :) k * h; % 只动 k形状由 h 决定 L(2:end, 1:m) diag(s);这样 k 和政策直接挂钩——调整结婚年龄、拉开两胎间隔改的是 h 的形状不是简单地把 k 调小。形状变了稳定年龄分布 v 跟着变老龄化路径也就变了这是固定矩阵做不到的。原文在模型改进里提的“引入控制变量、控制结婚年龄和生育间隔”落到底就是这句话。6.3 验证一个人口模型我一般看哪几件事回测排第一。拿 1990 年的数据建模型预测 2000 年和实际值比误差在 2% 以内才算过关。回测比任何拟合指标都诚实因为它用的是模型没见过的时间段。第二看残差结构不只看大小。残差如果全在零轴一侧说明系统性偏差如果在零轴上下对称分布才说明模型形式没选错。原文三段拟合的残差图价值就在这里。第三看参数敏感性。k 上下浮动 0.2峰值年份和峰值高度各变多少如果变动幅度超过你愿意接受的决策误差这个模型就不该拿来定政策。第四看独立来源。同一个指标用完全不同的数据、不同的模型算出来量级对不对得上。对不上不一定是自己错了但一定要能解释为什么。这几条我在早期做一个模拟项目时吃过亏当时只看 R²0.997汇报时很自信结果换一段数据重跑K 直接掉了 15%全场沉默。从那以后我每次碰人口类的递推模型都强制走一遍回测加 k 扫描哪怕结论看起来很明显。希望帮到你。本文还有配套的精品资源点击获取