人口增長模型

来源:https://www.cnblogs.com/hznudmh/archive/2022/06/15/16343159.html
-Advertisement-
Play Games

1. 指數增長模型 設第今年的人口為 \(x_0\)，年增長率為 \(r\)，預測 \(k\) 年後的人口為 \(x_k\)，則有 \[ x_k = x_0(1+r)^k \tag{1} \] 這個模型的前提是年增長率 \(r\) 在 \(k\) 年內保持不變. 利用 (1) 式可以根據人口估計年增 ...


目錄

1. 指數增長模型

設第今年的人口為 \(x_0\)，年增長率為 \(r\)，預測 \(k\) 年後的人口為 \(x_k\)，則有

\[x_k = x_0(1+r)^k \tag{1} \]

這個模型的前提是年增長率 \(r\) 在 \(k\) 年內保持不變.

利用 (1) 式可以根據人口估計年增長率，例如經過 \(n\) 年後人口翻了一番，那麼這期間的年平均增長率約為 \((70/n)\%\).
證明：

\[\begin{aligned} &x_n = x_0(1+r)^n \Longrightarrow 2x_0 = x_0(1+r)^n \Longrightarrow (1+r)^n = 2 \Longrightarrow \ln(1+r)=\frac{\ln 2}{n}\\ \Longrightarrow & 1+r=e^{\frac{\ln 2}{n}} \Longrightarrow r = e^{\frac{\ln2}{n}}-1 \sim \frac{\ln2}{n} \approx (\frac{70}{n})\% \end{aligned} \]

1.1 人口增長模型的建立

將人口看作連續時間 \(t\) 的連續可微函數 \(x(t)\). 於是有

\[\frac{\mathrm{d}x}{\mathrm{d}t} = rx,\quad x(0) = x_0 \tag{2} \]

求解上述微分方程：

\[\begin{aligned} &原式 \ = \frac{1}{x}\mathrm{d}x = r\mathrm{d}t\\ \\ \Longrightarrow & \ln x = rt+C\\ \\ \Longrightarrow & x = e^{rt}\cdot C\\ \end{aligned} \]

由於 \(x(0) = x_0\)，代入上式可得 \(x_0 = C\)，可以解出

\[x(t) = x_0e^{rt} \tag{3} \]

可以用 R 語言定義指數增長模型方程

## 指數增長模型方程
fun1 <- function(t) x0 * exp(r*t)

1.2 參數估計

以書本 P142 美國人口數據為例

1.1.1 線性最小二乘估計

對 (3) 取對數，有 \(\ln x = \ln x_0 + rt\)，於是 \(\ln x\) 與 \(t\) 存線上性關係，令10年為一個單位.
下麵用 R 進行參數估計.

> ## 美國人口數據的導入
> x <- c(3.9, 5.3, 7.2, 9.6, 12.9, 17.1, 23.2, 31.4,
+        38.6, 50.2, 62.9, 76.0, 92.0, 105.7, 122.8,
+        131.7, 150.7, 179.3, 203.2, 226.5, 248.7, 281.4)
> ## 線性最小二乘估計
> y <- log(x)
> t <- sequence(length(x), 0, 1)
> lm(y~t) -> sol
> summary(sol)

Call:
lm(formula = y ~ t)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.43901 -0.19764  0.00609  0.23405  0.32145 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.79999    0.10538   17.08 2.14e-13 ***
t            0.20201    0.00859   23.52 4.80e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2556 on 20 degrees of freedom
Multiple R-squared:  0.9651,	Adjusted R-squared:  0.9634 
F-statistic: 553.1 on 1 and 20 DF,  p-value: 4.8e-16
> exp(1.79999)
[1] 6.049587

結果分析：

\[\ln x = 0.2020t + 1.79999 \]

於是\(r = 0.2020/10年, \ x_0 = e^{1.79999}=6.0496\)

通過圖像看看擬合效果：

# 擬合圖像
x0 <- exp(1.79999)
r <- 0.20201
y_hat1 <- fun1(t)
plot(t, x)
lines(t, y_hat1, col="red")

1.1.2 基於數值微分的參數估計

對人口數據作數值微分，進而計算增長率，將增長率的平均值作為 \(r\) 的估計；\(x_0\) 直接採用原始數據.
理論依據：

\[\frac{\mathrm{d }x}{\mathrm{d}t} = rx \Longrightarrow \frac{\mathrm{d}x/\mathrm{d}t}{x} = r \]

數值微分可用中點公式算近似值：

\[x'(t_k) = \frac{x_{k+1}-x_{k-1}}{2\Delta t}, \ k = 1, 2, \cdots, n-1\\ x'(t_0) = \frac{-3x_0+4x_1-x_2}{2\Delta t}, \ x'(t_n) = \frac{x_{n-2}-4x_{n-1}+3x_n}{2\Delta t} \]

下麵用 R 語言計得出 \(r\) 的估計值

> ## 數值微分估計
> dx <- c()
> dx[1] <- (-3*x[1]+4*x[2]-x[3])/2 # 數值微分中點公式的應用
> for (i in 2:(length(x)-1)){
+   dx[i] <- (x[i+1]-x[i-1])/2 # 數值微分中點公式的應用
+ }
> dx[length(x)] <- (x[length(x)-2]-4*x[length(x)-1]+3*x[length(x)])/2 # 數值微分中點公式的應用
> R <- dx/x # 增長率序列
> r <- mean(R) # 其平均值就是r的估計
> r
[1] 0.02051861

通過圖像看看擬合效果

# 擬合圖像
x0 <- x[1]
y_hat2 <- fun1(t)
plot(t, x)
lines(t, y_hat2, col="red")

1.3 改進的指數增長模型

在前面的模型中，我們假定了每年人口的增長率是固定的，顯然這是存在缺陷的，每年的增長率應是時間的函數，即增長率為 \(r(t)\).

我們首先畫出 \(r(t)\) 與 \(t\) 的散點圖（這裡我們的增長率用的是1.2.2中做數值微分得到的增長率序列），看看他們存在什麼函數關係

## 畫出 r~t 散點圖
plot(t,R)

可以看出兩者疑似存線上性關係，不妨做線性擬合

> ## 線性擬合
> lm(R~t) -> sol2
> summary(sol2)

Call:
lm(formula = R ~ t)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.059303 -0.007904 -0.001237  0.015108  0.051549 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.3252461  0.0117048   27.79  < 2e-16 ***
t           -0.0114343  0.0009541  -11.98 1.39e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02839 on 20 degrees of freedom
Multiple R-squared:  0.8778,	Adjusted R-squared:  0.8717 
F-statistic: 143.6 on 1 and 20 DF,  p-value: 1.391e-10

結果分析：

\[r = 0.3252461 - 0.0114343t \]

且回歸方程的檢驗統計量\(F=143.6, \ p = 1.391e-10 << 0.01, \ R^2 = 0.8778\)，認為回歸方程的擬合程度不錯，兩者存在非常顯著的線性關係；常數項檢驗 \(t = 27.79, \ p << 0.01\)，一次項繫數檢驗 \(t = -11.98, \ p = 1.39e-10 << 0.01\)，認為回歸方程繫數顯著不為0.

看看擬合圖像

# 擬合圖像
plot(t,R)
abline(sol2, col = "red")

記 \(r = r_0 - r_1t\)，於是 \((2)\) 式變為

\[\frac{\mathrm{d}x}{\mathrm{d}t} = (r_0-r_1t)x, \quad x(0) = x_0 \tag{4} \]

用分離變數的方法求解上述微分方程，\((4)\) 可化為

\[\begin{aligned} &\frac{1}{x}\mathrm{d}x = (r_0-r_1t)\mathrm{d}t\\ \Longrightarrow &\int \frac{1}{x}\mathrm{d}x = \int (r_0-r_1t)\mathrm{d}t\\ \Longrightarrow &\ln x = r_0-\frac{1}{2}r_1t^2 + C\\ \Longrightarrow &x = C'\exp(r_0t-\frac{1}{2}r_1t^2)\\ \end{aligned} \]

將 \(x(0) = x_0\) 代入，即可得 \(C' = x_0\)，於是解得

\[x(t) = x_0\exp(r_0t-\frac{1}{2}r_1t^2) \]

可以用 \(R\) 語言定義上述函數

## 定義函數
fun2 <- function(t) x0 * exp(r0*t - r1*t^2/2)

下麵用 \(R\) 語言做指數增長模型的擬合

## 利用模型估計人口並做擬合圖像
r0 <- 0.3252461
r1 <- 0.0114343
x0 <- x[1]
t <- sequence(length(x), 0, 1)
y_hat3 <- fun2(t)
plot(t, x)
lines(t, y_hat3, col = "red")

可以看出改進的指數增長模型的擬合程度明顯好於改進前的指數增長模型.

2. logistic 模型

考慮自然資源、環境條件等因素對人口增長起阻礙作用，並且隨著人口的增加，阻礙作用越來越大.

2.1 logistic 模型的建立

隨著人口的增長，資源和環境的阻礙作用體現在對增長率 \(r\) 的影響上，於是 \(r(x)\) 是 \(x\) 的函數，主觀上我們認為 \(r(x)\) 是減函數，我們可以畫出兩者之間的散點圖來大致判斷兩者之間滿足何種函數關係（這裡我們的增長率用的是1.2.2中做數值微分得到的增長率序列）

## r~x 散點圖
plot(x, R)

散點圖顯示兩者之間存線上性關係，不妨設 \(r(x) = a + bx\)，為了賦予繫數 \(a,b\) 實際含義，引入兩個參數：

  • 內稟增長率 \(r_0\) \(\quad r_0\) 是 \(x=0\) 的增長率，於是 \(a = r_0\);
  • 人口容量 \(x_m\) \(\quad x_m\) 是資源和環境所能容納的最大人口數量，也就是當 \(x=x_m\) 後人口不再增長，即 \(r(x_m) = r + bx_m = 0\)，解得 \(b = -\frac{r}{x_m}\).

於是

\[r(x) = r_0(1-\frac{x}{x_m}) \]

那麼 \((2)\) 變為

\[\frac{\mathrm{d}x}{\mathrm{d}t} = r_0x(1-\frac{x}{x_m}), \quad x(0) = x_0 \tag{5} \]

上式稱為 logistic 方程.

將上式化為

\[\begin{aligned} &\frac{x_m}{x(x_m-x)} \mathrm{d}x = r\mathrm{d}t\\ \Longrightarrow &\int(\frac{1}{x} + \frac{1}{x_m-x})\mathrm{d}x = \int r_0\mathrm{d}t\\ \Longrightarrow &\ln{x} - \ln{(x_m-x)} = r_0t + C\\ \Longrightarrow &\frac{x}{x_m-x}=C'e^{r_0t}\\ \Longrightarrow &x = \frac{x_m}{1+C''e^{-r_0t}}\\ \end{aligned} \]

將點 \((t,x) = (0,x_0)\) 代入上式可得

\[C'' = \frac{x_m}{x_0}-1 \]

綜上就可求得

\[x(t) = \frac{x_m}{1+(\frac{x_m}{x_0}-1)e^{-r_0t}} \tag{6} \]

上式稱為 logistic 曲線. \((5),(6)\) 統稱為 logistic 模型.

利用 R 語言定義上述函數

## 定義函數
fun3 <- function(t) xm/(1+(xm/x0-1)*exp(-r0*t))

2.2 參數估計

2.2.1 線性最小二乘估計

\((5)\) 式可改寫為

\[\frac{\mathrm{d}x/\mathrm{d}t}{x} = r_0 - \frac{r_0}{x_m}x, \quad x(0) = x_0 \]

下麵用 R 語言做線性估計

> ## 線性最小二乘估計
> lm(R~x) -> sol3
> summary(sol3)

Call:
lm(formula = R ~ x)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.076819 -0.021568  0.000946  0.023945  0.078560 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.805e-01  1.243e-02  22.573 1.06e-15 ***
x           -7.968e-04  9.763e-05  -8.162 8.55e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03902 on 20 degrees of freedom
Multiple R-squared:  0.7691,	Adjusted R-squared:  0.7575 
F-statistic: 66.61 on 1 and 20 DF,  p-value: 8.546e-08

結果分析：

\[\frac{\mathrm{d}x/\mathrm{d}t}{x} = 0.2805 - 7.968e-04 x \]

且回歸方程的檢驗統計量\(F=66.61, \ p = 8.546e-08 << 0.01, \ R^2 = 0.7691\)，認為回歸方程的擬合程度不錯，兩者存在非常顯著的線性關係；常數項檢驗 \(t = 22.573, \ p << 0.01\)，一次項繫數檢驗 \(t = -8.162, \ p << 0.01\)，認為回歸方程繫數顯著不為0.

於是 \(r_0 = 0.2805, \ x_m = 0.1805/7.968e-04 = 352.0331\)，畫出擬合圖像

# 擬合
r0 = 0.2805
xm = 0.2805/7.968e-04
y_hat4 <- fun3(t)
plot(t, x)
lines(t, y_hat4, col = "red")

2.2.2 非線性最小二乘法

考慮 Levenberg-Marquadt 方法求解，\((6)\) 式可化為

\[x(t) - \frac{x_m}{1+(\frac{x_m}{x_0}-1)e^{-r_0t}} = 0 \]

帶估計參數為 \((x_m,x_0,r_0)\)，給定初始值 \((400,4,0.25)\)，用 R 求解

> ## 非線性最小二乘估計
> library(pracma)
> fun4 <- function(z) {
+   f <- numeric(22)
+   for (i in 1:22){
+     f[i] <- x[i] - z[1]/(1+(z[1]/z[2]-1)*exp(-z[3]*(i-1)))
+   }
+   f
+ }
> sol4 <- lsqnonlin(fun4, x0 = c(400,4,0.25))
> sol4
$x
[1] 443.9948201   7.6962896   0.2154783

$ssq
[1] 458.2113

$ng
[1] 0.0002470027

$nh
[1] 3.847784e-06

$mu
[1] 7.431693e-06

$neval
[1] 22

$errno
[1] 1

$errmess
[1] "Stopped by small x-step."

於是可得參數的估計值為 \(x_m = 443.9948201, \ x_0 = 7.6962896, \ r_0 = 0.2154783\).

下麵將估計的參數代入 \((6)\)，得出擬合圖像

# 擬合
xm = 443.9948201
x0 = 7.6962896
r0 = 0.2154783
y_hat5 <- fun3(t)
plot(t, x)
lines(t, y_hat5, col = "red")

2.3 比較

為了比較線性最小二乘估計和非線性最小二乘估計的擬合程度，引入均方誤差，均方誤差越小說明擬合程度越高

\[MSE = \frac{1}{n}\sum_{i=1}^n(x_i-\hat{x}_i)^2 \]

下麵用 R 計算兩種方法的均方誤差

> ## 比較
> MSE1 <- sum((y_hat4-x)^2)/length(x)
> MSE2 <- sum((y_hat5-x)^2)/length(x)
> print(MSE1)
[1] 127.1399
> print(MSE2)
[1] 20.82779

發現 \(MS1 = 127.1399 > 20.82779 = MSE2\)，因此非線性最小二乘參數估計的擬合效果更好.

不妨再計算一下改進的指數增長模型的均方誤差

> MSE3 <- sum((y_hat3-x)^2)/length(x)
> MSE3
[1] 51.51463

可以發現，logistic模型的非線性最小二乘估計參數的結果要好於改進的指數增長模型.


您的分享是我們最大的動力!

-Advertisement-
Play Games
更多相關文章
  • tree-shaking
    來源 tree-shaking 最早由 Rich Harris 在 rollup 中提出。 為了減少最終構建體積而誕生。 以下是 MDN 中的說明： tree-shaking 是一個通常用於描述移除 JavaScript 上下文中的未引用代碼(dead-code) 行為的術語。 它依賴於 ES201 ...
  • 【ES6】迭代器
    簡介 在 JavaScript 中，迭代器是一個對象，它定義一個序列，併在終止時可能返回一個返回值。 更具體地說，迭代器是通過使用 next() 方法實現迭代器協議的任何一個對象，該方法返回具有兩個屬性的對象： value，這是序列中的 next 值；和 done ，如果已經迭代到序列中的最後一個值 ...
  • node前後端交互（Express）
    1. Express框架是什麼 1.1 Express是一個基於Node平臺的web應用開發框架，它提供了一系列的強大特性，幫助你創建各種Web應用。我們可以使用 npm install express 命令進行下載。 1.2 Express初體驗 // 引入express框架 const expr ...
  • Java遞歸查找層級文件夾下特定內容的文件
    遞歸查找文件 引言 或許是文件太多,想找某個文件又忘記放哪了;又或者是項目改造,需要將外部調用介面進行改造,項目太多,又無法排查。那麼怎麼快速找到自己想要的內容就是一件值得思考的事情了。 根據特定內容尋找文件位置 package com.lizi.globalexception.Utils; imp ...
  • MyBatis流式查詢
    大家好，本文我將繼續來剖析SpringCloud中負載均衡組件Ribbon的源碼。本來我是打算接著OpenFeign動態代理生成文章直接講Feign是如何整合Ribbon的，但是文章寫了一半發現，如果不把Ribbon好好講清楚，那麼有些Ribbon的細節理解起來就很困難，所以我還是打算單獨寫一篇文章 ...
  • phpcurl函數類模擬Curl get post header refer攜帶Cookie模擬訪問來源Refer模擬UseaAgent
    phpcurl函數類模擬Curl get post header refer攜帶Cookie模擬訪問來源Refer模擬UseaAgent ...
  • SpringBoot集成feign的方法分享
    一個工作5年的粉絲找到我，他說參加美團面試，遇到一個基礎題沒回答上來。 這個問題是：“資料庫連接池有什麼用？以及它有哪些關鍵參數”？ 我說，這個問題都不知道，那你項目裡面的連接池配置怎麼設置的？ 你們猜他怎麼回答。懂得懂得啊。 好的，關於這個問題，我們來看看普通人和高手的回答。 普通人： 資料庫連接 ...
  • 流量管制-令牌桶與漏桶
    Principle of token bucket 隨著互聯網的發展，在處理流量的方法也不僅僅為 first-come，first-served，而在共用網路中實現流量管理的基本機制就是排隊。而公平演算法則是實現在優先順序隊列中基於哪些策略來排隊的”公平隊列“。Token Bucket 則是為公平排隊提 ...
一周排行
    -Advertisement-
    Play Games
  • 移動開發（一）：使用.NET MAUI開發第一個安卓APP
    移動開發（一）：使用.NET MAUI開發第一個安卓APP 對於工作多年的C#程式員來說，近來想嘗試開發一款安卓APP，考慮了很久最終選擇使用.NET MAUI這個微軟官方的框架來嘗試體驗開發安卓APP，畢竟是使用Visual Studio開發工具，使用起來也比較的順手，結合微軟官方的教程進行了安卓 ...
  • wpf ToggleButton選中效果和一個登錄界面
    前言 QuestPDF 是一個開源 .NET 庫，用於生成 PDF 文檔。使用了C# Fluent API方式可簡化開發、減少錯誤並提高工作效率。利用它可以輕鬆生成 PDF 報告、發票、導出文件等。 項目介紹 QuestPDF 是一個革命性的開源 .NET 庫，它徹底改變了我們生成 PDF 文檔的方 ...
  • 彈幕樹洞項目功能新增篇
    項目地址 項目後端地址： https://github.com/ZyPLJ/ZYTteeHole 項目前端頁面地址： ZyPLJ/TreeHoleVue (github.com) https://github.com/ZyPLJ/TreeHoleVue 目前項目測試訪問地址： http://tree ...
  • 第27篇 sqlserver2022詳細安裝步驟
    話不多說，直接開乾 一.下載 1.官方鏈接下載： https://www.microsoft.com/zh-cn/sql-server/sql-server-downloads 2.在下載目錄中找到下麵這個小的安裝包 SQL2022-SSEI-Dev.exe，運行開始下載SQL server； 二. ...
  • .NET 開源高性能 MQTT 類庫
    前言 隨著物聯網（IoT）技術的迅猛發展，MQTT（消息隊列遙測傳輸）協議憑藉其輕量級和高效性，已成為眾多物聯網應用的首選通信標準。 MQTTnet 作為一個高性能的 .NET 開源庫，為 .NET 平臺上的 MQTT 客戶端與伺服器開發提供了強大的支持。 本文將全面介紹 MQTTnet 的核心功能 ...
  • Serilog文檔翻譯系列（六） - 可用的接收器、增強器、格式化輸出
    Serilog支持多種接收器用於日誌存儲，增強器用於添加屬性，LogContext管理動態屬性，支持多種輸出格式包括純文本、JSON及ExpressionTemplate。還提供了自定義格式化選項，適用於不同需求。 ...
  • 警惕 Visual Studio 屬性求值副作用導致邏輯不符合預期
    目錄簡介獲取 HTML 文檔解析 HTML 文檔測試參考文章 簡介 動態內容網站使用 JavaScript 腳本動態檢索和渲染數據，爬取信息時需要模擬瀏覽器行為，否則獲取到的源碼基本是空的。 本文使用的爬取步驟如下： 使用 Selenium 獲取渲染後的 HTML 文檔 使用 HtmlAgility ...
  • [使用目前最新版]HybridCLR6.9.0+YooAsset2.2.4實現純C# Unity熱更新方案 (一)
    1.前言 什麼是熱更新 游戲或者軟體更新時，無需重新下載客戶端進行安裝，而是在應用程式啟動的情況下，在內部進行資源或者代碼更新 Unity目前常用熱更新解決方案 HybridCLR，Xlua，ILRuntime等 Unity目前常用資源管理解決方案 AssetBundles，Addressable， ...
  • 在 ASP.NET Core Web API 中使用操作篩選器統一處理通用操作
    本文章主要是在C# ASP.NET Core Web API框架實現向手機發送驗證碼簡訊功能。這裡我選擇是一個互億無線簡訊驗證碼平臺，其實像阿裡雲，騰訊雲上面也可以。 首先我們先去 互億無線 https://www.ihuyi.com/api/sms.html 去註冊一個賬號 註冊完成賬號後，它會送 ...
  • 第28篇 如何.net中實現高效可靠數據同步api
    通過以下方式可以高效，並保證數據同步的可靠性 1.API設計 使用RESTful設計，確保API端點明確，並使用適當的HTTP方法（如POST用於創建，PUT用於更新）。 設計清晰的請求和響應模型，以確保客戶端能夠理解預期格式。 2.數據驗證 在伺服器端進行嚴格的數據驗證，確保接收到的數據符合預期格 ...
所有分類