MATLAB 阻滞增长模型

人口增长模型

Malthus模型（指数模型）

最早的人口增长模型是 Malthus 于 1798 年提出的指数模型，基本假设是人口增长率 r 是常数。

N(t)=N0ertN(t)=N_0e^{rt}

这个模型的问题在于，它没有考虑到人口增长的阻滞因素，即人口增长的上限。因此，当预测时间 t 较短时，模型的预测结果是合理的，但是当预测时间 t 较长时，模型的预测结果就会出现较大的误差。

Logistic模型（阻滞增长模型）

由于人口不可能无限制的增长，当人口达到一定数量后，那么增长率就会下降。因此， Verhulst 于 1838 年提出了阻滞增长模型，基本假设是人口增长率 r 随着人口数量 N 的增加而减小。

人口增长率 r 是一个关于人口数量 N 的线性函数：

dNdt=r(N)N\frac{dN}{dt}=r(N)N

人口增长率 r 的函数形式为：

r(N)=r−sN,(r,s>0)r(N)=r-sN,(r,s>0)

其中， r 是人口增长率的最大值， s 是人口增长率的下降速率。

而人口增长率的下降速率与人口容量 K 有关，即：

s=rKs=\frac{r}{K}

结合上述公式，得到人口的增长速度：

dNdt=rN(1−NK)\frac{dN}{dt}=rN(1-\frac{N}{K})

最后得出阻滞增长模型：

N(t)=K1+(KN0−1)e−rtN(t)=\frac{K}{1+(\frac{K}{N_0}-1)e^{-rt}}

MATLAB 实现阻滞增长模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
function [t,N]=logistic(r,K,N0,t0,tf,dt)
    %LOGISTIC 阻滞增长模型
    %   [t,N]=logistic(r,K,N0,t0,tf,dt)
    %   r:人口增长率的最大值
    %   K:人口容量
    %   N0:初始人口数量
    %   t0:初始时间
    %   tf:终止时间
    %   dt:时间步长
    %   t:时间序列
    %   N:人口数量序列
    t=t0:dt:tf;
    N=zeros(1,length(t));
    N(1)=N0;
    for i=2:length(t)
        N(i)=N(i-1)+r*N(i-1)*(1-N(i-1)/K)*dt;
    end
end

例

根据美国人口数据进行预测，美国人口数据如下：

年份17901800181018201830184018501860187018801890
人口/百万3.95.37.29.612.917.123.231.438.650.262.9
年份19001910192019301940195019601970198019902000
人口/百万76.292.2106.5123.2132.2151.3179.3203.2226.5249.6281.4

数据来源：美国人口数据

MATLAB 中提供了 nlinfit 函数，可以用来拟合非线性模型，其语法格式如下：

1
beta = nlinfit(X,Y,modelfun,beta0)

其中， X 是自变量， Y 是因变量， modelfun 是模型函数， beta0 是模型参数的初始值， beta 是模型参数的估计值。

使用 nlinfit 函数进行拟合：

1
2
3
4
5
6
7
8
9
10
% 原始数据
t=1790:10:2000;
N=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.2,92.2,106.5,123.2,132.2,151.3,179.3,203.2,226.5,249.6,281.4];
% 拟合
f=@(beta,t)beta(1)./(1+(beta(1)./N(1)-1).*exp(-beta(2).*t));
beta0=[300,0.5];
beta=nlinfit(t-t(1),N,f,beta0)
% 参数
K=beta(1)
r=beta(2)

得到参数估计值：

1
2
3
4
5
6
K =
  340.0598


r =
    0.0274

预测并绘图：

1
2
3
4
% 预测
[t_hat,N_hat]=logistic(r,K,N(1),t(1),t(end),1);
% 绘图
plot(t,N,'o',t_hat,N_hat,'-','LineWidth',2),grid on,legend('原始数据','预测结果'),xlabel('年份'),ylabel('人口/百万')

图1

可以看出，预测结果与原始数据非常接近。

预测 2080 年的美国人口数量变化趋势：

1
2
3
4
% 预测
[t_hat,N_hat]=logistic(r,K,N(1),t(1),2080,1);
% 绘图
plot(t,N,'o',t_hat,N_hat,'-','LineWidth',2),grid on,legend('原始数据','预测结果'),xlabel('年份'),ylabel('人口/百万')

图2