群馬大学 | 医学部 | サイトトップ | Software Tips
最終更新： 2004年 4月 15日 (木曜日) 07時42分 （MTの図も追加）
このページは，GPLに従って配布されているRを使ってシミュレーションをした例をメモしておくために作成した。Rについての詳細は，統計処理ソフトウェアRについてのTipsを参照されたい。
Rでは，1.7.0からデフォルトのRNGがメルセンヌツイスターとなった。制御構造も簡単に書ける上，グラフ作成までできるので，構成要素が単純な配列で示せるようなシミュレーションに使うにはCよりも向いていると思う。
なお，RNGはシミュレーションで確率的なイヴェントを表すために使うわけだが，メルセンヌ・ツイスターの利点は，周期がほぼ10の6000乗くらいと長く，多次元疎結晶構造もほとんど問題にならないところにある。古典的には乗算合同法を含む線型合同法がよく用いられてきたが，周期の短さと多次元疎結晶構造が線型合同法の欠点であることは既にいろいろなところで指摘されている。Rにも入っているMarsaglia-Multicarryが生成する擬似乱数列を4個ずつ組にして４次元空間に割り付け，その１つの軸の値をきわめて狭い範囲に限った，いわば３次元スライスを描画してみると，下図のように直線状の点の並びができてしまうのが一目瞭然であろう。
ちなみに，上図を描くプログラムは下記のとおり。乱数をいっぺんに発生させようとするとWindows環境ではメモリが足りなくなるので，いったん短い配列で３次元散布図を書かせておいて，重ね書きしていく。
# ライブラリ呼び出し
library(scatterplot3d)
# 重ね書き用の関数定義
addpoints <- function(s3d) {
U<-runif(400000,0,1)
x<-1:100000; y<-x; z<-x; q<-x
for (i in 1:100000) {
x[i]<-U[i*4-3]; y[i]<-U[i*4-2]; z[i]<-U[i*4-1]; q[i]<-U[i*4]
}
s3d$points3d(x[q>0.9999],y[q>0.9999],z[q>0.9999],pch=".")
}
# 擬似乱数列の指定
RNGkind("Marsaglia-Multicarry")
set.seed(1)
# １回目の描画
U<-runif(400000,0,1)
x<-1:100000; y<-x; z<-x; q<-x
for (i in 1:100000) {
x[i]<-U[i*4-3]; y[i]<-U[i*4-2]; z[i]<-U[i*4-1]; q[i]<-U[i*4]
}
# 出力用デバイスを開く
png("mmbad.png",width=480,height=480)
# ３次元プロット
graph<-scatterplot3d(x[q>0.9999],y[q>0.9999],z[q>0.9999],grid=F,
pch=".",angle=50,scale.y=0.8,main="Marsaglia-Multicarry")
# 続く2500回分を重ねて描画
for (j in 1:2500) { addpoints(graph) }
# デバイスを閉じる
dev.off()
ちなみに，上のプログラムで"Marsaglia-Multicarry"となっているところ（２ヵ所）を"Mersenne-Twister"と変えて実行すると下図が得られる。ほぼ均一にプロットされているように見える。なお，これだけ長いループを回すとかなり時間がかかる。Athlon XP 2200+，nForce2，512 MBメモリ，Windows2000SP4という環境でも数時間かかったので，実行するなら気分的に余裕があるときにした方がよいだろう。
2003年4月13日と27日に統一地方選があり，投票してきたときに感じた引っかかりが，このシミュレーションを思いついたきっかけである。現在の投票は，複数の候補者から複数の当選者を選ぶのに１票しか投票できないシステムになっている。これは，一人一人の有権者の全権を誰か一人の議員に代表して貰うという考え方に基づいていると思われるが，当選しない人に投票した人の意思はまったく反映されないことになってしまうから，明らかに不合理である。現実には個々の有権者は誰か１人に全幅の信頼をおいて代議を任せるわけではなく，この政策についてはA氏がいいし，この政策についてはB氏を支持したいということもあるだろうから，この人は何％支持，この人は何％支持という具合に振り分けた投票ができた方が，個々の有権者の意思を正しく反映できるのではないかと考えられる。
そこで，個々の有権者が複数の候補者に対して少しずつ支持をshareするという内部構造を考え，どういう投票システムをとったときにそれが正しく反映されるのかを調べてみることにした。
まず15人の候補者がいて，1000人の有権者が，それぞれの候補者を何％かずつ支持しているという状況を想定した。内部構造として平均何％の支持を得ているという条件を割り当てるわけである。しかし，支持しているなら必ず投票するかといえばそうではない。投票数の上限を超えては投票できないし，棄権したり白票を投じる人もいる。そこで，有権者１人ずつについて候補者１人ごとに区間(0,1)の一様乱数を発生させ，それが仮想支持率より低ければ上限を超えないうちは投票するというプログラムを書いた（複数当選の場合は仮想支持率を当選人数倍したものと比較した）。候補者をチェックする順序が結果に影響しないように，有権者１人ずつ，候補者をチェックする順番はランダムに並べ替えた（runifで一様乱数を発生させてrankをとったものを配列の添え字として使うと簡単である）。
投票可能数が１票，３票，15票のそれぞれの場合について検討し，仮想支持率と比べて，どの場合がもっとも良く合致するかを検討した。用いたプログラムはelection.Rである（forやwhileを多用していてRらしくないプログラムだが，たぶん環境に依存せず実行できると思う。ただ，実行するとpng形式のグラフをファイルとしてカレントディレクトリに書き出すので，書き込み不可能なメディアがカレントになっているとエラーを起こすと思う）。本当は何度も実行させてばらつきをみるなどすべきなのだが，複数当選時に１人にしか投票できないと逆転が起こりうることさえ示せればいいので，１回の実行結果しか示さない。
上の図には１つの仮想支持率を与えたときに，３通りの投票システムがどのような得票数分布を生み出すかを示した。左上のパネルは仮想支持率の分布である。15本の棒は，それぞれの候補者に割り当てた支持率を示している。右上のパネルは，15票まで投票していいという条件での投票結果である。有権者各個人が，投票したいと思った人全員に投票できるとすれば，ベースラインよりも少しでも高い仮想支持率が割り当てられた候補者（左から２番目，８番目，10番目，11番目）が多くの得票を集める結果となるのは予想できる。
一方，左下のパネルのように，各有権者が最初に投票を決意した１人にしか投票できないとして，その決意のしやすさが仮想支持率に依存している場合は，他よりも抜きん出て支持率の高い２番目と10番目の候補者は仮想支持率に見合った多くの得票を集めるけれども，それほどベースラインと差がない８番目と11番目の候補者の得票は13番目の候補者の得票より低くなってしまった。これは，８番目や11番目にも投票する可能性があった有権者が２番目や10番目に入れてしまい，13番目の候補者に投票した有権者は，たまたま仮想支持率が高い候補者には入れる気がなかったという状態を示していると解釈できるだろう。現実の選挙でもこれはよくあることだと思う。不合理ではないか？
右下のパネルのように複数人への投票を許すと，状況は改善される。これは３人までの投票を許した場合の図だが，11番目の候補者（つまり３番目に仮想支持率が高い候補者）は，明らかにベースラインよりも多くを得票している。当選が支持率を反映するためには，当選者数分の投票をする方が良さそうだといえると思う。投票システムの改善を求めたいが，いかがだろうか。
なお，別の内部構造を仮定したプログラムelection_alt.Rも作ってみたが，たぶんここに示したものの方が現実に近いと思う。
平成13年度に高崎経済大学経済学部で非常勤講師として生態学の講義をしたとき，後半は，松田裕之「環境生態学序説」（共立出版）を使って具体例を取り上げた。第14回に扱った「魚種交替現象」の中で，松田さんが提案している三すくみモデルは，Rでシミュレーションをする対象として向いていると思うのでここに示す。
「魚種交替」は，簡単に言えば，日本近海の漁業で獲れるマイワシ，マサバ，カタクチイワシという３種類の魚の漁獲高は経時的に上下動し，1980年代には豊漁だったマイワシが，現在では少ししか獲れなくなり，代わってカタクチイワシが増えている現象のように，主に獲れる魚の種類が替わる現象をさす。この現象を説明するモデルは，種間競争と被食捕食の組み合わせ，種内競争，自己振動モデル，プランクトン食害説など多数提案されているが，中でも松田の三すくみ説はモデルとしてエレガントだと思う。
簡単に三すくみ説を紹介すると，３種間の種間競争としてロトカ＝ヴォルテラの競争方程式を拡張したものである。常に微小量の聖域（他種の影響を受けない繁殖場所）からの供給（c）があるとして，
という３本の差分方程式によって３種の魚種の現存量N1, N2, N3が表現できると考えるのである。Rのプログラムとしては，3sukumi.Rとなる。
上の図で３本の線は３つの魚種の現存量の変化を示している。実行させるたびに乱数の違いでパタンは変わるが，必ずピークをずらして振動するグラフになっていることがわかる。
このモデルについてのRのプログラムを群馬大学の青木先生が改良してくださり，掲載許可をいただいたので3sukumi2.Rとして公開する。ぼくのプログラムよりずっとエレガントだと思う。forループの中の式には感動した。
既にNakazawa, M., H. Ohmae, A. Ishii and J. Leafasia (1998) Malaria infection and human behavioral factors: A stochastic model analysis for direct observation data in the Solomon Islands, American Journal of Human Biology, 10: 781-789.として発表済みの論文（やや詳しい概要）で使ったマラリア感染と行動防御による介入のモデル（Cで書いたもの）をRに移植すること。Rの方が乱数生成アルゴリズムがいいことと，グラフ出力が簡単なことから，この場合のようにデータ構造が単純な場合はCで書くよりも向いているように思う。
モデルの詳細は論文を参照されたい。基本構造はヒト個体群についてSEIR，蚊の個体群についてSEIで，ヒトが日々確率的に行動防御をするというものである。Rで予めいくつかの条件を設定ファイルに書いておいて，条件ごとに別々のファイル（個々のファイル名を設定ファイルに書いておく）に結果を出力するように設定したプログラム[malariasim.R]と設定ファイル例[malsimauto.dat]をダウンロードできるようにしたので試されたい。
なお，このファイルでは，シミュレーション実行期間をmaxdateという変数に２年間（日単位なので365*2）と設定しており，シミュレーション回数を50回（メインルーチンにforループを2:50としている部分）と設定している。また，軸の名前に日本語を表示するために，なかまさんの日本語化パッチを使っている。日本語化していない環境では，軸の名前を設定している部分（plot関数のxlabとylab引数）を英数字に書き換えないと文字化けする。
上記ファイルを実行した結果できるグラフの例（乱数は実行するたびに違うので）を下に示す。
| 防御実行率 | 防御効率30％ | 防御効率80％ |
|---|---|---|
| 平均12％ |  |  |
| 平均70％ |  |  |
| 平均95％ |  |  |
Correspondence to: minato@phi.med.gunma-u.ac.jp.