JC69の定常分布
JC69で求めた遷移確率は、\(t\)→∞の極限操作で全て\(\frac{1}{4}\)に収束することが述べられています。極限操作をmaximaで試してみましょう。
1. 極限操作
関数limitを用いてこの処理を行うことができます。前回のJC69で作成した遷移確率行列\(P\)の要素を指定して、極限を撮ってみます。\(P\)の対角要素はすべて同じ形式であり、同様に非対角要素も全て同じ形式となっています。まず、対角要素について、\(P\)の１行1列の要素を用いて計算してみます。
limit(P[1][1],x,inf);
テキストでも説明されているように、\(t\)とλを分離できないので、時間\(t\)ではなく距離\(x(=λt)\)に対して極限操作を行っています。 \(\frac{1}{4}\)が出力されます。同様に、非対角要素として１行2列を用いて計算します。
limit(P[1][2],x,inf);
これも\(\frac{1}{4}\)が出力されます。
遷移確率行列\(P\)そのものの極限をとることもできます。
limit(P,x,inf);
これを実行すると、次のように遷移確率行列はその要素が全て\(\frac{1}{4}\)になることがわかります。
\begin{equation} \left[ \begin{array} {rrrr} \frac{1}{4} &\frac{1}{4} &\frac{1}{4} &\frac{1}{4} \\ \frac{1}{4} &\frac{1}{4} &\frac{1}{4} &\frac{1}{4} \\ \frac{1}{4} &\frac{1}{4} &\frac{1}{4} &\frac{1}{4} \\ \frac{1}{4} &\frac{1}{4} &\frac{1}{4} &\frac{1}{4} \end{array} \right] \end{equation}ここで、任意の塩基組成\((a,b,c,d)\)を考えます。塩基組成ですから、\(a + b + c + d = 1\)です。
v:[a,b,c,d];
\(P\)と\(v\)の積をとると、初期状態の塩基組成\([a,b,c,d]\)から出発して, \(t\)時間後の塩基組成が得られます。 \(t\)→∞での塩基組成を極限操作でもとめてみましょう。
ratsimp(limit(P.v, x, inf));
\((a + b + c + d)\)は、初期状態の塩基頻度の和ですから１になります。従って、初期状態がどのような塩基組成から 出発しても、\(t\)→∞では塩基組成は\((1/4, 1/4,1/4,1/4)\)になります。この極限分布を定常状態分布あるいは定常分布とよんでいます。 遷移確率の(1/4)への収束の様子が図1.3に示されています。次にこれを作図してみましょう。
2. 遷移確率の収束の様子
maximaのplot2dを使って作図します。
plot2d([P[1][1],P[1][2]],[x,0,2]);
plot2dの最初の[]内にプロットする関数が記述されています。次の[]内にプロットする\(x\)の範囲が記述されています。図をファイルに保存する場合は、次のようにオプションを加えてください。
plot2d([P[1][1],P[1][2]],[x,0,2],[gnuplot_term, "postscript eps color"],[gnuplot_out_file,"fig1.3a.eps"]);
ここで、gnuplot_termオプションをcolorあるいはmonoにすることで、カラーかモノクロかを指定できます。また、gnuplot_out_fileオプションで出力ファイルを指定できます。
P[1][1]が\(p_{TT}(t)\)に、またP[1][2]が\(p_{TC}(t)\)に相当します。これを実行すると、次のような図が得られます。大方の傾向は似ていますが、図1.3では距離=2のところで、0.25にこの図に示すほど近接していません。
極限操作で見たように、P[1][1]やP[1][2]は、距離\(x\)の関数になっており、ここまでの処理では\(x=λt\)としています。しかし、図1.3では距離は\(3λt\)とおかれています。\(3λt\)がJC69の進化距離に相当するからです（次の「JC69の距離計算」で説明します）。そこで前回から使ってきた遷移確率行列Pの代わりに、新しく\(p_{TT}(t)\)と\(p_{TC}(t)\)に相当する関数で、\(x\)を\(x/3\)に変更したものを定義し、それを使って作図してみましょう。maximaでは、関数名:=式 の形で関数を定義できます。
ptt(x):=0.75*exp(-4*x/3) + 0.25;
ptc(x):=0.25 - 0.25*exp(-4*x/3);
plot2d([ptt(x),ptc(x),0.25],[x,0,2]);
すると,次の図が作成されて、図1.3が再現されたことがわかります。先ほどの図では入れ忘れていた収束値の0.25のプロットも含めました。
前ページ JC69 | 次ページ JC69の距離計算