Rei Frontier Tech Blog

人工知能を活用した位置情報分析プラットフォーム「SilentLog Analytics」を運営する、レイ・フロンティア株式会社のエンジニアメンバーで運営する技術ブログです。

安定人口モデルについて

レイ・フロンティア株式会社でデータアナリストをしている齋藤です。これから、考えていることや興味を持って勉強していること、業務のなかで学んだ技術などについて定期的に発信していくことになりました。私個人としては数学や応用数理に興味を持って勉強を続けているので、そのような記事が中心になるかと思います。
お付き合い頂けると幸いです。

さて、現実の現象を理解するための手段として、現象を数理モデルによって記述し、モデルの性質を数学の問題として調べ、得られた結果を現実の言葉に翻訳するという数理的アプローチは基本的なものであり、今日ではあらゆる分野において数理モデルを用いた研究が行われています。今回はその中から、人口現象を説明するモデルである安定人口モデルを紹介します。

モデルの導出

歴史的経緯

人口学的現象を数理モデルによって説明しようとする試み自体は古くから存在しています。最も古いものは、おそらく13世紀にフィボナッチ(Leonardo Fibonacci)が彼の著書 "Liber Abaci"の中で紹介したフィボナッチ数列でしょう。
いま、ウサギのひとつがいを単位として、最初に1つがいのウサギが生まれてから\( n \)ヶ月後のつがいの数を\(P_n\)とします。また、

  • ウサギは生まれてから2ヶ月後に1つがいのウサギを産む。
  • ウサギが死ぬことは考えない。

とします。すると、漸化式 $$\left\{\begin{align} & P_{n+2} = P_{n+1} + P_n \\ & P_0 = 1 \end{align}\right.$$ が成り立ちます。\(n \to \infty\)とすると、比\(P_{n+1} / P_n\)は黄金比\(\phi = 1.6180339887...\)に収束します。これは、ウサギのつがい数は十分に時間が経つと公比\(\phi\)の等比数列に近づいていくことを意味します。このことを、つがい数は漸近的にMalthus径数\(\phi\)の指数的成長を行うと言い表します。後にも触れますが、人口が漸近的に指数的成長をすると仮定した場合の成長率を調べるというのは、人口モデルの解析における常套手段です。

人間の個体数の議論に数理モデルが広く受容されるようになったのは、マルサス(T. R. Malthus)が著書『人口論』において差分方程式モデル $$ P_{n+1} = (1 + r) P_n $$ を用いたことがきっかけでした。マルサスは、人口は指数関数的に増大するが食料は一次関数的にしか増大しないことから、人類は将来深刻な生活資源不足に陥ると警告しました。
人口モデルに年齢という概念がはじめて導入されたのは、ロトカ(A. J. Lotka)とシャープ(F. R. Sharpe)の論文([4])で提起された再生方程式(renewal equation) $$ B(t) = \int_0^\infty B(t-a) \ell(a) \beta(a)\, da$$ でした。ここで、\(B(t)\)は時刻\(t\)における単位時間あたりの出生数、\(\ell(a)\)は\(a\)歳までの生残率、\(\beta(a)\)は\(a\)歳の人口の年齢別出生率です。ロトカらは、再生方程式の解が漸近的に指数的成長をすることと、その漸近的な成長率が\(\ell\)と\(\beta\)によって決定されることを(数学的な厳密さは欠いていたものの)示しました。すなわち、外部との人口の流入出がない封鎖人口において、年齢ごとの死亡率と出生率が時間によらず一定ならば、時間が経つにつれて人口構造(ピラミッド型・壺型の図でおなじみの、年齢ごとの人口の比率のことです)は一定となり、総人口は指数的に変化するようになることを発見したのです。この理論は安定人口理論(stable population theory)と呼ばれ、人口学における基本的理論のひとつとなっています。彼らの結果はフェラー(W. Feller)によって数学的に正当化されました([1])が、その証明はLaplace変換を利用した極めて技巧的なものであり、それ以上の発展は到底見込めませんでした。それから40年以上経って、ウェブ(G. F. Webb)によってより一般の非線形人口モデルについて、それに対応する非線形半群による解を構成し、その振る舞いを調べるという関数解析的手法が打ち立てられました([5])。その後、メッツ(J. A. J. Metz)やディークマン(O. Diekmann)らによって年齢だけでなく位置やサイズなど様々な構造をもつモデルが研究され(その成果は論文集 [3] にまとめられています)、構造化人口動態学(structured population dynamics)という応用数学の一分野が確立されるに至ります。本記事で紹介するのは、構造化人口動態学の枠組みで考えられる最も単純なモデルです。あくまで紹介記事ですので、定理の証明などを細かく追うことはしません。きちんとした数学的な議論を見たい方は、たとえば[2], [3], [5]を参照するとよいでしょう。

パラメータの定義

まず、モデルに登場するパラメータたちについて説明します。
年齢構造をもつ人口集団において、時刻\(t\)における年齢密度関数(age-density function)を\(p(t,a)\)と書きます。積分$$ \int_{a_1}^{a_2} p(t,a)\, da $$は時刻\(t\)における年齢\(a_1\)以上\(a_2\)以下の人口になります。本記事で紹介するモデルは、\(p(t,a)\)に関する初期値境界値問題として定式化されます。
関数\(\mu (a)\)を、年齢\(a\)における個体の死亡の瞬間的な発生率とし、死亡力(force of mortality)あるいは年齢別死亡率(age-specific death rate)と呼びます。積分 $$\int_{a_1}^{a_2} \mu(a) p(t,a) \, da$$は単位時間あたりの年齢\(a_1\)以上\(a_2\)以下の死亡数となります。確率変数\(X\)を個体の死亡年齢とし、\(X\)の分布関数を\(F(a)\), 密度関数を\(f(a)\)とすると、\(\mu (a)\)は $$ \begin{align}\mu (a) &= \lim_{h \to 0} \frac{\mathbb{P}(a \leq X < a + h \mid a \leq X)}{h} = \lim_{h \to 0} \frac{\mathbb{P}(a \leq X < a + h)}{\mathbb{P}(a \leq X) h} \\ &= \frac{f(a)}{1 - F(a)} = -\frac{d}{da} (1 - F(a)) \end{align}$$ です(ここで、\(\mathbb{P}(A)\)は事象\(A\)がおこる確率を表します)。 このとき、 $$ \ell (a) := 1 - F(a) $$は個体が\(a\)歳まで生存する確率を意味し、生残率(survival rate)と呼ばれます。\(\mu(a)\)の定義から、微分方程式 $$ \frac{d\ell (a)}{da} = -\mu(a)\ell(a) $$が成り立ちます。これと\(\ell(0)\)(0歳の時点では確率1で生き残っている)から、$$ \ell (a) = \exp \left( - \int_0^a \mu (b) \, db \right) $$を得ます。いま、各個体は有限の最大年齢\(\omega\)をもつと仮定しましょう。すると\(\ell(\omega) = 0\)でなければならないので、$$\int_0^\omega \mu(a)\, da = \infty $$が必要です。たいていの場合\(\mu\)は上の式を満たす局所可積分な非負関数と仮定されますが、実際のデータでは年齢が\(0, \omega\)に近いところで急激に増加するU字型のグラフを示します。
関数\(\beta(a)\)を、単位時間あたりに微小な年齢区間\([a, a + da]\)にいる一個体から産まれる新生児の数とし、これを年齢別出生率(age-specific birth rate)と呼びます。積分$$\int_{a_1}^{a_2} \beta(a) p(t,a) \, da$$は単位時間に年齢\(a_1\)以上\(a_2\)以下の個体から産まれる新生児の数です。単位時間に産まれる新生児の総数\(B(t)\)は、$$B(t) = \int_0^\omega \beta(a) p(t,a) \, da$$によって与えられます。再生方程式は、この\(B(t)\)についての積分方程式です。実際のモデル解析においては、\(\beta\)はある有界閉区間\([\beta_1,\beta_2]\)の外で0である非負有界可積分関数と仮定されるのが通例です。区間\([\beta_1,\beta_2]\)は子どもを産める年齢の区間を意味し、再生産年齢(reproduction age)と呼ばれます。

方程式の導出

時刻\(t\)における\(a\)歳以下の人口を\(F(t,a)\)とします。すなわち、$$F(t,a) := \int_0^a p(t,s)\, ds$$とします。すると、時間間隔\([t, t+h]\)における人口の増減を考えることにより、人口方程式
$$ \begin{align}F(t+h,a+h) - F(t,a) &= ( [t, t+h] で生まれた人数 ) - ( [t, t+h] で亡くなった人数 ) \\ &= \int_t^{t+h} B(s)\, ds - \int_0^h \!\!\! \int_0^{a+s} \mu(b) p(t+s,b)\, dbds \end{align}$$を得ます。両辺を\(h\)について微分して\(h=0\)とおくと、積分記号下での微分を形式的に用いて、連続方程式$$p(t,a) + \int_0^a \frac{\partial p(t,s)}{\partial t}\, ds = B(t) - \int_0^a \mu(t,s)p(t,s)\, ds$$が得られます。特に\(a=0\)とおけば、$$p(0,t) = B(t) = \int_0^\omega \beta(a)p(t,a)\, da$$を得ます。連続方程式を\(a\)について微分すれば、$$\frac{\partial p(t,a)}{\partial t} + \frac{\partial p(t,a)}{\partial a} = -\mu(a) p(t,a)$$となります。これをマッケンドリック方程式(McKendrick equation)と呼びます。初期条件を\(p(a,0) = p_0(a)\)と設定することにより、年齢密度関数\(p(t,a)\)は以下の初期値境界値問題にしたがいます:$$\left\{ \begin{align} & \frac{\partial p(t,a)}{\partial t} + \frac{\partial p(t,a)}{\partial a} = -\mu(a) p(t,a), \quad t>0,\ a>0 \\ & p(0,t) = \int_0^\omega \beta(a)p(t,a)\, da,\quad t>0 \\ & p(0,a) = p_0(a),\quad a \in [0,\omega] \end{align} \right.$$これが安定人口モデル(stable population model)と呼ばれるものです。

モデルの解析

解の存在と一意性、人口作用素の性質

以下では、1980年代半ばに大きく発展を遂げた安定人口モデルへの関数解析的アプローチを概観します。安定人口モデルを\(L^1\)における抽象的な微分方程式として捉え直せば、Hille-Yosidaの半群の生成定理を用いて解作用素を導くことができます。本記事では線形の安定人口モデルのみを扱っていますが、この手法が真に効力を発揮するのは、モデルの右辺がより複雑な非線形問題を考えるときです。実際、1980年代以降の構造化人口動態学の発展は半群の理論を駆使した関数解析的手法の開発によるところが大きいのです。
いま、安定人口モデルにおける死亡力と年齢別出生率は\(\mu, \beta \in L_+^\infty (0,\omega) \)をみたすと仮定します。すなわち、$$ \bar\mu := \sup_{a \in [0,\omega]} \mu(a) < \infty , \ \bar\beta := \sup_{a \in [0,\omega]} \beta(a) < \infty $$が成り立つとします。このとき、安定人口モデルは次のようなBanach空間\(L^1(0,\omega)\)上の抽象的Cauchy問題として考えることができます: $$\left\{ \begin{align} & \frac{dp(t)}{dt} = Ap(t) \\ & p(0) = p_0 \end{align} \right.$$ ここで、\(A\)は以下のように定義される人口作用素(population operator)です: $$ \left\{ \begin{align} & (A\phi)(a) = -\frac{d\phi (a)}{da} - \mu(a)\phi(a) \\ & \mathcal{D}(A) = \left\{ \phi \in L^1(0,\omega) \,\middle| \, \phi は絶対連続かつ \phi(0) = \int_0^\omega \beta(a) \phi(a) \, da \right\} \end{align} \right. $$ \(\phi \in \mathcal{D}(A)\)ならば、ほとんどいたるところ微分係数\(\phi^\prime\)が存在して\(\phi^\prime \in L^1(0,\omega)\)となります。
人口作用素は以下の性質をもちます:

定理1.
人口作用素\(A\)に関して、以下の(1)-(6)が成り立つ。ただし、$$\hat\Psi := \int_0^\omega e^{-\lambda a} \beta(a) \ell(a)\, da , \quad \underline{\mu} := \inf_{a \in [0,\omega] \mu(a)}$$と定義しておく。
(1) \( \sigma(A) = P_\sigma(A) = \{ \lambda \in \mathbb{C} \mid \hat\Psi(\lambda) = 1 \}\) である。ここで、\(\sigma(A), P_\sigma(A)\)はそれぞれ\(A\)のスペクトル集合、点スペクトルを意味する。
(2) レゾルベント \(R(\lambda, A)\)はコンパクト作用素であって、$$\begin{align}(R(\lambda,A)\phi)(a) &= \int_0^a e^{-\lambda (a-z)} \frac{\ell(a)}{\ell(z)}\phi(z)\, dz \\ & + \frac{e^{-\lambda a}\ell(a)}{1-\hat\Psi(\lambda)} \int_0^\omega \!\!\! \int_z^\omega e^{-\lambda (a-z)}\frac{\ell(a)}{\ell(z)} \beta(a)\, da \phi(z) \, dz \end{align}$$ (3)\(\lambda \in \sigma(A)\)であれば\(\lambda\)の代数的重複度は有限であって、\(\lambda\)に対応する固有空間は\(N(\lambda-A) = \{ce^{-\lambda a} \mid c \in \mathbb{C} \}\)で与えられる。
(4) \( \rho(A) \supset \{\lambda \in \mathbb{C} \mid \Re \lambda < \bar\beta - \underline{\mu}\}\)であり、\(\lambda \in \mathbb{R}\)で\(\lambda>\bar\beta-\underline{\mu}\)であれば\(R(\lambda,A)\)は正値作用素である。
(5) \(\lambda\in\sigma(A)\)ならば\(\bar\lambda\in\sigma(A)\)である。
(6) \( ^\forall \alpha > -\infty \)について、\(A\)は半平面\(\Re \lambda \geq \alpha\)に高々有限個の固有値しかもたない。

定理2.
人口作用素\(A\)は閉作用素であり、かつ\(\mathcal{D}(A)\)は\(L^1(0,\omega)\)で稠密である。

以上の準備のもとで、Hille-Yosidaの定理により次が成り立ちます:

定理3.
人口作用素\(A\)は\(C_0\)半群 \(T(t)\ (t\leq 0) \)を生成する。\(T(t)\)は以下をみたす:$$ \begin{align} & \|T(t)\| \leq e^{(\bar\beta - \underline{\mu}) t} \\ & T(t)(L_+^1) \subset L_+^1 \end{align}$$しかも、\(p_0 \in \mathcal{D}(A)\)ならば\(p(t)=T(t)p_0\)は抽象的Cauchy問題の古典解を与える。

人口作用素のスペクトルとSharpe-Lotka-Fellerの定理

\(C_0\)半群とその無限小生成作用素のスペクトルに関する一般的な性質について少し述べます。作用素\(A\)のスペクトル上限\(s(A)\)とスペクトル半径\(r(A)\)はそれぞれ $$\begin{align} & s(A) := \sup\{ \Re \lambda \mid \lambda \in \sigma(A) \} \\ & r(A) := \{ |\lambda| \: \mid \lambda \in \sigma(A) \} \end{align}$$ また、半群\(T(t)\)の成長上限(growth bound)は$$\omega_0 := \inf \{w \in \mathbb{R} \mid ^\forall t\leq 0,\ ^\exists M(w)\in \mathbb{R}_+,\ s.t.\ \|T(t)\| \leq M(w)e^{wt} \}$$ と定められます。このとき、$$\omega_0 = \lim_{t\to \infty} \frac{\log\|T(t)\|}{t}$$であり、\(-\infty \leq \omega_0 < \infty\)が成り立ちます。

定理4.
\(\omega_0\)を\(C_0\)半群\(T(t)(t\geq 0)\)の成長上限、\(A\)を無限小生成作用素とすれば、\(^\forall t\geq 0\)に対して以下が成り立つ:$$e^{s(A) t}\leq r(T(t)) = e^{\omega_0 t}$$ \(A\)が有界作用素であれば、$$e^{s(A) t} = r(T(t)) = e^{\omega_0 t}$$とくに、\(\omega_0 = s(A)\)である。

\(\omega_0 = s(A)\)であるならば、無限小生成作用素のスペクトルの実部の上限が負であれば、\(\lim_{t\to \infty}\|T(t)\| = 0\)となります。これは、連立微分方程式\(x^\prime (t) = Ax(t)\)の解の漸近挙動を調べるには、行列\(A\)の固有値の実部の正負に注目すればよい……という議論の無限次元バージョンとなっています。
\(t_0\)が存在して\(T(t_0)\)がコンパクト作用素になるとき、半群\(T(t)\)は終局的コンパクト(eventually compact)であるといいます。終局的コンパクトな半群に関して、以下の定理が成り立ちます:

定理5. (スペクトル写像定理) $$\sigma(T(t)) \setminus \{0\} = e^{t\sigma(A)} := \{e^{t\lambda} \mid \lambda \in \sigma(A)\} , \quad t\geq 0 .$$ 定理3によって存在が保証された半群は人口半群(population semigroup)と呼ばれますが、これは終局的コンパクトであることが知られています。よってスペクトル写像定理が成り立ち、そのとき\(\omega_0 = s(A)\)なので、抽象的Cauchy問題のゼロ解が大域的安定である(すなわち、\(t\to\infty\)とすると人口がゼロに近づいていく)ためには、人口作用素\(A\)のスペクトル上限が負であれば十分であることがわかります。一般に、そうなっているかどうかは\(\mu\)と\(\beta\)のバランスによって決定されます。
ゼロ解の大域安定性について、より精密な条件を与えるのがSharpe-Lotka-Fellerの定理です。この定理によって、特性方程式\(\hat\Psi(\lambda) = 1\)のただ一つの実根\(\lambda_0\)によって安定人口モデルの漸近挙動が支配されることが判明するのです。ここでは、関数解析的な言葉を使って記しておきます:

定理6.
\(\lambda_0\)を人口作用素の実固有値とする。このとき、1次元の射影子\(P_0\)と正数\(\varepsilon >0, M(\varepsilon)\geq 1\)が存在して、$$ \| e^{-\lambda_0 t} T(t) - P_0\| \leq M(\varepsilon)e^{-\epsilon t}$$となる。ここで、\(u_0,v_0\)を\(A\)の固有値\(\lambda_0\)に対応する固有ベクトルおよび共役固有ベクトルとすれば、$$P_0\psi = \frac{\langle v_0,\psi \rangle}{\langle v_0,u_0 \rangle}u_0$$と表される。

つまり、半群\(T(t)\)は初期値によって決まる定ベクトルの\(e^{\lambda_0 t}\)倍に近づいていきます。\(\lambda_0\)は自然成長率(intrinsic rate of natural increase)あるいは安定人口成長率、内的成長率などと呼ばれています。\(\lambda_0<0\)ならば人口はゼロ解に近づいていき、\(\lambda_0>0\)ならば増大していきます。

Lotkaの再生方程式と基本再生産数

安定人口モデルを特性線\(t-a=Const.\)に沿って積分することにより、以下の表示を得ます:$$p(t,a) = \left\{ \begin{array}[ll] \, B(t-a) \ell(a), & t-a>0 \\ p_0(a-t) \frac{\ell(a)}{\ell(a-t)}, & a-t>0 \end{array} \right.$$したがって、出生数\(B(t)\)が求まれば年齢密度関数も求まります。これを境界条件に代入することにより、以下の積分方程式を得ます:$$B(t) = G(t) + \int_0^t \Psi(a)B(t-a)\, da$$ ここで、\(G(t), \Psi(a)\)は以下によって与えられる既知関数です:$$G(t)= \int_t^\omega \beta(a) \frac{\ell(a)}{\ell(a-t)}p_0(a-t)\, da , \quad \Psi(a) = \beta(a)\ell(a)$$特に、\(G(t)\)は時刻\(0\)に生存している個体から時刻\(t\)に産まれる単位時間あたりの新生児数と解釈されます。これが、Lotkaの積分方程式、あるいは再生方程式と呼ばれているものです。定理6から、\(B(t)\)に関するSharpe-Lotka-Fellerの定理が帰結します:

定理7. (Sharpe-Lotka-Fellerの定理)
安定人口モデルにおける再生方程式の解\(B(t)\)に関して、以下が成り立つ:$$B(t) = q_0 e^{\lambda_0 t}(1+O(e^{\eta t}))$$ここで\(\eta\)はある正数であって、$$q_0 = \frac{\int_0^\infty e^{-\lambda_0 t}G(t)\, dt}{\int_0^\infty a e^{-\lambda_0 a}\Psi(a)\, da}$$である。

さて、関数\(B_n(t)\ (n=0,1,2,\dots)\)を次のように定義してみましょう:$$\left\{ \begin{align} &B_0(t) = G(t)\\ &B_n(t) = \int_0^t \Psi(t-a) B_{n-1}(a)\, da ,\quad n\geq 1 \end{align} \right.$$関数\(B_n(t)\)は、初期人口の第\(n\)世代の子孫の出生率を与えています。各世代にわたる和\(\sum_{n=0}^\infty B_n(t)\)はLotkaの積分方程式の解となっています。
各世代の総出生数を計算してみましょう。漸化式から、$$\begin{align} \int_0^\infty B_n(t)\, dt &= \int_0^\infty \, dt \int_0^t \Psi(t-a) B_{n-1}(a)\, da \\ &= \int_0^\infty B_{n-1}(a)\, da \int_a^\infty \Psi(t-a) \, dt = \mathcal{R}_0 \int_0^\infty B_{n-1}(t)\, dt \end{align}$$ を得ます。ここで、$$\mathcal{R}_0 = \int_0^\omega \Psi(a) \, da$$は基本再生産数(basic reproduction number)といい、一個体が生涯に産む個体数と解釈される量です。人口学的な意味から、\(\mathcal{R}_0>1\)ならば人口は増大し、\(\mathcal{R}_0 <1\)ならば人口は減少すると予想されますが、実際、自然成長率は特性方程式$$\int_0^\omega e^{-\lambda_0 a}\Psi(a)\, da = 1$$をみたす実数であったので、符号関係$$\mathrm{sgn}(\mathcal{R}_0-1) = \mathrm{sgn} \lambda_0$$が成立します。基本再生産数は安定人口モデルだけではなく、その拡張のモデルや感染症モデルなどの様々な個体群動態モデルに現れており、構造化人口動態学における最も重要な概念の一つとして知られています。

参考文献

[1] W. Feller (1941), "On the integral equation of renewal theory", Ann. Math. Stat. 12: 243-267.

[2] 稲葉 寿 (2002), "数理人口学", 東京大学出版会.

[3] J. A. J. Metz, O. Diekmann (1986), "The Dynamics of Physiologically Structured Populations", Lect. Notes in Biomath. 68, Springer, Berlin.

[4] F. R. Sharpe, A. J. Lotka (1911), "A problem in age-distribution", Philosophical Magazine, Series 6, Vol. 21: 435-438.

[5] G. F. Webb (1985), "Theory of Nonlinear Age-Dependent Population Dynamics", Marcel Dekker, New York and Basel.