ベータ分布は現象例としては少ないが、ベイズ統計学での役割は大きいという（1）。Yang本でも、系統解析のベイズ推定のところで出てきていた。
ベータ分布
ベータ分布（Beta distribution）は、０から１までの上の確率分布で、確率密度関数は以下のように定義される。
ここで、α > 0、β > 0である。ベータ分布はBe(α, β)と表す。
ベータ分布の平均と分散
確率変数Xがベータ分布Be(α, β)に従う場合の期待値や分散は以下のように示される。
規格化する定数として使われているベータ関数
ベータ分布の確率密度関数で、B(α, β)が分母にあるが、これは、積分して1にするための定数となる。定義は以下の通りで、
となる。これをベータ関数という。ベータ関数とガンマ関数との間の関係は、以下の通りになる。
ベータ分布の特徴：αとβの値で自由に形が選べる
ベータ分布のαとβの値をちょっとずつ変えていくと、いろいろな確率密度関数が作れる。「いろいろな 確率的主観を表すのに都合がいい」という表現をされているが、それがベイズ統計学で使われている理由なのだろうか。αとβの値を変えるとどんな風に分布が変化するのか実際にRでプロットしてみた。スクリプトとして簡単に以下のものを書いてプロットしている。
# Plot Beta distribution
betadist <- function(i, x, a, b, c1, c2){
# First plot
if (i==1) {
plot(x, dbeta(x, a, b), type = "l", col = i, ylim = c(0,5),
xlab = "x", ylab = "")
# i>2 plot
} else {
par(new=T)
a <- a + c1*(i-1)
b <- b + c2*(i-1)
plot(x, dbeta(x, a, b), type = "l", col = i, ylim = c(0,5),
xlab = "", ylab = "", axes=F)
}
}
x <- seq(0.01, 1.0, len = 500)
# Example 1
for (i in 1:5) {
betadist(i, x, a=1, b=1, c1=-0.1, c2=1)
}
最初のプロットはα=β=1で、一様分布になる。下のプロットの黒い線が該当する。αは0.1ずつ減らしていき、βは1つずつ増やしていく感じで変化していくと、どんどん急勾配のプロットになっていく。
αとβを一つずつ増やしていくと以下のように変化する。
for (i in 1:20) {
betadist(i, x, a=1, b=1, c1=1, c2=1)
}
右側が大きくなるように形を変えることもできる。下は、αを増やしながら、βを減らした場合。
for (i in 1:10) {
betadist(i, x, a=1, b=1, c1=1, c2=-0.1)
}
α=0.3、β=0.3に指定すると、鍋のような形を作ることもできる。
このようにαとβの値によって形がいろいろと変化するので、確率分布の形に関数を当てはめたい時に用いられる。
先に述べたように、α=β=1の時は、一様分布となる。 この時の期待値E(X)=1、分散V(X)=1/12となる。
Site modelでomegaを表すために使う
random-sites modelでは、M7やM8のモデルにおいて、dN/dSの分布についてベータ分布を適用している。M7はベータ分布を使い、M8では、ベータ分布＋Ws (positive selectionを受けているサイト)となる。M8では、０から１までの値のOmegaに加えて、正の自然淘汰を受けているサイトを追加のクラスとして足していることになる。上で示したように、ベータ分布はαとβの指定により形が柔軟に変わる。ベータ分布を用いることで、帰無仮説のOmegaの分布が柔軟に設定できることが利点となる。codemlを使って推定すると、このαとβに当たる値が推定される。例えば、Yang本のp277のTable 8.3には192のMHC allelesについてのパラメータの推定値が書かれているM7とM8のモデルではベータ分布のパラメータの推定値が出ている。M8モデルでは、ベータ分布に従うサイトと、positive selectionを受けているサイトの割合(P0, P1)も推定されている。
| Model | #p | l | Estimates of parameters | Positively selected sites |
| M7 | 2 | -7498.97 | p = 0.103, q = 0.354 | Not allowed |
| (beta) | ||||
| M8 | 4 | -7238.01 | p0 = 0.915 (p1 = 0.085) | 9F, 24A, 45M, 63E, 67V, 69A, |
| (beta&W) | p = 0.167, q = 0.717, | 70H, 71S, 77D, 80T, 81L, | ||
| Ws = 5.079 | 82R, 94T, 95V, 97R, 99Y, | |||
| 113Y, 114H, 116Y, 151H, | ||||
| 152V, 156L, 163T, 167W |
#pがパラメータ数、lはlog likelihoodの値、pとqがベータ分布のパラメータに該当する。正の自然淘汰の影響を受けていると思われるサイトが右の列に示されている。M8モデルでは、Omegaの値がWs=5.079と推定されるpositively selected sitesの割合がp1=0.085と推定されている。
パラメータの推定値をもとに、omegaの分布を(0,1)の範囲でプロットしてみると、以下のようになる。M8の場合には、実際にはこの右っ側に1以上の値を持つサイトがあると推定されていることになる。
M7で推定されたパラメータによる分布を黒、M8のを赤で示している。M8の場合、これに加えて、positive selectionのサイトが0.085の割合で存在することになる。
参考文献
1. 東京大学教養学部統計学教室 編 統計学入門 I pp. 126-1262. Ziheng Yang (2006) Computational Molecular Evolution. Oxford University Press.