読み書きプログラミング

日常のプログラミングで気づいたことを綴っています

(44) 種の絶滅

遺伝子ドリフトと類似の問題である種の絶滅について考えます。
2種類の種A, Bがあり、もし全個体数が一定nで、次の世代の占める割合の確率がそれぞれの個体数比に比例するとします。
Aの個体数がiの時、次の世代の個体数それぞれの確率は、

Tは遷移確率行列と呼ばれ、二項分布成分を持ちます。後で使うので、i=1...n-1, j=1...n-1の範囲でTと同じ値を持つn-1行n-1列行列をQとしておきます。

t世代後の確率分布は、

ここでp(t)はインデックス0からnまでのn+1の成分を持つ確率分布のベクトルです。後で使うので、i=1...n-1の範囲でp(t)と同じ値を持つn-1列ベクトルをqとしておきます。

絶滅状態はi=0, またはi=nです。従って、初期状態がiの状態から絶滅までの平均世代数は、以下のように書けます。

のMaclaurin展開ですから、結局以下の式を得ました。

さて、n=5として、これをMaximaで計算してみましょう。

/* 7.3.19m */
n : 5$
T[i, j] := binomial(n, j)*(i/n)^j*(1 - i/n)^(n - j)$
Q : genmatrix(T, n - 1, n - 1)$
print(Q)$
M : invert(ident(n - 1) - Q)$
print(M)$
list : apply("+", args(transpose(M)))$

使用した機能


各個体比の場合の、絶滅までの世代数の厳密な期待値が計算できました。