読者です 読者をやめる 読者になる 読者になる

Take a Risk:林岳彦の研究メモ

自らの研究に関連するエトセトラについてのメモ的ブログです。主にリスク学と統計学を扱っています。

進化学者のためのMCMCのアナロジカルなアヤシイ解説/その3

進化 統計

今回は本解説シリーズ(第1回第2回)の最終回として、「無性生物の進化シミュレーションとしてのMCMC(マルコフ連鎖モンテカルロ法)」について解説していきます。今回は長いっす。

無性生物集団の進化は「マルコフ連鎖」&「モンテカルロ」

まずは「マルコフ連鎖」「モンテカルロ」という語句の説明から入りたいと思います。

まずはマルコフ連鎖ですが、Wikipediaから引用すると:

マルコフ連鎖は、未来の挙動が現在の値だけで決定され、過去の挙動と無関係である

と書いてあります。つまり、ざっくり言うと「時間t+1における状態が時間tの状態のみに依存する」ということになります。一方、「モンテカルロ」はざっくり言うと何らかの形で「確率的に生成される乱数*1を使う」ことを表します*2


これらの用語を「無性生物の進化」のアナロジーで説明すると以下の図ように表せるかと思います:

それぞれの丸は無性生物の個体を、数列はその遺伝子型を表しています。上の図では5遺伝子座を想定しています。最初の個体の遺伝子型は"00000"ですが、世代ごとに突然変異が起こることにより系統の遺伝子型が変化していきます*3

図を見てわかる通り、「子の遺伝子型(世代t+1の遺伝型)」は祖父母やその遺伝子型ではなく「親の遺伝子型(世代tの遺伝型)」にのみ依存して決まります。これは「マルコフ連鎖」に相当します。また、世代ごとにランダムに生じる突然変異により遺伝子型は親のものから変化します。ここの部分は「モンテカルロ」に相当します。

つまり以上のような無性生物の進化は、「マルコフ連鎖」「モンテカルロ」的な系であると理解できます。

進化シミュレーションとしてMCMCのアルゴリズムを「読む」

それではMCMCのアルゴリズムとして一般的なメトロポリス-ヘイスティング法について書いていきます。メトロポリス-ヘイスティング法についてはWikipediaの記事がわりとまとまっていますので、以下、Wikipediaの記事を引用しながらその解説という形で書いていきます。ちょっととっつきにくいかもしれませんが、あとでまた具体的なRコードをみながら復習していきます。

M-Hアルゴリズムの手順(1):突然変異による新しい遺伝子型の生成

Wikipediaの記事のM-Hアルゴリズムの手順を見ると、その最初には以下のように書いてあります。(以下、Wikipediaの記事の方を参照しながらお読みください)

最新のサンプル値はx^tであるとする。

ここでx^tは世代tにおけるパラメータセットの一揃いの値のサンプルを指します。「パラメータ=座位」「パラメータの値=対立遺伝子」という進化アナロジーで解釈すると、x^tは世代tにおける遺伝子型を表すことになります。

*尚、以下では無性生物「集団」の遺伝子型の進化を想定して解説していきます。なので「突然変異の集団内での固定」という要素が新しく入ってきますのでご注意ください。また以下での「世代」の意味は、多少曖昧になりますが、集団内で突然変異が生じて固定するくらいの長さをもった時間単位としての「集団にとっての世代」という意味で解釈していただければ幸いです。

次の

次に確率密度関数Q(x'; x^t)を用いて新しい状態の候補を提案し

というところは突然変異の計算に相当します。突然変異の大きさとバイアスは関数Q(x'; x^t)によって決まります。突然変異による新しい遺伝子型(x')は現在の遺伝子型(x^{t})のみに依存しながら(=マルコフ連鎖)、突然変異は確率密度関数Q(x'; x^t)に従いランダムに生じます(=モンテカルロ)。
一般に*4、MCMCにおいてはパラメータが複数ある場合にはパラメータを1つづつサンプルしていきますが、これはAdaptive Dynamicsなどの解析で良くみられる「集団内多型として一度に考慮するのはsingle mutationによって生じたone mutant genotypeのみ」という仮定に相当するとも解釈できます*5

また

が標本候補 x^{'}とその一つ前の標本x^{t}の比率であり、
 a_1 = \frac{P(x')}{P(x^t)}

ここは尤度比( a_1)の計算に相当します。ここで用いられている"比率"という語は、正確には"尤度比"ですかね?*6。「尤度=適応度」という進化アナロジーに従うと、ここは「現在集団内に固定している遺伝子型(x^{t})」と「突然変異により生じた新しい遺伝子型(x')」の"適応度の比"を計算していることに相当します。

 a_2 = \frac{Q( x^t; x' )}{Q(x';x^t)}
が提案密度の両方向( から へとその逆方向)の比率となるような値
 a = a_1 a_2
を計算する。提案密度が対称の場合は aは1 となる。

ちなみに最後の行のaは" a_2"の誤植と思われます*7。以降、説明のために突然変異にバイアスがない場合( a_2=1)についてのみ考えていきます。 a_2=1の場合、 aは尤度比=適応度比そのものになります(つまり a = a_1)。

M-Hアルゴリズムの手順(2):適応度の比に応じた集団内での固定の計算

次はサンプルされたパラメータの採択の計算になります。

次に、以下の規則を用いて新しい状態  x^{t+1} を決定する。

 \displaystyle a \geq 1 の場合:

 \displaystyle x^{t+1} = x'

 \displaystyle a < 1 の場合:

[tex
\displaystyle a] の確率で  \displaystyle x^{t+1} = x'
 \displaystyle 1 - a の確率で \displaystyle x^{t+1} = x^t

この部分は、進化シミュレーションとして解釈した場合には、新しい遺伝子型の集団内での固定を計算している部分にあたります。

ここれでの固定のルールは、 \displaystyle a \geq 1 の場合、つまり「突然変異により生じた新しい遺伝子型(x')の適応度のほうが現在の遺伝子型(x^{t})よりも高い」場合には、突然変異遺伝子型が集団内に固定する( \displaystyle x^{t+1} = x')というものです。これは強い正の自然選択による決定論的過程を扱っているものと捉えられます。

一方、 \displaystyle a < 1 の場合、つまり「突然変異により生じた新しい遺伝子型(x')の適応度のほうが現在の遺伝子型(x^{t})よりも低い」場合には、突然変異遺伝子が集団内に固定する確率は \displaystyle aとなります。これは、弱有害遺伝子の遺伝的浮動による固定を扱っているものと捉えられます。当然、適応度比( a)が1に近いほど固定する確率は大きくなります。

この前者( \displaystyle a \geq 1 )の決定論的過程だけでは、進化は適応度地形(=尤度関数or事後分布の全体)における局所的最大値のところにトラップされてしまうため、適応度地形全体(=尤度関数or事後分布)を探索していくためにはたまには"山を下る"からくりも必要となり、後者の確率的過程( \displaystyle a < 1 )が必要となってくるわけです。

また、"山を下る"ときに"適応度比( a)に応じて固定する"という制約を設けているため、あまりにも適応度が低い遺伝子型(のクラスター)の方へは進化しない(MCMCのサンプルパスが伸びない)という性質もあります。このような性質により、前回に見たような”ほとんどが穴のチーズ"のようなすっかすかの適応度空間(=尤度関数or事後分布)であっても適応度(=尤度or事後確率)がゼロとならない部分だけを効率的に探索して計算していくことができるわけです。

進化シミュレーションとしてのMCMCを実装してみる

かなり抽象的な話が続いたので、実際に「正規分布の平均値 \thetaを推定する」という例題について、RでのMCMC計算を実装してみたいと思います。

ここで分散は既知( \sigma =1)とし、データは適当にData = c(6,9,10,12)とします。

このとき、尤度関数は:

 L(Data|\theta) = \frac{1}{\sqr{2\pi}} e^{-\frac{(\theta-6)^2}{2}} \frac{1}{\sqr{2\pi}} e^{-\frac{(\theta-9)^2}{2}} \frac{1}{\sqr{2\pi}} e^{-\frac{(\theta-10)^2}{2}} \frac{1}{\sqr{2\pi}} e^{-\frac{(\theta-12)^2}{2}}

となります。ここで、ちょっと苦しいかもしれませんが、データを「環境要因」と解釈し、 \thetaをその環境要因のもたらす選択圧に係る形質の遺伝子型と解釈すると、 L(Data|\theta) はある「環境要因のセット」に対する遺伝子型 \thetaの「適応度地形」を表す関数であると解釈できます。

ここで、平均値 \thetaの事前分布として「平均10・分散無限大」の正規分布を仮定してみましょう。この場合、事前分布は実質的にフラットな無情報分布となるため、平均値 \thetaの事後分布の形は尤度分布(=適応度地形)と同じ形になります。

 P(\theta|Data) \propto L(Data|\theta)P(\theta)
 P(\theta|Data)  \propto L(Data|\theta)

では、MCMCにより事後分布(=適応度地形)の計算をしてみましょう。

以下がそのRコードになります。正直言いますと、私はいつもWinBUGSしか使っていないので、自分でMCMCのコードを書くのは今回が初めてです。ビクビクしながら書いておりますので、もし間違っていたらご指摘よろしくお願いいたします!>識者がた

Data <- c(6,9,10,12)


#対数尤度(Log likelihood)の計算
#あとでまた指数変換して戻して比をとる(対数尤度のまま比をとっちゃだめ!)
#指数変換して戻すときに比を取るので正規分布の規格化定数(1/root(2)pi)の部分は無視してよい
#でもexpの中の分母の2は無視すると計算がオカシクなるので無視しちゃだめ!
LL <- function(theta){
(1/2)*(-(theta-Data[1])^2 - (theta-Data[2])^2 -(theta-Data[3])^2 -(theta-Data[4])^2)
}


#MCMCのパラメータ行列の初期化(50000世代計算する)
theta.MCMC <- numeric(50000)


#初期値(初期遺伝子型)の設定
theta.MCMC[1] <- 10


#MCMCの計算
for(t in 1:50000){
#更新前の対数尤度(適応度)の計算
LL.current <- LL(theta.MCMC[t])



#標本候補(突然変異遺伝子型)の作成
#提案密度(突然変異生成関数)には正規分布N(0,1)を仮定
theta.dash <- theta.MCMC[t] + rnorm(1,0,1)

#候補値の対数尤度(適応度)の計算
LL.dash <- LL(theta.dash)
#尤度比(適応度比)の計算:尤度を対数からもとに戻してから比をとっています
LLR <-exp(LL.dash)/exp(LL.current)
#尤度比 > 1のとき
if(LLR > 1) {
theta.MCMC[t+1] <- theta.dash
}


#尤度比 < 1のとき
if(LLR>runif(1)){
theta.MCMC[t+1] <- theta.dash
}else{
theta.MCMC[t+1] <- theta.MCMC[t]
}
}


#解析的に求めた平均値thetaの事後分布からのサンプリング
theta.post <- rnorm(50000,mean=mean(Data),sd=sqrt(1/length(Data)))


#MCMCサンプルのdensityの表示(黒)
plot(density(theta.MCMC),xlim=c(6.5,11.5),ylim=c(0,0.85),ann=F)
#解析的に求めた事後分布を重ねあわせる(赤)
par(new=T)
plot(density(theta.post),xlim=c(6.5,11.5),ylim=c(0,0.85),xlab="Theta",ylab="Posterior distiburion of Theta",main="",col="red")

#サンプル/進化パスの表示(5000世代まで)
plot(theta.MCMC,type="l",xlim=c(0,5000),xlab="generation",ylab="Sample/Evolutionary path of Theta", col="red")

このRスクリプトのMCMCのサンプルパスを描くと、こんな感じになります:

ここで、このサンプルパス=進化パスにおいて現れた遺伝子型(パラメータセット)ごとの頻度を拾い集めて最終的に総計してみると、「事後分布=適応度地形」が描かれているということになります。つまり、1次元的な「進化の路」を網羅的に線描していくことにより、最終的に「地形」の姿が顕れてくるというわけです。

上の例では、サンプルを集計した後に顕れる最終的な「地形」はこんな感じになります:

赤線がMCMCによって書いた「事後分布/尤度分布/適応度地形」、青線が解析的に得られる事後分布から描いた「事後分布/尤度分布/適応度地形」です*8。MCMCによって「地形」がちゃんと再現できていることが分かります。適応度/尤度が最大となるのは \theta=9.3くらいの遺伝子型/パラメータの値のときのようです。


というわけで、MCMCによる適応度地形/尤度分布/事後分布の計算について見てきました。以上の単純な例は解析的にも解けるので、あまりMCMCのメリットは実感できないかもしれません。しかし、多次元遺伝子(パラメータ)空間を扱う場合には、MCMCによる「進化的」な計算をしないとなかなかその全貌を描くことは難しくなるわけです(詳しくは前回参照)。

考察:インテリジェント・MCMC・デザイン説

さて、上記の論考を眺めてみますと、逆に言えばそもそも「生物の進化」というものは余りにMCMC的であるように見えます*9

何か、怪しい気がしませんか?

立ち止まってよくよく考えてみましょう。「ある種の知性体」の存在を想定せずに、このようなMCMC的なアルゴリズムがこの世界に埋め込まれていることを果たして説明できるでしょうか?

もしかしたら、われわれの『この世界』自体が、「神のイコンとしての多次元空間内での適応度地形」を描くための多次元積分を計算するために神によって造られた『マルコフ連鎖モンテカルロ計算機の実装』である、のかもしれません。

よく分かりませんが、おそらく、きっとそうに違いありません。
我ながら思いがけないことに、進化とMCMCの考察を通して、非常に斬新な進化理論および神の存在証明にまで辿りつくことができました。何でも考えてみるものですね。

今後、この学説をとりあえず「インテリジェント・MCMC・デザイン説(I-MCMC-D説)」と呼んで普及していきたいと思います。きっとこの説は、空飛ぶスパゲッティ・モンスター説と並んで現代進化論に対する有力なオルタナティブとして今後注目されていくことでしょう*10

最後は真面目に:フィッシャーの掌の上で

さて、賢明な読者の皆様がたは薄々感づかれていたかもしれませんが、本解説シリーズの全体*11は偉大なるR. A. フィッシャーに対するオマージュを意図して書かれております*12

結局のところ、我々は未だにフィッシャーの掌の上で右往左往しているだけなのかもしれない、と私はときどき思うことがあります。


関連文献

R.A. Fisher - The Life of a Scientist

R.A. Fisher - The Life of a Scientist

Fitness Landscapes and the Origin of Species (Monographs in Population Biology)

Fitness Landscapes and the Origin of Species (Monographs in Population Biology)

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

*1:古典的な例でいえば「サイコロを振って出た目」とか

*2:たぶん

*3:説明の便宜のため、本記事では常に「1倍体」かつ「表現型=遺伝子型」であるケースを想定して説明していきます

*4:私の見聞きした範囲では

*5:たぶん

*6:英語版の記事では"Likelihood Ratio"と書いてある/文脈によっては事後確率の比を示すこともあるから単に"比率"って書いてあるのかな?

*7:Wikipediaの訂正ってどうやるのかな?

*8:詳しくは上記のRのスクリプトを参照

*9:伊庭先生の前回の記事へのコメントによると生物の進化は寧ろ逐次モンテカルロ的らしいです。また改めてその辺りはチェックしてみたいと思います。

*10:というかこの考察のくだり、単に「スパゲッティ・モンスター教」を引用してみたかった、というただそれだけのために書いています。ごめんなさい

*11:I-MCMC-D説を除く(もちろん)

*12:フィッシャーの頭のなかでは「Likelihood of survival and reproduction」と「Likelihood of data」の概念はどういう案配で繋がっていたのでしょうかね