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

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

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

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

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

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

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

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

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

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

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

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

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

M-Hアルゴリズムの手順（１）：突然変異による新しい遺伝子型の生成

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においてはパラメータが複数ある場合にはパラメータを１つづつサンプルしていきますが、これは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アルゴリズムの手順（２）：適応度の比に応じた集団内での固定の計算

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

次に、以下の規則を用いて新しい状態 $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による「進化的」な計算をしないとなかなかその全貌を描くことは難しくなるわけです（詳しくは前回参照）。