人口問題 - 驚異のアニヲタ社会復帰の予備

MikuHatsune2011-12-20

微分方程式で数学モデルを作ろう　第一章
　
 Thomas Robert Malthusの人口論を例にとって考えてみる。
$t$ ：ある時刻
$N(t)$ ：ある時刻 $t$ での人口数
$\alpha$ ：出生数を規定する比例定数
$\beta$ ：死亡数を規定する比例定数
$\gamma$ ： $\alpha - \beta$
$\delta$ ：微小変化量
　
ある短い時刻 $\delta t$ の間の全人口の増加は
$\begin{align}\delta N&=\alpha N\delta t\hspace{3}-\hspace{3}\beta N\delta t\\&=\gamma N\delta t\\\frac{\delta N}{\delta t}&=\gamma N\end{align}$

これはオレでも解くことができて、初期状態を $N_0$ とでもおけば
$N(t)=N_0 \hspace{3}e^{\gamma t}$
$\gamma$ の影響が大きい。

time <- seq(from=0, to=10, length=20)
gamma <- c(0.1, 0, -0.1)
N_0 <- 1000
# 地味に計算する場合
N <- matrix(0, nr=length(gamma), nc=length(time))
N[, 1] <- N_0
for(j in 1:length(gamma)){
for(i in 2:length(time)){
	N[j, i] <- N[j, i-1] + gamma[j] * N[j, i-1]
}
}
matplot(t(N), type="l", lty=1, ylim=c(0,10000))
legend(1, 10000, paste("gamma", c(">", "=", "<"), 0), col=1:3, lwd=3)

# 関数が解ける場合
curvelist <- lapply(mapply(rep, 0, 1:length(gamma)), unique)
for(i in 1:length(gamma)){
	curvelist[[i]] <- function(x){
		return(N_0 * exp(gamma[i] * x))
	}
}
xlim <- range(time)
ylim <- c(0, 10000)
for(i in 1:length(curvelist)){
	plot(curvelist[[i]], col=i, xlim=xlim, ylim=ylim)
	par(new=TRUE)
}
legend(1, 10000, paste("gamma", c(">", "=", "<"), 0), col=1:3, lwd=3)

これを実データに当てはめてみると
$\begin{matrix}year&population \hspace{3}in\hspace{3} US\hspace{3} (10^6)\\1790&3.9\\1800&5.3\\1810&7.2\end{matrix}$
1790年を時刻0とし、10年刻みで $t$ が増加するとすると、
$\begin{align}N(1)&=5.3*10^6\\&=3.9*10^6\hspace{3}e^{\gamma}\\\gamma &=log\frac{5.3}{3.9}\\&=0.307\end{align}$
これらより、 $N(2)=7.3*10^6$ となる。
　
用量反応曲線のように推定することも可能である。

EXD.3

で、c、d、eの3つのパラメータを推定して以下の曲線を描ける。
$f(x)=c+(d-c)e^{-\frac{x}{e}}$

library(drc)
x <- 0:2
y <- c(3.9, 5.3, 7.2)*10^6
plot(drm(y ~ x, fct=EXD.3()), axes=FALSE)
axis(1, 0:2, 1790 + (0:2)*10)
axis(2)