代码收藏家技术教程 2023-05-14

【续】数学模型——人口增长模型

很抱歉写的不够详细，让一些初学者无法得知数据的来源。本篇文章基于之前的数学模型——人口增长模型（基于python）进行注释，主要解释上一篇文章中，个别需要计算的参数是如何得到的。

需要导入的库如下：

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

上篇及本篇的数据如下所示。其中t表示年份。

data = pd.DataFrame({
    "x": [3.9, 5.3, 7.2, 9.6, 12.9, 17.1, 23.2, 31.4, 38.6, 50.2, 62.9, 76, 92, 105.7,                         
          122.8, 131.7, 150.7, 179.3, 203.2, 226.5, 248.7, 281.4],
    "t": range(22),
    "r": [0.2949, 0.3113, 0.2986, 0.2969, 0.2907, 0.3012, 0.3082, 0.2452, 0.2435, 0.242, 
           0.2051, 0.1914, 0.1614, 0.1457, 0.1059, 0.1059, 0.1579, 0.1464, 0.1161, 0.1004, 
           0.1104, 0.1349]
})

一、对于初始的增长模型

这里假定了，人口增长率是一个定值（当然现实中肯定不是，在之后会一步步放开这个假设）。根据普通最小二乘回归，可得参数 $r$ 。算法函数如下。（这里不放算法的具体推到过程，后续会专门做一篇）

def ols(x, y):
    x_mean = np.mean(x)  # x的均值
    y_mean = np.mean(y)
    x_mean_square = np.array([(i - x_mean) ** 2 for i in x])  # 计算xi-x_bar的平方
    y_mean_square = np.array([(i - y_mean) ** 2 for i in y])  # 计算yi-y_bar的平方
    xy_mean = np.array([((i - x_mean) * (j - y_mean)) for i, j in zip(x, y)])

    b = xy_mean.sum() / x_mean_square.sum()
    a = y_mean - b * x_mean

    fitvalue = np.array([])
    fitvalue = np.append(fitvalue, [a + b * i for i in x])  # 计算拟合值

    y_pred_square = np.array([(i - j) ** 2 for i, j in zip(y, fitvalue)])  # 残差平方和SSR
    SE = np.sqrt(y_pred_square.sum() / (len(y) - 2))
    R_square = 1 - np.array(y_pred_square.sum() / np.sum([(i - y_mean) ** 2 for i in y]))

    return a, b, fitvalue, SE, R_square

这个函数输入一个自变量x和一个因变量y，输出结果中，a表示截距，b表示回归系数，fitvalue表示模型拟合的因变量的值，SE表示标准误，当需要使用此模型进行预测或者解释回归系数的显著性时（显著性：该自变量是否真的能影响因变量），就得好好观察这个数值。R_square表示拟合模型的可决系数，这个系数越接近于1，表明模型拟合的越好，也意味着，x和y之间存在着明显的线性关系。

根据以下代码，即可得出 $r_0=0.2020$

x = data.x
y = np.log(x)
x0 = 6.0496
t = data.t

a, b, fitvalue, SE, r2 = ols(t, y)
print("截距：{}\n回归系数：{}\n标准误：{}\nR方值：{}".format(a, b, SE, r2))

二、对增长率进行改进的增长模型

在此处，已经放开了前面“增长率是不变”的假设，上述数据集data中的变量r即为人口增长率的变化。通过下面的代码，可知 $r_0=0.3252,r_1=0.0114$ ，此时会发现模型的输出结果，会发现 $r_1$ 是负的，但并不是算错了，在上篇文章中，假定的变化率方程是 $r(t)=r_0-r_1*t$ ，已经带有一个负号，所以结果是没有错的。

r = data.r
r0, r1, f1, se, R2 = ols(y=r, x=t)
print(r0, r1)

三、 Logistics模型

很遗憾，由于初始值等原因，我无法很好的复刻出书中给的数值，但这边依旧给出非线性最小二乘的估计过程，望读者见谅。

代码如下：

from scipy.optimize import curve_fit

def func(x, R, xm):
    return R * (1 - x / xm)

x = data['x']
r = data['r']

params = curve_fit(func, x, r)
r0 = params[0][0]
xm = params[0][1]
print(r0, xm)

def logistics(x, t):
    return np.array(r0 * x * (1-x/xm) + 0 * t)

T = np.arange(0, 30, 1)
x = odeint(logistics, 3.9, T)
plt.scatter(data['t'], data['x'], c='r')
plt.plot(x)
plt.show()

根据上述代码计算得到的两个参数分别为 $r=0.2805$ ， $x_m=352.0477$ ，代入模型的结果如下图所示。