Last updated on December 25, 2002 (WED) 23:49
このページは作成途中です。また,数理生物学懇話会ニューズレター第26号(1998年9月)に書いた記事をベースにしています。 (This page is under construction. And, the content of this page is based on the article in the Newsletter No. 26 (Sept. 1998) of the Japanese Association of Mathematical Biology.) Sorry, mostly written in Japanese.
シミュレーションとは,日本語でいえば模擬実験,つまり,ある事象を実際にやる代わりにモデル(模型)で実験してみることをいう。中でも人口学では,ふつう,人口モデルの数値計算による予測あるいは過去の推計をさす。これらは,メカニズムよりも数値間に現れる関係式に基づいてアプリオリに立てたものが多く,観測不可能な変数も多い。この推計においてもっとも重視されるアウトプットは,幅はついているにせよ「推計値」である。
しかし,考えてみれば,アプリオリに立てたモデルのアウトプットというものは,正しい(あるいは現実に即している)という保証がない。モデルが正しいかどうかを評価する手段がないのに推計値に意味をもたせるのは危険である。そうすると,人口学におけるシミュレーションには意味がないのだろうか?
この方面の先駆者であるBennett Dykeは,次の2つの用法でシミュレーションを使うなら意味があるといっている。一つは,シミュレーションによるデータの生成である。本質的に観測不可能な事象でも,シミュレーションによってデータを生成すれば,完全なデータセットを得ることが可能である。
たとえば,聞き取りによって200年以上にわたる個人レベルの詳細な人口動態を完全に復元することは不可能である。したがって,信頼性のないデータが解析モデルにどういうバイアスを与えるかを評価する必要がでてくる。こういう場合にシミュレーションを使えば,出生,死亡などについていろいろなレベルを設定した上で,完全な人口動態データを短時間のうちに得ることができる。この完全なデータに対して解析モデルを適用して得られる結果と,何らかの脱落を仮定して得られるデータに対して同じモデルを適用して得られる結果がどの程度一致するかによって,聞き取りの不完全さに対する解析モデルの頑健性を評価できるわけである。さらに,乱数の初期値だけを変えて何度もシミュレーションを実行すれば,与えた初期条件の下でのデータの分布を得ることができる。これは一般にコンピュータ集約型統計学と呼ばれる方法で,解析的に求めにくい分布関数を,シミュレーションならば容易に求めることができる。
もう一つは,単純な解析的モデルが存在しないような,複雑な現象についてシミュレーションを用いることである。たとえば,結婚過程はそうである。なるほど,稲葉(1997)が指摘するとおり,Fredrickson (1971)以来,男女のペア形成の分布を考慮した非線形微分方程式モデルはいくつか提案されているが,たとえば,交叉イトコ婚を数学的に表現するには,まず交叉イトコという親族関係の変数を作り,年齢構造とは別にモデルに組み込んで,親族関係ごとに関数を定義しなければならないわけで,この手口ではきわめて難しい。その点,個人レベルのシミュレーションならばルールベースという形で結婚を表現することは容易である。他の方法ではまったく解が得られないなら,シミュレーションの利用は大いに価値があるのだ。
表題の「シミュレーション人口学」はぼくの造語であり,上記2つの目的でおこなうシミュレーションを含む人口分析をさすものとご理解いただきたい。ペア形成を含めて人口再生産のシミュレーションを個人レベルでおこなうことは,2つの目的のどちらにも利用でき,シミュレーション人口学の基本である。
個人レベルで人口再生産を表現するにはどうしたらよいだろうか? まず,シミュレーションの単位となる「個人」とはいかなるものか(どういう属性をもつか),決める必要がある。次に,人口再生産を表現するには,基本的に,出生と死亡の発生規則を決めなければならない。死亡は個人単位で決まるイベントと考えれば,死亡の発生確率は個人レベルで決めることができる(年齢の関数とすることが多い)。出生はペア単位で起こるイベントなので,ペア形成のルールを規定することと,ペアの出生力を決めることの両方が必要である。これら4段階を規定した後で,組み合わせることによって集団全体の人口再生産がおこなわれる。ただし,Goldmanらの一連の研究が既婚未婚が死亡に影響を与えることを示唆しているので(たとえばHu and Goldman., 1990),死亡を個人レベルのイベントとするこの枠組み自体も検討の余地を残している。
いいかえれば,ヒトの属性のうち,どの部分をピックアップするかということである。最低限必要なのは,個体識別子,親子関係,年齢(月齢,日齢),配偶関係であるが,居住地や疾病抵抗性遺伝子といった属性を扱うこともできる。どの属性を要素に含めるかは,死亡やペア形成のルールベースを決める際の必要性から決定する。C言語では構造体として定義するのが自然である。
個人レベルでの死亡発生確率を決めるには,大きく分けて4通りの方法がある。モデル生命表を使う方法,死亡曲線のモデル(Silerモデルなど)を使う方法,Brassのロジットシステムを使う方法,及び,雪崩モデルである。個人単位の人口シミュレーションには雪崩モデルがもっとも適していると思う(ただし乳幼児死亡は別に考えなくてはならないが)。
これを扱えることがシミュレーションの最大の利点である。たとえば,年齢差何歳未満で,異なるクランに属していて,配偶者が何人以内といった条件を満たす配偶者を探して,候補がいればその中からある確率で結婚する,という複雑な配偶者選択ルールでも,シミュレーションプログラム上では,条件式のかたまりとして,わずか数行で記述できる。
まともに扱うことがきわめて難しく,まだ成功した研究はないと思う。Bongaartsの近成要因モデルに男性側の振る舞いを組み合わせることで,シミュレーションは可能かもしれないが,まだ行われていない。
イベントをプログラム上で「確率的に」起こすために,乱数を用いる必要がある。つまり,確率0.3で起こるイベントであれば,区間(0,1)の一様乱数を発生させて,その値が0.3未満だったらイベントを発生させ,そうでなければイベントを発生させない。数学的に厳密に乱数列を扱うことはきわめて難しいが,シミュレーションで用いるには,次の4つの性質を満たしていれば,擬似乱数列で十分である。(1)連続する値の間に線形の関係がない,(2)値の出現する期待値がとりうる値全部について等しい(と期待される),(3)短い周期の繰り返しにならない,(4)算法で生成される(初期値を決めれば再現できる)。実用上は,Lehmerの乗算合同法で生成される,次のアルゴリズムで十分である(Park and Miller, 1988)。
double frand() { unsigned long const a=16807L; unsigned long const m=2147483647L; unsigned long const q=127773L; /* m div a */ unsigned long const r=2836L; /* m mod a */ extern unsigned long rnd_seed; unsigned long lo,hi; signed long test; double xf; hi=rnd_seed/q; /* / equals to div for int operand */ lo=rnd_seed%q; test=a*lo - r*hi; if (test>0) { rnd_seed=test; } else { rnd_seed=test+m; } xf=(double)rnd_seed/(double)m; return xf; }
また,「準数値算法 乱数」の著者Knuthも,WEB上でC言語による乱数生成器を公開しているし,慶応大学の松本さんと西村さんが開発したメルセンヌ・ツイスターという素晴らしい乱数発生器がCやFORTRANのコードとともにBSDライセンスで公開されているので,これを使う手もある(というか,これが現在利用できる最高の擬似乱数列だと思う)。
cflowやlint,GNUで公開されているChecker,あるいは小さなデータでのテスト実行を通してのプログラムの評価とは別に,モデルの評価(妥当性検証と感度分析)を行う必要がある。結果をどの変数で評価するかということと同時に,表現方法が問題となる。人口再生産シミュレーションで評価対象としてよく選ばれる変数は,人口あるいはその増加率である。目的によってはもっと特異的な変数を選んでもよいが,少なくとも妥当性検証については,実測データが得やすいことが選択基準となる。表現方法としては,平均と標準偏差であらわされることが多いが,注意すべきなのは,一例でも「偶然変わったことが起こりうる」なら重要だということである。つまり,たとえば,同じパラメータで乱数列の初期値だけが異なるシミュレーションを100回実行した場合に,100回中99回は人口が単調増加したとしても,残りの1回で集団が絶滅するならば,その1回の試行には大きな意味がある。それゆえ,結果の表現としては,全部のシミュレーション結果のグラフ表示を含むべきである。
妥当性検証については,「1度でも実測値と完全に一致するならば,現実がそのシミュレーションで想定したメカニズムの通りに起こった可能性を否定できなくなる」という意味で,完全一致するまでパラメータや乱数を変えて試行を繰り返すのが王道である。Dyke (1981)もMacCluer (1980)も,これを「シミュレーションならではの妥当性検証」として強く支持している。しかし現実には,平均値の変化のグラフを書いて,実測データの変化と重ねて「だいたい合っている」といって済ましている論文が多くみられるのは,完全一致を達成するのが難しいからである(不可測なパラメータを増やせば容易だが,それでは意味がない)。次善の策としては,平均値に標準偏差や標準誤差をつけたり,中央値に上下5パーセンタイルまでの幅をつけ,その範囲内に実測データが入ることを確認するとか,カイ二乗適合度検定で「適合していないとはいえない」ことで消極的に妥当性検証を行うことである。
感度分析とは,モデルが妥当であるといえた上で,個々のパラメータを変えたときの結果の変化を検討することである。「単純な解析モデルが存在しないような複雑な現象のシミュレーション」モデルを作る目的はここにある場合が多い。分布を求めることが目的なら必須ではないが,パラメータの性質を知るためにも有効なので,常に感度分析はすることを奨める。
参考プログラムfedのソースと英文ドキュメントを,ここからダウンロードすることができる。コンパイルするには,gccなどANSIレベルのCをサポートした処理系が必要で,実行するには広大なメモリ空間とディスクの空きが必要だが,おそらく大抵のUnix系OS上で動くと思う。DOSやWindows95/98上での実行はあまりお薦めしないが,djgppなど(Windows 95/98/2000/NTなどのWin32環境では,CygnusのGNU-Win32プロジェクトまたはMinGWがお薦め)でコンパイルすれば不可能ではない。実際,Windows 2000/NTでMinGWを使えば,unix系OSに拘る必要はなくなったと言っていいと思う。
このプログラムは,パプアニューギニア低地に居住するギデラという人たちの個人レベルの人口再生産をシミュレートしたものである。個人の定義には配偶関係(一夫多妻婚がある)とマラリア感染,それに対する抵抗性の遺伝子を含み,マラリア流行度の異なる2地域間での移動があるものとし,死亡の発生についてはBrassのロジットモデルを用い,ペア形成のルールはある年齢差以内のクラン外婚かつ一夫多妻を許すもので(ただし可能なペアについては等確率),ペアの出生力は扱わず従来通りの女性の出生力にした,甚だプリミティヴなものであるが,単純なので,シミュレーションの理解には役に立つと思う。
Correspondence to: minato@ypu.jp.