简介：《Leslie人口增长模型》是一份面向人口统计、数学建模学习者与相关竞赛读者的doc文档，围绕中国人口总量及年龄结构的中长期预测展开。包内仅1个doc文件，压缩包约988KB，篇幅集中，便于直接查阅完整建模过程与结论。资源已有5157人学习，说明其在人口预测与MATLAB建模场景中具有一定参考价值。文档以Logistic阻滞增长模型和Leslie矩阵模型为主线，给出参数设定、数据选取、拟合与预测结果，并延伸到老龄化、人口抚养比、育龄妇女人数及劳动力结构分析，还包含残差检验、模型优缺点评价与推广讨论。读者可据此理解两类模型的适用条件、预测口径和结果解释方法，并借鉴其用MATLAB完成数据处理、矩阵构造与趋势分析的完整思路。

1. 一份跑不通的人口预测代码包，卡在哪几行

把一份人口预测的建模代码包解压出来，最常见的翻车不是公式看不懂，而是 MATLAB 一跑就报错、拟合曲线明明穿过了每个数据点、往外推三十年却给出一个人口腰斩的结果。这份 Leslie人口增长模型 资源把两条路线塞进了同一套 MATLAB软件 脚本：一条是只盯总量的 Logistic人口模型，一条是按年龄分组的 Leslie人口模型，中间还夹着老龄化、抚养比、育龄妇女这几组派生指标。适合正在做人口增长预测课程设计、建模复现，或者需要一套能改参数跑灵敏度的人口推算骨架的人。下面按「总量拟合 → 年龄结构 → 指标计算 → 排错 → 灵敏度」的顺序拆，重点放在参数怎么设、报错怎么看、结论什么时候不能信。

2. Logistic 阻滞增长模型：三个起点与 0.9987 的可决系数

2.1 为什么全段一次拟合会失真

Logistic 人口模型的形式不复杂：增长率随人口数量线性下降，解出来是一条 S 形曲线

N(t) = Nm / (1 + (Nm/N0 - 1) · exp(-r·t))

Nm 是资源环境能容纳的人口上限，r 是内禀增长率，N0 是初始年份人口。很多人拿到 1954 到 2005 年那 52 个总数点，直接一次性丢进拟合函数，跑出来的曲线看着挺顺，可决系数也不低，但一推到 2033 年就明显偏高。原因不在算法，在数据本身分了三段完全不同的生成机制。

1954 到 1962 年这段，先是出生高峰把人口规模从缓推快，接着连续几年人口增长几乎停滞，这一段残差既不服从正态分布也不对称，最小二乘的前提直接不成立。1963 到 1979 年，增长基本按自然规律走，城市因为收入水平生育率低一些，但整体没有大的外部扰动，这一段残差可以认为是正态的。1980 年之后生育水平进入受控阶段，机制又变了，不过干预强度本身稳定，残差同样近似正态。

所以这份资源的做法是把同一套 Logistic 形式拟合三次，起点分别是 1954、1963、1980，终点统一取 2005。三段拟合出来的参数差异很大，这不是模型不稳定，而是数据段本身的含义不同：

三个 Nm 差出近 30 千万，直接决定了 2033 年预测值落在 147 千万还是 158 千万。这也是为什么这份资源最终选了 1980 段：不是因为它的可决系数最高（0.9987 这个数看着漂亮，但另外两段也不会太低），而是它的残差结构最干净，可决系数才有解释力。

2.2 用三点反解给 nlinfit 一个靠谱初值

Logistic 拟合的第二个坑是初值。nlinfit对初值极其敏感，初值给偏了，要么报「Jacobian 矩阵奇异」，要么收敛到一个物理上没意义的解（比如 Nm 小于最大观测值）。原文附录里用符号解方程先估一组初始参数，思路是对的，但那段solve在新版 MATLAB 里语法已经变了，直接跑会挂。

我一般会用三个后期预测点反解 (Nm, r) 来定初值，比盲猜稳得多。原理是把 Logistic 解写成 Nm/N(t) - 1 = A·exp(-r·t)，两个时刻相除就能消掉 A，得到只含 Nm 的方程：

% 用三个已知预测点反解 Nm 与 r，仅用于给 nlinfit 提供初值 N1 = 135.5357; N2 = 141.8440; N3 = 147.0172; % t = 30, 40, 53 对应的预测值（千万） % 由 Nm/N1-1 与 Nm/N2-1 之比得 exp(10r)，再由 N2、N3 得 exp(13r)，联立消 r g = @(Nm) ((Nm/N1-1)./(Nm/N2-1)).^(13/10) - (Nm/N2-1)./(Nm/N3-1); Nm0 = fzero(g, [147.1, 300]); % 区间端点需异号，147.1 略大于 N3 r0 = log((Nm0/N1-1)/(Nm0/N2-1)) / 10; % 反解增长率 fprintf('初值 Nm=%.4f, r=%.5f\n', Nm0, r0);

fzero的区间不能乱给：左端点必须略大于最大预测值 N3，否则Nm/N3-1为负，开 1.3 次方会出现复数；右端点给到 300 千万足够，再大人群上限就荒唐了。跑出来 Nm0 落在 153.6 附近、r0 约 0.0478，把它俩填进拟合的初值，收敛基本一次到位。

有了初值再拟合 1980 段，代码可以完全不依赖 Curve Fitting Toolbox，只用nlinfit：

t = (0:25)'; % 1980 年为 t=0 N = [98.705 100.100 101.654 103.008 104.357 105.851 107.300 109.300 ... 111.026 112.704 114.333 115.823 117.171 118.517 119.850 121.121 ... 122.389 123.626 124.761 125.786 126.743 127.627 128.453 129.227 ... 129.988 130.756]'; % 单位：千万 N0 = N(1); logistic = @(p,t) p(1) ./ (1 + (p(1)/N0 - 1) .* exp(-p(2).*t)); p0 = [153.6, 0.0478]; % 上一步反解的初值 p = nlinfit(t, N, logistic, p0); % p(1)=Nm, p(2)=r fprintf('Nm=%.4f r=%.5f\n', p(1), p(2));

函数句柄logistic里的 N0 是闭包捕获的，不要写成p(3)再让nlinfit去估——初值年份人口是确定量，多估一个参数只会让雅可比矩阵条件数变差。数据向量 N 必须是列向量，写成行向量nlinfit会报维度不匹配，这个错误新手一天能撞三次。

2.3 可决系数和残差绝对值之和要一起看

拟合完别只看 R²。原文那段R1=r1.^2; R2=(x-mean(x)).^2; R=1-R1/R2是有问题的：r1是残差向量，r1.^2是逐元素平方后得到一个向量，再拿向量去除标量，MATLAB 会做逐元素除法，最后R也是个向量，不是标量。正确写法是先求和：

yhat = logistic(p, t); SSE = sum((N - yhat).^2); % 残差平方和 SST = sum((N - mean(N)).^2); % 总离差平方和 R2 = 1 - SSE/SST; % 可决系数，标量 W = sum(abs(N - yhat)); % 残差绝对值之和 fprintf('R2=%.4f W=%.4f\n', R2, W);

两个指标要配对使用。R² 靠近 1 只说明曲线解释了大部分方差，不说明残差结构合理；W 是绝对偏差的累加，对系统性偏移更敏感。1980 段拟合出来 R² 在 0.9987 量级、W 在 10.17 量级，对应到 26 个点上平均每个点偏 0.39 千万，也就是 390 万人，这个量级对全国总人口预测是可以接受的。而 1954 段 W 到了 58，平均每点偏 1.1 千万，即便 R² 也接近 0.99，这条曲线也不能用。

预测时把年份换算成相对 t，别直接丢年份进去：

ty = [2010 2020 2033]' - 1980; Npred = logistic(p, ty); % 单位：千万 disp(Npred); % 约 135.5 / 141.8 / 147.0

3. Leslie 矩阵：91 维递推里的人口年龄结构

3.1 从抽样男女数摊出分年龄女性向量

Logistic人口模型 只回答「总共有多少人」，回答不了「有多少人 65 岁以上」。Leslie人口模型 把人口按年龄切成一岁一组，只跟踪女性，靠生育率和存活率往前推。初始化这一步最容易错，因为原始数据给的是抽样调查的市、镇、乡男女合计数，不是现成的分年龄女性向量。

第一步先算市镇乡的人口占比。假设抽样里市男 147907、市女 147465、镇男 80279、镇女 77976、乡男 394690、乡女 372242，那么：

urban = 147907 + 147465; town = 80279 + 77976; rural = 394690 + 372242; tot = urban + town + rural; w = [urban town rural] / tot; % 三类地区人口权重 N2001 = 127.627; % 2001 年全国总人口（千万） N_region = w * N2001; % 市/镇/乡总人口（千万）

第二步把每个地区的总人口，按该地区分年龄女性的占比摊成人数。这里的关键是分母：占比是「该年龄女性占该地区总人口的比例」，不是「占该地区女性的比例」，两者差一个性别比，混用会让总和生育率凭空偏大一截。

% femaleRatio(i,r)：第 r 类地区、第 i 年龄组女性占该地区总人口比例（91×3） % 下标 i 从 1 开始对应 0 岁 x0 = zeros(91,1); for r = 1:3 x0 = x0 + N_region(r) * femaleRatio(:,r); end % x0 即 2001 年分年龄女性人口向量，单位千万

超过 90 岁的统一并入第 91 组，这是原文假设里写明的。不合并的话向量长度会随最高寿命浮动，Leslie 矩阵维度不固定，循环推到 2050 年时会报矩阵乘法维度不一致。

3.2 生育率、存活率和那条次对角线

Leslie 矩阵的结构是固定的：第一行全是分年龄女性生育率 b(i)，表示第 i 组每个女性平均生育的女儿数；次对角线是存活率 s(i)，表示第 i 组活到下一年的比例；其余元素全零。

n = 91; L = zeros(n); L(1,:) = b; % 第一行：生育率 for i = 2:n L(i,i-1) = s(i-1); % 次对角线：存活率 end

s(i) = 1 - d(i)，d 是分年龄女性死亡率。生育率 b 要经过两次换算才能直接用，这两步是整份资源里最值得抄的部分。

第一次换算是把「包含男孩的生育率」变成「只生女孩的生育率」。原始数据给的是育龄妇女总生育率，要除以出生性别比：

% sr：以女婴为 100 计的男婴数，正常区间 103~107 b_girl = b_total ./ (1 + sr/100);

忽略这一步，等于假设生出来的全是女孩，一代人之后女性人口会翻倍，推到 2050 年总人口能奔着 20 亿去。

第二次换算是总和生育率校准。把 b_girl 全加起来就是模型隐含的总和生育率 TFR，这份资源里算出来只有 1.42，而专家给的更替水平是 1.8。差这么多不是算错了，是抽样数据本身的生育率偏低。处理办法是按比例整体放大：

TFR_raw = sum(b_girl); scale = 1.8 / TFR_raw; % 约 1.268 b_fix = b_girl * scale; L(1,:) = b_fix; % 用校准后的生育率重建矩阵

按比例放大是个务实选择，它保留了年龄别生育率的相对形状，只调整水平。更讲究的做法是按年龄组分别加权，但抽样数据支撑不了那么细的校正。要注意放大之后矩阵最大特征值会变，人口从增长转为下降的拐点年份会整体后移，这正是后面灵敏度分析要扫的东西。

3.3 一年一步推到 2050 年

矩阵搭好之后递推本身只有三行：x(k+1) = L · x(k)。写成一个预分配的三维数组，方便后面按年龄切片算指标。

T = 50; % 2001 → 2051 X = zeros(n, T+1); X(:,1) = x0; for k = 1:T X(:,k+1) = L * X(:,k); end

循环里不要写X = [X, L*X(:,end)]，每次迭代都重新分配内存，50 步看不出来，改成 200 步做压力测试时能慢一个数量级。预分配之后再算总人口：

female_ratio = 0.4875; % 女性占总人口比例，由抽样数据算 total = sum(X,1) / female_ratio; % 1×(T+1)，单位千万

这里用固定性别比是个简化。真实性别比随年龄变化，0 岁组和 80 岁组差得远，用常数折算会让老年人口被低估几个百分点。要精确一点就按年龄组分别除，代价是得多准备一张分年龄性别比表。

按这套流程跑出来的结果是：总人口 2010 年 142.609 千万、2020 年 149.513 千万、2023 年达到峰值 149.85 千万，2033 年回落到 147.455 千万。劳动年龄人口（15 到 64 岁）2010 年 104.421 千万，2013 年到峰 104.852 千万，之后缓慢下行。峰值年份比参考预测提前了大约十年，原因后面第 5 章会拆。

4. 指标层：老龄化、抚养比与育龄妇女的算法

4.1 60 岁和 65 岁两条老龄化曲线差在哪

拿到 X 这个 91×(T+1) 的矩阵之后，老龄化就是切片求和。但 60 岁和 65 岁两个口径经常被混着讲，导致同一份预测能得出「老龄化严重」和「老龄化可控」两个结论。

age = (0:90)'; % 行下标对应的年龄 idx60 = age >= 60; idx65 = age >= 65; old60 = sum(X(idx60,:), 1) / female_ratio; % 60 岁及以上总人口 old65 = sum(X(idx65,:), 1) / female_ratio; % 65 岁及以上总人口 ratio60 = old60 ./ total; ratio65 = old65 ./ total;

判定老龄化的国际惯例是 65 岁以上占比超 7%，或者 60 岁以上占比超 10%，两条线不是一回事。这份资源算出来 2001 年 60 岁以上 15.538 千万、占比 11.5693%，65 岁以上占比从 2000 年的 8.009% 起步。到 2020 年，60 岁以上到 29.07 千万、占比 19.443%，65 岁以上到 20.628 千万、占比 13.797%。到本世纪 40 年代中后期形成高峰平台，60 岁以上 44.5 千万、占比 33.277%，65 岁以上 35.1 千万、占比 25.53%。

两条曲线的形状不一样。60 岁那条上升更陡，因为它的分母里包含了大量即将进入老年的 55 到 59 岁人群，这波人一旦跨过线，曲线会有个明显的抬升台阶。65 岁那条更平滑，但它在 2040 年之后的平台期更长。对外汇报时用哪条口径必须写清楚，否则两个团队拿同一份预测能吵起来。

4.2 抚养比的分母口径决定了结论方向

抚养比的定义是「非劳动年龄人口 / 劳动年龄人口」。看着简单，实际有两个口径：劳动年龄算 15 到 59 岁，还是 15 到 64 岁。这两个口径算出来的抚养比能差十几个百分点。

child = sum(X(1:15,:), 1); % 0-14 岁 laborA = sum(X(16:60,:), 1); % 15-59 岁 laborB = sum(X(16:65,:), 1); % 15-64 岁 aged = sum(X(66:end,:), 1); % 65 岁及以上 depA = (child + aged) ./ laborA; % 15-59 口径 depB = (child + aged) ./ laborB; % 15-64 口径

X的行下标是年龄加 1，所以 15 岁对应第 16 行，60 岁对应第 61 行。切片写成X(16:60,:)而不是X(15:59,:)，这个偏移在指标计算里极其常见，错了之后抚养比整体偏低一岁的人，误差不大但很难查。

按 15 到 64 岁口径算出来，总抚养比在 2001 年之后先降后升，2010 年前后到谷底，之后持续爬升到 2040 年前后接近 0.6。这个「先降后升」的形状本身是有信息量的：下降段对应劳动年龄人口占比上升的人口红利期，上升段是红利消退。如果只画一条上升的曲线，说明分母口径选错了，或者劳动年龄人口切片切错了。

提示：抚养比的分母如果是 15 到 64 岁，分子里的老年人口也必须用 65 岁及以上，两边口径要配套。分子用 60 岁以上、分母用 15 到 64 岁，是常见的口径错配。

4.3 育龄妇女和生育旺盛期妇女的峰值错位

控制人口这件事，光看总人口没用，得看未来还有多少人生孩子。这份资源算了两组：

fertile = sum(X(16:50,:), 1); % 育龄妇女 15-49 岁 peakAge = sum(X(21:30,:), 1); % 生育旺盛期 20-29 岁

育龄妇女（15 到 49 岁）人数在 2010 年左右到高峰，之后缓慢下降。生育旺盛期妇女（20 到 29 岁）在 2012 年达到高峰，2025 年前后进入一个小低谷，2037 年前后又出现一个小高峰。

这两个峰值的错位很有讲究。2037 年那个小高峰不是新出生人口变多了，而是 2012 年出生高峰那批女婴到 2037 年正好进入 25 岁。也就是说，一个人口高峰会在 25 年后自动复制出下一个生育高峰，这是 Leslie 模型里最直观的滞后效应。超过 49 岁的不计入育龄，所以这条曲线只受 15 到 49 岁这个窗口内的人口存量影响，和总人口曲线完全不同步。

如果育龄妇女曲线和总人口曲线同步下降、没有滞后错位，八成是年龄切片又切偏了，或者生育率向量的非零区间没对准 15 到 49 岁。

5. 常见问题排查：五次翻车的现象、成因与补丁

5.1 拟合曲线完美穿过数据，预测值却离谱

现象：nlinfit正常返回，R² 到 0.999，画图看着严丝合缝，但把 t 推到 53 年（2033 年）时总人口跳到 158 千万以上。

原因：nlinfit收敛到了局部最优，Nm 被估到一个偏大的值。Logistic 曲线在数据区间内对 Nm 不敏感，Nm 从 150 变到 180，区间内的曲线形状几乎不变，但外推几十年差距就拉开了。1990 年以前的拟合区间短、S 形曲线的拐点还没出现，Nm 的可辨识性最差。

解决：拟合完必须做一次外推检查。把预测值和独立来源的参考值（比如附件 1 里给的 2010、2020 预测数）对一下，偏差超过 2% 就换初值重拟合，或者改用 2.2 节的三点反解定初值。另外把拟合区间截短做交叉验证：用 1980 到 2000 年拟合，看能不能预测准 2001 到 2005 年，这个检验比 R² 有用得多。

5.2 报「矩阵维度必须一致」

现象：搭 Leslie 矩阵那行L(1,:) = b;报维度错误，或者L(i,i-1) = s(i-1)报下标越界。

原因：生育率向量 b 和存活率向量 s 的长度对不上。原始数据里生育率只给到 49 岁，死亡率给到 90 岁，直接拿来用就是一个 50 维向量往 91 列的矩阵第一行里塞。另一边s(i-1)在 i=91 时取 s(90)，如果 s 只有 90 个元素就会越界。

解决：进矩阵之前先补齐长度，缺的补零：

b = zeros(1, n); % 先建满长度 b(16:50) = b_raw; % 只填 15-49 岁 s = ones(1, n); % 默认存活率 1 s(1:numel(s_raw)) = s_raw; % 填入实际死亡率换算结果

补零对生育率是合理的（50 岁以上不生育），补 1 对存活率也合理（90 岁以上已并入最后一组）。但补完之后要再检查一次sum(b)的总和生育率有没有被补零破坏。

5.3 总和生育率怎么算都是 1.4 出头

现象：sum(b_girl)反复算是 1.42 左右，和预期的 1.8 对不上，怀疑数据读错了。

原因：不是读错，是抽样数据本身的生育水平就偏低。这份资源里明确写了这一点，中外专家对 90 年代中期以来的更替率估计是 1.8，抽样算出来只有 1.42，中间差的是漏报和抽样偏差。

解决：按比例校准，但校准之后一定要重新算一遍峰值年份。校准系数 1.8/1.42 ≈ 1.268，整体放大后人口峰值年份会从 2023 往后推，推进多少取决于死亡率结构。这份资源没做校准前的对照，实际操作时建议把校准前后的两条总人口曲线画在一张图上，峰值年份的偏移量就是校准带来的不确定性区间。

5.4 峰值年份比参考预测早了十年

现象：Leslie 模型跑出来总人口在 2023 年见顶，而附件里的参考预测说峰值在 2033 年前后。

原因：两个因素叠加。一是出生性别比里女性占比偏小，同样的总生育率下生出来的女孩更少，下一代育龄妇女基数被压小；二是 5.3 里提到的总和生育率偏低，虽然做了比例校准，但校准的是水平不是年龄结构。

解决：先分别关掉两个因素各跑一次，看谁的影响大。把出生性别比设成 103（正常下限）跑一遍，再把 TFR 设成 1.8 跑一遍，两次的峰值年份差就是各自的贡献。这个分解过程比直接给一个提前十年的结论有价值得多，因为它告诉你在哪一环上数据最不可靠。原文其实已经点到了：用 2001 年数据算出的更替率是 1.42，而专家口径是 1.8，「两者相差甚远」，说明原始数据不够真实。这句话是整个模型里最该被引用的一句。

5.5 残差图画出来是一条斜线

现象：残差散点图不是围绕横轴上下波动，而是从左上到右下一条斜线，或者呈现明显的 U 形。

原因：拟合区间跨越了机制不同的两段数据。用 1954 到 2005 全段拟合，1959 到 1961 年那三个点会形成巨大的负残差，后面的点全是正残差，图上看就是一条斜线。这不叫「拟合还行」，这叫模型形式和数据不匹配。

解决：按 2.1 节分三段拟合，分别看残差图。判断标准很简单：残差应该围绕零线随机波动，正负交替，不出现连续五个以上同号。1980 段拟合出来的残差绝对值之和 10.17，三个阶段的比值大约是 58 : 28 : 10，这个量级差就是分段必要性的直接证据。

6. 灵敏度分析：把 k 从 1.6 扫到 2.2 看什么

模型跑通只是第一步，真正决定这套代码能不能用来支撑判断的，是参数扫一遍之后结论稳不稳。Leslie 模型里最值得扫的参数就是总和生育率 k，因为它同时影响总人口、老龄化、抚养比三条曲线，而且方向相反。

ks = [1.6 1.8 2.0 2.2]; res = zeros(numel(ks), T+1); for j = 1:numel(ks) Lj = L; % 复制基础矩阵 Lj(1,:) = b_girl * (ks(j) / TFR_raw); % 按目标 k 重标生育率 Xj = zeros(n, T+1); Xj(:,1) = x0; for k = 1:T Xj(:,k+1) = Lj * Xj(:,k); end res(j,:) = sum(Xj,1) / female_ratio; % 存总人口曲线 end plot(2001:2051, res', 'LineWidth', 1.2); legend('k=1.6','k=1.8','k=2.0','k=2.2','Location','northwest');

跑出来的四条曲线形状差异很大。k=1.6 时人口增长缓慢，见顶之后下降很快，出现明显的负增长；k=2.2 时增长速度很快，见顶之后下降平缓，总人口在高位停留很久；只有 k 在 1.8 附近时曲线才相对稳定。这个结论和 Logistic 模型那边的 Nm 敏感性

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