日本の人口をデータから見た👀

みなさん、こんばんは！若葉のマフィンです！

今日は少し天気がよくなかったですね...こんなとき私は気分がどうしても沈んでしまいます💦

ですが、あまりに元気がないのもあれなんでデータ解析で元気出していきましょー！

というわけで、今回は前回の記事

wakaba-mafin.hatenablog.com

で手元のChromebookに落としてきたデータを簡単に解析してみます！ただPythonのライブラリ一発でデータ解析しましたという記事ではなく、個体群の数理モデルのお勉強もかねまして、データの解析の一連の手順をなぞって今後の日本の人口がどうなるのか（増えるのか/減るのか）について考えたいと思います！

（※この記事の内容をご自身のPCで試す場合は自己責任でお願いします）

1種類の個体群の人口動態モデル
データの準備
モデルの推定
- Malthusian Model
- Logistic Model
考察
まとめ

1種類の個体群の人口動態モデル

まずは人口動態に関する数理的な知識を確認します！今回はこちらの本を参考にしました。

[丸善出版ニュース] 「マレー数理生物学入門」* 数理生物学テキストとして世界的な権威書。* 2014年1月発売

この本では今回のデータ解析のような単純な個体群の成長モデルから、弱肉強食のように多種の個体群が相互作用しながらどのように個体数が変わるのかや感染症がどのように広がっていくのかまで様々な数理モデルが記述されています。もちろん、数学的な解析も一部含まれはしますが、まさに数理生物の入門書にいい一冊だと思います！

特に今回はこちらの書籍の紹介ではないので、早速人口動態のモデルに触れていきましょう！まず、人口 $N$ のダイナミクスをモデルで表す上で最も重要な式があるので、まずはそちらを確認しましょう。それは

$\frac{dN}{dt} = \mbox{出生} - \mbox{死亡} + (\mbox{移住})$

です。これは単位時間あたりに人口が出生で増えて死亡で減って、そして移住によって変わることを表しています。これが最も基本的な式で、出生・死亡・移住をどうやってモデリングするのかによって様々なバリエーションが出てきます。

人口のダイナミクスを表すモデルで最も簡単で有名なものは1978年にMalthusさんが作ったモデルです。これは定数 $r=$ （単位時間あたり出生率） $-$ （単位時間あたり死亡率）を用いて

$\frac{dN}{dt} = rN$

で表されるモデルです！昔（文明が発展している途中）の人口動態を表すにはこのモデルでもよかったらしいのですが、このモデルは基本的に無限に増加し続けるか、絶滅するかの２つのパターンしか表せません💦つまり長期的には絶対に予測できないという問題点があります。これを改善するための１つのモデルが1838年にVerhulstによって提案された人口の抑制効果も考えたロジスティックモデルです。

$\frac{dN}{dt} = rN ( 1 - \frac{N}{K}) = rN -\frac{r}{K}N^{2}$

ここで $r$ は内的自然増加率といって抑制が働かない場合の増加率で、 $K$ は環境収容力といって食べ物の量など環境によってどれくらいの人口を保つことができるかを表す定数です。つまり、人口が $K$ に比べて小さい場合には $dN/dt \approx rN$ となって指数関数的に増加し、ある程度増えてくると二項目の二乗の引き算の項が効果を発揮して人口の増加が抑制されます。これはまさに人口が増えすぎると食糧難などが起こり、それ以上人口が増えにくくなることを表していると思ってもらえればいいと思います！

以上、人口動態のモデルとして最も基本的な２つのモデル

マルサスモデル
ロジスティックモデル

を紹介しました。ここからこれらのモデルを使ってe-Statから持ってきたデータを解析しましょう！

データの準備

まずはデータの確認です。データはすべてe-Stat（政府統計の総合窓口）からダウンロードしてきます。詳しくは

wakaba-mafin.hatenablog.com

まずはライブラリ

import scipy as sp
from scipy import optimize
import matplotlib.pyplot as plt
import pandas as pd

今回はこれくらい使えばOKです！そしてデータのロードとplot

# load
df = pd.read_csv("population.csv") # ファイル名は以前の記事参照
gender = df[" 性別 "].values
year = df[" time_code "].values[gender==" 総数 "][::-1]/1e+6
N = df[" value "].values[gender==" 総数 "][::-1]
N0 = N[0]
# plot
plt.plot(year, N, "k-", label="population")
plt.xlabel("year")
plt.legend()
plt.show()

f:id:wakaba-mafin:20181109230616p:plain — fig.1 e-Statからダウンロードした人口動態

ここまでは前回の記事とほとんど同じ。

そして、データ解析には必須の前処理。今回は1939年から1950年あたりのデータを使わないようにします。上のグラフを見てもわかるのですが...少し人口動態の質が違うというか...理由はご自身で考えてみてください！

# preprocess
dt = 1 # sampling period
target_loc = (year<1939) + (year>1950)
year_edit = year[target_loc]
dNdt_edit = sp.diff(N)[target_loc[1:]]/dt
N_edit = N[target_loc]

これで省きたいところを省けました！（ついでに微分の値も差分によって計算しています💦）

f:id:wakaba-mafin:20181109230812p:plain — fig.2 省きたい部分を除去したデータ

こんな感じで。というわけで、データの準備ができました！

モデルの推定

ここからは上で紹介したマルサスモデルとロジスティックモデルを日本の人口動態に対して推定します！モデルの推定の方法としては２つのアプローチを取ります！

アプローチ１：システムノイズの仮定
- 説明変数： $t$ 年における人口 $N(t)$
- 従属変数： $t$ 年における人口増加率 $\frac{N(t+1)-N(t)}{1\mbox{年}}$ （※データはすべて一年おきでしかないので差分を取るだけで実はOK）
アプローチ２：観測ノイズの仮定
- 説明変数：年次 $t$
- 従属変数： $t$ 年における人口

の２つです。前者は人口のダイナミクスの中にノイズが入っていると考え、データにおけるその分散、すなわち二乗誤差を最小にする推定です。逆に後者は人口のダイナミクス自体にはノイズなどは入っておらず解析的にモデルを解いた後で、人口データを観測する際にノイズが入っていると考えて、その二乗誤差を最小にする推定です。後者の方が見た目上フィッティングできているように見えるはずです...💦

Malthusian Model

では早速、日本の人口データに対してマルサスモデルを推定してみましょう！

# Malthusian Model
## approach 1
def deriv_MM(x, r):
    return r*x

dMMparam, dMMcov = sp.optimize.curve_fit(deriv_MM, N_edit[:-1], dNdt_edit, bounds=[0,1]) # 推定

## approach 2
def int_MM(t, r):
    return N_edit[0]*sp.exp(r*(t-year[0]))

iMMparam, iMMcov = sp.optimize.curve_fit(int_MM, year_edit, N_edit, bounds=[0,1]) # 推定

するとそれぞれのアプローチで自然増加率が

アプローチ１： $0.00608526\pm0.00055964$
アプローチ２： $0.01066609\pm0.00011088$

と推定できました。

そして、実際にどれくらいデータに合っているか確認すると...

f:id:wakaba-mafin:20181109232821p:plain — Fig.3 マルサスモデル推定結果

全然ですね笑　特にアプローチ１なんか全く外れてますね💦 （一応、青い点線で1900年から1980年までのデータだけ使って推定した結果も載せておきまーす！）

Logistic Model

まぁまぁ、気を落とさずにロジスティックモデルもやりましょう！同じようにすればOKです。

# Logistic Model
## approach 1
def deriv_LM(x, r, K):
    return dt*r*x*(1-x/K)

dLMparam, dLMcov = sp.optimize.curve_fit(deriv_LM, N_edit[:-1], dNdt_edit, bounds=[0,(1, 2*N.max())]) # 推定

## approach 2
def int_LM(t, r, K):
    return N_edit[0]*K*sp.exp(r*(t-year_edit[0]))/(K+N_edit[0]*(sp.exp(r*(t-year_edit[0]))-1))

iLMparam, iLMcov = sp.optimize.curve_fit(int_LM, year_edit, N_edit, bounds=[0,(1, 2*N.max())]) # 推定

すると、内的自然増加率と環境収容力は

アプローチ１
- 内的自然増加率： $0.0263658963\pm0.00180462122$
- 環境収容力： $1.40852585\times 10^8 \pm 3.45960046\times 10^6$
アプローチ２
- 内的自然増加率： $0.0224541033\pm0.00056251535$
- 環境収容力： $1.59150353\times 10^8 \pm 3.30948118\times 10^6$

そして、実際にデータに合ってるか確認すると...

f:id:wakaba-mafin:20181109233947p:plain — Fig4. ロジスティックモデルの推定結果

こちらは見た目上はかなり合ってるように見えますね！

考察

というわけで考察です。

まずマルサスモデルに関してですが、現在の日本の人口が増加から減少に辿っているという事実を考えると、モデル自体に問題があるためうまく推定できていないということが考えられますねー💦というのも、1980年以降人口はほぼ横ばいになっていますが、このような挙動をマルサスモデルでは取ることができないからですね（マルサスモデル触る必要なかったか...💦）ですが、1980年までのデータのみを使って推定した結果を見ると、かなりよくフィッティングできていることがわかります。このことから1980年代くらいまではマルサス的な人口増加をしてきたと考えられます！

で、続いてロジスティックモデルですが、こちらはどちらのアプローチでもほとんど同じ推定結果に落ち着き、見た目もかなりフィッティングできています。さらに内的自然増加率の推定値を見てみても、0.02程度。ここで内的自然増加率は一年あたりの一人あたりの出産数ですから、ここに性別の存在と女性が子供を産むことができると考えられる年齢15歳〜49歳を考慮に入れればー

$0.02[\mbox{年}^{-1}] \times 2 \times (49-15)[\mbox{年}] = 1.36$

日本の近年の合計特殊出生率と近くなって割と妥当っぽい！（？） ですが、環境収容力がそれぞれ1.4億人以上になってますねー...この単純なモデルでは人口まだ増えるみたいです笑

まぁ、ここでは経済状況やその他諸々の要素を完全に切り捨ててますので、こんなもんすね！笑

政府の方はさすがにすごいです！モデル立てて必要なパラメータは推定して、エビデンスがないパラメータや値については複数のシナリオを作ってシミュレーションをしてますね！

https://www.mhlw.go.jp/file/05-Shingikai-12601000-Seisakutoukatsukan-Sanjikanshitsu_Shakaihoshoutantou/0000138825.pdf

まとめ

というわけで、今回はe-Statの人口動態のデータを数理的な観点から解析してみました！

特に今回は

マルサスモデル
ロジスティックモデル

の２つのモデルをデータに対してフィッティングし、なんかちゃっちい考察にはなりましたが一応考察も与えました！これからの人口どうなるんでしょうね？経済に与える影響も心配です...💦

間違いのご指摘やご意見等ございましたら、ぜひコメントの方よろしくお願いします！

ではでは、長くなってしまいましたが今回はこれにて！

次回以降もよろしくお願いしまーす！

若葉の技術メモ

若葉の技術メモ

コンピュータやプログラミング・数理に関して調べたり、取り組んだりしたことをまとめる若葉のノート。

コンピュータ・プログラミング・数理が初めてって方も一緒に勉強していきましょう！

初心は大事！いつでも若葉☘のような意気込みで！

日本の人口をデータからじっくりと考えてみた...

1種類の個体群の人口動態モデル

データの準備

モデルの推定

Malthusian Model

Logistic Model

考察

まとめ