人口動態論

近年の日本では少子化が社会問題となっています。そこで， 今回は人口学 demography の問題を考えてみたいと思います。 人口の増加，あるいは個体数の増加と齢構成の関係はどうなっているのでしょうか。 それが人口学の問題なのです。

人口の行列モデル

それぞれの年齢における死亡率と出産率が与えられたとき， 人口の増加，もしくは減少が予測できるでしょうか。 ある年 t において年齢 x の雌個体の数を n_x,t とします。 x 歳での生存率を p_x とおけば x+1 歳のときの個体数は

と書けます。

また x 歳の雌の出産率を m_x で表わすことにすると， 次の年 t+1 での新生児の個体数 n_0,t+1 は次のように表わせます。

0 歳での出産率 m₀ を考えるのは奇妙な気がしますが， 一年性の草木や昆虫などのことも考えているので，このような事態が問題になるのです。 人間の場合には当然 0 歳での出産率は 0 ですねー。 だから m₀=0 となります。 人口学では雄の存在や役割は問題にしないことになっています( したがって雄である私は人口問題に関して何かをいえる立場にない T_T)。

ここで次の行列を考えてみましょう。

この行列は，1 行目にそれぞれの年齢での出産率， 2 行目以降の対角行列に生存率をおいた行列です。 ある年 t における雌の個体数 (n_0,t, n_1,t, n_2,t, …, n_w,t)^t を要素とする 縦ベクトルを n_tと書くと (1), (2)は

と書けます。ある年の雌の人口を n₀ とし， 繰り返しこの行列を適用することで x 年後の人口構成は，

と行列を x 回掛け合わせた，行列の累乗の形で表現できます。

この式の安定性は，行列の固有値と固有ベクトルによって表わされます。 実際，数値計算の分野では， 行列の固有値を求めるときに冪乗法といって行列を繰り返しかけ合わせることで， 固有ベクトルを求める方法が知られています。

安定年齢分布と指数増殖

従って，(4) を繰り返し用いて計算する場合，十分大きな t に対して， その年齢を持つ雌の個体比率が一定の値となり，年齢分布の形が収束することになります。 n_t = u とおくと， 次の年にはベクトルの方向はそのままで，長さだけが λ 倍になるので n_t+1 = λ u とおけます。 すると (4) は

となります。λ の満たすべき方程式は，

となります。すなわち固有ベクトルの性質から，u は， 方向だけに意味があるのであり，定数倍の自由度が残っています。 この収束先の年齢分布は安定齢分布と呼ばれています。(4)から， 1 年当たりの増加率が固有値 λ であり， 行列 M はフロベニウスの定理から， 複素数まで数えると w+1 個の根があることになります。 正の実数解はそれらの中で最大の絶対値を持ちます。

実習

適当な行列 M を作って上記のことを確かめてください。 matprd というプログラムは 行列の積を計算してくれます。 エディタなどで n 次の正方行列， M を作り ./matprd M M > M2 などとすれば M² が M2 というファイルに格納されます。 ./matprd M2 M > M3 とすれば M³ が M3 というファイルに作られます。 以下同様にして M10 くらいまで作って， 適当な n 行 1 列の行列 u を作って ./matprd M10 u として人口の動態を観察してみてください。 この場合 u の 10 年後の人口動態（各年齢の人口分布） が観察されることになります。

行列 M の固有値は simple-eigen というプログラムで求めることができます。 ./simple-eigen -p 1 M とすれば， 行列 M の最大固有値 1 つと 対応する固有ベクトルが表示されます。 最大固有値 λ が 1 以上であれば， ./matprd M10 u の結果は増大しているはずです。 なお u がどんな値であっても， 固有ベクトルの方向に近づくので，初期値 u はどんな値であっても良いです。 たとえば，すべての要素が 1 である u = (1, 1, …, 1)^T であっても，一定の比率（固有ベクトルの方向）に収束するはずです。 このことを確かめてみましょ。

たとえば，生存率 m が一定で m₀ = m₁ = … = m₈ = m = 0.96 で 0 代から 80 代までの出産率を m=2 の ポアソン分布に従うとしてみると，

としてみます。M の最大固有値は λ_max=1.11292 であるので 増加するはずです。matprdによって M¹⁰ を作り u = (1, 1, 1, 1, 1, 1, 1, 1, 1)^T を右からかけてみると

一方のこの M の出産率を少し下げて

としてみると M の最大固有値は λ_max=1.0003 と ほぼ 1 なので人口は安定しているはずです。 実際，matprd によって M¹⁰ を作って M¹⁰u を求めてみると，

となって，すべてを足し合わせてもほぼ 9 で安定していることが分かります。 両者を比較してみるためにグラフ化してみると

固有値 λ が 1 を超えている M においては出産率が高いために 0-10 代の人口がかなり増大していることが分かります。

人口構成の連続モデル

人口やサイズなどは実際には連続量となります。ある年齢 x の時刻 t における 分布は n(x, t) という分布関数によって表わされるでしょう。

ここで J(x,t) という関数を考えます。 これは x が増加する方向，図では 右方向への流れの量を表しているものとしましょう。 すなわち単位時間あたりに特性値 x がより小さい値から大きい値へと変わった個体数と，逆方向へ変化した個体数との差を表わします。 そうすると n(x,t) の時間変化は

という偏微分方程式で表わされます。 ∂n(x,t)/∂t は 2 変数関数 n(x,t) の t についての偏微分(変な微分じゃないよー，偏った微分のことね) と呼ばれています。もう一方の変数 x を定数のようにみなして t についてだけ微分したものです。

このときの時間変化は x+a からこの区間を出て行く流れと x から入ってくる 流れとの差し引きになります。よって

となります。この式を a で割って a を 0 に近づけると(13)を得ます。

死亡などで個体変動が起こらず，すべての個体が一定速度 v で変化する場合を考えると

となり，x の分布 n(x,t) は

となります。さらに v=1 とおけば，

となります。

一方，出産による新生児の加入速度 n(0, t) は， x 歳の女性が女の子を出産する 単位時間あたりの速度を m(x) として

として表わされます。

単位時間あたりの x 歳での死亡率を u(x) としてみましょう。 生まれてから x 歳まで 生存できる率 l(x) は

です。t > x とすると， 時刻 t において x 歳の個体は 時刻 t - x で生まれた。 すなわち

が成立するはずです。(20)を(18)に代入すれば，

であり，時刻 t での新生児の数 n(0,t) を

と表してみると

この式から，人口構成が収束した時の増加率 r を求めることができます。

は

に対応し(ここで，x 歳での死亡率を u(x) としました。 分布 n(x,t) がこの死亡率により u(x) n(x,t) の速度で減少するとして右辺第二項が加えられました），

(2)は

(18)は

に対応することがわかります。

また行列 M の固有方程式は

に対応し，(27)は

に対応します。最大固有値 λ が e^r に対応します。

生物の拡散過程

生物が離散的生息場所に棲んでいて， 単位時間あたりに場所 i から i+1 へ移動する個体の数を， 場所 i にいた個体数 n_i に比例した D n_i と書くことにします。 逆に場所 i+1 から i に移動した個体数を Dn_i+1 と書くことにしましょう。 この差のことを，個体群の流れといい，i の増加する方向を正として計ると

に代入すると

を得ます。これを拡散方程式 diffusion equation と呼んでいます。 時刻 t=0 において N 個体が場所 x=0で放たれ， ランダムな運動によって生息場所を広げていく過程を表していると考えられています。 この式の解は，

で与えられます。

実習

java DiffusionEquation

とすると拡散方程式の様子を描くプログラムが起動します。 数値をいろいろ変えて遊んでみてください。

拡散方程式の解き方

ランダムウォークの極限としてのブラウン運動

まず最初に，一次元のブラウン運動を考えてみます。これは ランダムウォーク， 酔歩の問題として知られているものです。 1 ステップの幅 a と 1 ステップ進むために要する時間 τ を 0 に近付けたときの極限を考えます。

と表現できます。 右に進んだ回数と，左に進んだ回数とをそれぞれ求めるためには，上式を s について解いて，

とすることができます。最右辺は n を t/τ で置き換えたものです。

左右に移動する確率を 1/2 で等しいとすれば n ステップ回で s 回右に行く確率は 2 項分布を用いて

として整理すると

上式第 1 因数は 1 です。 第 2 因数と第 4 因数を使って

などに注意して，n, s で表された (45) を x と t で書き換えると

を得ます。 (44)式第 3 因子については

を得ます。まとめると

座標 x₁ と x₂ との間， もしくは元の座標 s₁ と s₂ の間に粒子を発見する確率を問題にする場合

を計算しなければなりません。 連続性の仮定が成り立つと考えて，和を積分に書き換えて

として、最終的に

となります。これは

を確率密度（正規分布）と同一視することができます。