人口ピラミッドの形を問う問題があったのでシミュレーションしてみる。
死亡率と出生率は適当なS字曲線に従い、死亡率の低下が先にきてから、出生率の低下が起こる。
出生率は全人口に対して何人生まれたかの‰
死亡率は全人口に対して何人死んだかの‰だが、新生児ほど死にやすい、かつ、高齢ほど死にやすい、というふうにしようと思って、ある年齢で死亡率が最小だとすると、人口ピラミッドの最高年齢の人が死亡率
でこれが最大とすると、新生児の死亡率は
、年齢
の最小の死亡率
として、線形関係にした。
# S字曲線を作成する関数 hill1 <- function (dose, parm){ parmMat <- matrix(c(parm, 1), 1, 5, byrow = TRUE) parmMat[, 2] + (parmMat[, 3] - parmMat[, 2])/((1 + exp(parmMat[, 1] * (dose - parmMat[, 4])))^parmMat[, 5]) } # 出生率早期・後期、死亡率早期・後期、出生率と死亡率のEC50の順 bdrate <- list(c(0.077, 0.00127), c(0.6, 0.05), c(180, 100)) tsrate <- vector("list", length(bdrate)) # 時系列 x <- 1:300 # 年代 for(i in seq(tsrate)){ parms <- c(0.08, bdrate[[i]][2], bdrate[[i]][1], bdrate[[3]][i]) tsrate[[i]] <- hill1(x, parms) } plot(x, tsrate[[1]], type="l", xlab="", ylab="Crude birth rate", cex.axis=1.3, cex.lab=1.6, lwd=3) par(new=TRUE) plot(x, tsrate[[2]], type="n", xlab="Year", ylab="", axes=FALSE, cex.lab=1.6) axis(4, cex.axis=1.5) lines(x, tsrate[[2]], col=4, lwd=3) mtext("Mortality rate", s=4, line=2.5, cex=1.6, col=4)
年齢に応じた死亡率に変換する。
start.max.age <- 50 prop1 <- -(2/start.max.age^2) * seq(start.max.age) + 2/start.max.age # 行が時代、列が年齢の人口行列 pop <- matrix(0, nr=length(x), nc=start.max.age+length(x)) pop[1, seq(prop1)] <- 100000*prop1 # 初期値 min.mortality.age <- 20 # 最も死亡率の低い年齢 min.mortality <- 0.1 # 最高年齢に対する、最も死亡率の低い年齢の死亡率 newborn.mortality <- 0.5 # 最高年齢に対する新生児(0歳)の死亡率 AMR <- function(j){ d <- d.max.age if(j <= min.mortality.age){ r <- (min.mortality-newborn.mortality)*d/min.mortality.age*j + newborn.mortality*d } else { r <- (1-min.mortality)*d/(Max.age-min.mortality.age) * (j-min.mortality.age) + min.mortality*d } r } adj.mortality <- mapply(AMR, seq(Max.age)) plot(adj.mortality, type="o", ylim=c(0, max(adj.mortality)), xlab="Age", ylab="Adjusted mortality")
人口ピラミッドの形は三角形から始まって、いまの日本のように老年人口が増えて年少が少ないみたいな紡錘形になっていく。シミュレーションでは150歳の人とかでてきてしまうけど、年齢による死亡率をもう少しいじったらたぶんもっとキレイになる。
人口爆発は死亡率の低下より出生率が低下し始めたころに起こるらしい。
for(y in 1:(nrow(pop)-1)){ Max.age <- which(pop[y,]==0)[1] - 1 d.max.age <- 2*tsrate[[2]][y]/((newborn.mortality+min.mortality)*min.mortality.age + (1+min.mortality)*(Max.age-min.mortality.age)) adj.mortality <- mapply(AMR, seq(Max.age)) No.death <- pop[y, seq(adj.mortality)] * adj.mortality pop[y+1, seq(adj.mortality)+1] <- pop[y, seq(adj.mortality)]-No.death pop[y+1, 1] <- sum(pop[y, seq(Max.age)])*tsrate[[1]][y] # 新しく生まれた人 } saveGIF({ #par(ask=TRUE, mar=c(5, 4.5, 4, 4), mfrow=c(2,2)) for(y in seq(nrow(pop))){ tmp <- pop[y ,seq(100)] par(mar=c(5, 4.5, 4, 4), mfrow=c(2,2)) barplot(tmp, horiz=TRUE, xlim=c(0, max(apply(pop, 1, max))), xlab="人口軸固定", ylab="年齢軸100歳まで") axis(2) barplot(pop[y,], horiz=TRUE, xlab="人口軸可変", ylab="年齢軸すべて") axis(2) plot(x, tsrate[[1]], type="l", xlab="", ylab="Crude birth rate", cex.axis=1.3, cex.lab=1.6, lwd=3) par(new=TRUE) plot(x, tsrate[[2]], type="n", xlab="Year", ylab="", axes=FALSE, cex.lab=1.6) axis(4, cex.axis=1.5) lines(x, tsrate[[2]], col=4, lwd=3) mtext("Mortality rate", s=4, line=2.5, cex=1.6, col=4) abline(v=y, lty=2) plot(rowSums(pop), type="l", xlab="Year", ylab="Total population", lwd=3) points(x[y], rowSums(pop)[y], pch=16, col=4, cex=2) } }, interval=0.01, ani.width=600, ani.height=600)
