読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

微分方程式を含むモデルのパラメータ推定(フェルフルストの人口モデル)

R 微分方程式 数値計算 最適化

要旨

  • R で微分方程式を含むモデルのパラメータ推定を行う.
  • まずは単純なモデル(解析的に解ける常微分方程式)で練習.
  • deSolve パッケージの ode 関数(4次ルンゲ=クッタ法)で微分方程式を数値積分しつつ, minpack.lm パッケージの nls.lm 関数(Levenberg–Marquardt アルゴリズム)でパラメータ推定.

フェルフルストの人口モデル

人口増加に関するモデルはマルサスが1798年に『人口論』で提案したものが有名だ. マルサスのモデルだと人口は未来に向かって際限なく増加する.

フェルフルスト(Verhulst;『微分方程式で数学モデルを作ろう』ではヴェアフルストと表記されている)は1837年にマルサス・モデルの修正を提案した.

フェルフルストは人口の変化量  dN/dt は,

  • 現在の人口レベル N
  • 未利用の人口資源の上限に対する比 (1- N/K)K は人口の上限)

のそれぞれに比例すると仮定した. つまり,

 \displaystyle  \frac{dN}{dt} = r N \left(1 - \frac {N}{K} \right).

この方程式をロジスティック方程式という.

N について解くと

 \displaystyle  N = \frac{K}{ 1 + C e^{-rt} },

ここで  C = K/N(0) -1 .

パラメータの推定

下図は年度ごとの米国の人口を表している. 『微分方程式で数学モデルを作ろう』(p.7)から持ってきた.

f:id:abrahamcow:20150221225230p:plain

t = seq(1820, 1930, by = 10)
length(t)
population = c(9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2)

par(mgp=c(2,0.8,0))
ylab <- bquote(.("population (") ~ 10^6 ~")")
plot(t, population, pch = 21, bg = "gray", xlab = "year", ylab = ylab)

このデータにロジスティック方程式の解を当てはめる.

nls 関数でやろうとしてもなかなかうまくいかない.

> nls(population ~ K/(1+(K/9.6-1)*exp(r*(1:12))),start = c(K = 3.9, r = 0.3))
 以下にエラー nls(population ~ K/(1 + (K/9.6 - 1) * exp(r * (1:12))), start = c(K = 3.9,  : 
   勾配が特異です 

そこで deSolve パッケージの ode 関数で微分方程式を数値積分しつつ, minpack.lm パッケージの nls.lm 関数でパラメータを推定する.

nls.lm 関数は Levenberg–Marquardt アルゴリズムを使っているらしい.

library(deSolve)
logistic = function(t, state, parameters) {
  with(as.list(c(state, parameters)), {    
    dN = r * N*(1 - N/K)
    list(c(dN))
  })
}

fn1 = function(params) {
  times = 1:12
  state = c(N = population[1])
  
  ## solve the ODE.
  out = ode(y = state, times = times, func = logistic, parms = params, method = "rk4")
  
  ## modeled - observed
  fn1 = out[, "N"] - population
}

library(minpack.lm)
fit1 = nls.lm(par =  c(r = 0.1, K=100), fn = fn1)

こんな感じでパラメータが求まる.

> summary(fit1)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
r 3.146e-01  1.978e-03  159.03  < 2e-16 ***
K 1.981e+02  3.729e+00   53.14 1.35e-13 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.627 on 10 degrees of freedom
Number of iterations to termination: 11 
Reason for termination: Relative error in the sum of squares is at most `ftol'. 

プロットしてみる.

f:id:abrahamcow:20150221231440p:plain

params = coef(fit2)
modeled = ode(y = state, times = times, func = logistic, parms = params, 
                    method = "rk4")
plot(t, population, pch = 21, bg = "gray", xlab = "year", ylab = ylab)
lines(t,modeled[,"N"],col="blue3", lwd=2)
legend("bottomright", legend = c("Observed data", "predicted value"), 
       lty = c(NA, 1, 1), pch = c(21, NA), col = c("black",  "blue3"), pt.bg = "gray", 
       bty = "n", lwd=2)

微分方程式で数学モデルを作ろう』によると, フェルフルストはパラメータを, r = 0.3134,  K= 197 \times 10^6 と選んだらしいが, 今回の推定値のほうが少しだけ当てはまりがいい.

f:id:abrahamcow:20150221233421p:plain

state = c(N = population[1])
times = 1:12
parameters = c(r = 0.3134, K=197)
modeled2 = ode(y = state, times = times, func = logistic, parms = parameters, 
               method = "rk4")
plot(t, population, pch = 21, bg = "gray", xlab = "year", ylab = ylab)
lines(t,modeled[,"N"],col="blue3", lwd=2)
lines(t,modeled2[,"N"],col="red3", lwd=2)
legend("bottomright", legend = c("Observed data", "predicted value", "predicted value (Verhulst)"), 
       lty = c(NA, 1, 1), pch = c(21, NA,NA), col = c("black","blue3","red3"), pt.bg = "gray", 
       bty = "n", lwd=2)

残差平方和はそれぞれこんな感じ.

> ###residual
> sum((population - modeled[,"N"])^2) #nls.lm
[1] 3.931386
> sum((population - modeled2[,"N"])^2) #Verhulst
[1] 7.203779