2013年8月13日火曜日
混合分布のパラメータ推定：その2（ガンマ混合モデル）
（このページは URL を2013ベントス・プランクトン学会大会の要旨に掲載したものからリンクしています）
# 4部構成にしてあります、こちらはガンマ混合モデルを使用の例。
# その1：正規混合モデル（オーソドックスな基本形）
# その2：ガンマ混合モデル（体長の増加に伴い分散も増加するモデル）
# その3：前後関係からユーザがグループ数を仮定して正規混合モデルを実行
# その4：体長と体幅など複数種類のデータで分ける（有効な最終手段）
*******************************************************
（その1からの続きです）
体サイズ頻度分布のヒストグラムでは、複数のコホートがオーバラップすることは珍しくなく、実用的には多数の重なったグループの分離はできてほしい。
まずは、その1で使用したプログラムを試したところ、平均の推定はできた。しかし、標準偏差の推定はうまくいかず、とくに真ん中のグループが狭められてしまう傾向にあった（最後のところに推定の図を載せておく）。
（注：推定の粗さが気にならなければ、"その1"の推定だけで十分だと思います）
（注：推定の粗さが気にならなければ、"その1"の推定だけで十分だと思います）
解決策としては、目的に応じてモデルに制約を掛けることが挙げられる。ひとつの方法は分散がグループ間で共通という仮定することである。これならばその1で使用したプログラムを微修正するだけで十分に推定が可能だった（mclust()にmodelNames="E"の項を加える）。
今回はコホートの分離が目的なので、体サイズの成長に伴いコホートの分散は増加するという仮定を加えてもよいだろう。それならばガンマ分布を使用すれば、平均の増加に伴いバラツキが増加するのでグループ毎に変更するパラメータを1つ減らすことができる。
その1と同様に、ガンマ分布の母集団を作り、そこから抽出したサンプルデータから元の母集団のパラメータを推定できるか確認し、手法の有効性を確かめた。
ガンマ分布ではrate、shapeといったパラメータを使用するが、平均や標準偏差とは以下のような対応関係にある。
# 母集団の設定
N <- 200 # サンプル全体のN数
mean1 <- 10
mean2 <- 15
mean3 <- 20
sd <- 1.2 # 標準偏差は共通
ratio <- c(0.5, 0.3, 0.2) # 5 : 3 : 2 の比率で混ぜる
N1 <- N*ratio[1]
N2 <- N*ratio[2]
N3 <- N*ratio[3]
# rate, shapeパラメータへ変換
rate1 <- mean1 / sd^2
rate2 <- mean2 / sd^2
rate3 <- mean3 / sd^2
shape1 <- mean1 * rate1
shape2 <- mean2 * rate2
shape3 <- mean3 * rate3
# shapeは平均 * rateとなっており、このため標準偏差をグループ間で共通にしても平均に比例して値が増加する。
set.seed(1) # 乱数のタネを指定
LengthG <- c(rgamma(N1, rate=rate1, shape=shape1), rgamma(N2, rate=rate2, shape=shape2),
rgamma(N3, rate=rate3, shape=shape3)) # ガンマ分布3グループを混ぜる
嫌な感じに混ざった…標準偏差は共通にしたにもかかわらず、山は次第に幅が増加するのがガンマ分布の性質
その1同様に、Mclustの推定結果を初期値に用いた。じつはガンマ混合分布も扱えるより汎用性の高いパッケージもあるのだが推定があまりよくなかったため使用しなかった（ページの末を参照）。Mclustは正規混合分布なので分布型が違うのだが、あくまで初期値としての使用なので構わないだろう。
library(mclust)
plot(mclustBIC(LengthG))
# グループ分け数のモデルをBICで比較
その1同様に、Mclustの推定結果を初期値に用いた。じつはガンマ混合分布も扱えるより汎用性の高いパッケージもあるのだが推定があまりよくなかったため使用しなかった（ページの末を参照）。Mclustは正規混合分布なので分布型が違うのだが、あくまで初期値としての使用なので構わないだろう。
library(mclust)
plot(mclustBIC(LengthG))
3グループ、分散共通のモデルがベストであるという解析結果が得られた。
mc <- summary(densityMclust(LengthG), parameters=T) # 推定されたパラメータの一覧
もし、前の時点の解析結果などから、今回のグループ数推定がおかしいと判断される場合は、グループ数をユーザ側で指定してもよいだろう。G=3のようにグループ数を指定できる。
# mc <- summary(densityMclust(LengthG, G=3), parameters=T)
Mclustの結果を初期値に用いて、ベイズのガンマ混合モデル推定をした。
mc <- summary(densityMclust(LengthG), parameters=T) # 推定されたパラメータの一覧
もし、前の時点の解析結果などから、今回のグループ数推定がおかしいと判断される場合は、グループ数をユーザ側で指定してもよいだろう。G=3のようにグループ数を指定できる。
# mc <- summary(densityMclust(LengthG, G=3), parameters=T)
Mclustの結果を初期値に用いて、ベイズのガンマ混合モデル推定をした。
n.g <- 3 # グループ数の設定
data <- list(n.g=n.g, N=length(LengthG), size=LengthG, # 上記のデータ LengthG を使用
alpha=mc$pro, # mclustの推定した混合比率を初期値に使用
M=mc$mean) # mclustの推定した平均値を初期値に使用
par <- list(ratio=mc$pro, # mclustの推定した混合比率を初期値に使用
tau=1/mean(mc$variance),F, # mclustの推定した分散の平均を初期値に使用
m.group=mc$mean, # mclustの推定した平均値を使用
group=sample(1:n.g, length(LengthG), replace=T),F, # 初期値不明のためランダム
sigma=NA, delta.m=NA, shape=NA, rate=NA)
data <- list(n.g=n.g, N=length(LengthG), size=LengthG, # 上記のデータ LengthG を使用
alpha=mc$pro, # mclustの推定した混合比率を初期値に使用
M=mc$mean) # mclustの推定した平均値を初期値に使用
par <- list(ratio=mc$pro, # mclustの推定した混合比率を初期値に使用
tau=1/mean(mc$variance),F, # mclustの推定した分散の平均を初期値に使用
m.group=mc$mean, # mclustの推定した平均値を使用
group=sample(1:n.g, length(LengthG), replace=T),F, # 初期値不明のためランダム
sigma=NA, delta.m=NA, shape=NA, rate=NA)
model <- function() {
for (i in 1 : N) { # iループ
size[i] ~ dgamma(shape[group[i]], rate[group[i]]) # 混合したサイズ分布（グループ毎に変化）
group[i] ~ dcat(ratio[1:n.g]) # グループ選択の事前分布、カテゴリカル分布を使用
} # iループ終わり
ratio[1:n.g] ~ ddirch(alpha[]) # 混合比率の事前分布：ディリクレ分布=ratioの合計が1になる
tau ~ dgamma(1, 0.01) # 分散の逆数、グループ間で共通
for (g in 1:n.g) { # gループ（n of group）
rate[g] <- m.group[g]*tau # rateパラメータ。tauは分散の逆数であることに注意
shape[g] <- m.group[g]*rate[g] # shapeパラメータ（平均 * rate）
m.group[g] ~ dnorm(M[g], 0.5) # 各グループの平均値の事前分布
} # gループ終わり
for (gg in 2:n.g) { delta.m[gg] <- m.group[gg] - m.group[gg-1] } # ggループ
sigma <- sqrt(1/tau) # 標準偏差の事後分布の計算
} # model終わり
source("WBUGS.R") # cf. WBUGS、R2WinBUGSの補助ツール
seg <- wbugs(data, par, model, n.iter=110000, n.burn=10000, n.thin=100) # 計算実行
# 実行結果
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
# ratio[1] 0.515 0.037 0.444 0.489 0.514 0.540 0.590 1.001 3000
# ratio[2] 0.289 0.036 0.221 0.264 0.289 0.313 0.361 1.002 1100
# ratio[3] 0.196 0.029 0.141 0.176 0.196 0.215 0.256 1.001 2600
# m.group[1] 10.141 0.145 9.853 10.040 10.140 10.240 10.430 1.001 3000
# m.group[2] 15.021 0.232 14.590 14.870 15.020 15.170 15.500 1.001 3000
# m.group[3] 19.984 0.242 19.510 19.820 19.980 20.150 20.470 1.003 970
# sigma 1.326 0.085 1.179 1.266 1.319 1.378 1.510 1.001 3000
# delta.m[2] 4.881 0.245 4.402 4.719 4.875 5.045 5.353 1.001 3000
# delta.m[3] 4.963 0.299 4.364 4.761 4.969 5.165 5.537 1.001 3000
# shape[1] 59.152 7.140 45.460 54.287 59.025 63.972 73.270 1.001 3000
# shape[2] 129.826 15.976 98.769 119.000 129.650 140.500 161.500 1.001 3000
# shape[3] 229.892 28.998 174.195 210.200 229.550 249.500 288.800 1.001 3000
# rate[1] 5.835 0.714 4.493 5.347 5.825 6.320 7.244 1.001 3000
# rate[2] 8.644 1.066 6.618 7.918 8.631 9.370 10.750 1.001 3000
# rate[3] 11.503 1.440 8.727 10.540 11.480 12.480 14.420 1.001 3000
# deviance 680.517 15.944 655.697 669.200 678.400 689.400 717.600 1.002 1500
source("WBUGS.R") # cf. WBUGS、R2WinBUGSの補助ツール
seg <- wbugs(data, par, model, n.iter=110000, n.burn=10000, n.thin=100) # 計算実行
# 実行結果
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
# ratio[1] 0.515 0.037 0.444 0.489 0.514 0.540 0.590 1.001 3000
# ratio[2] 0.289 0.036 0.221 0.264 0.289 0.313 0.361 1.002 1100
# ratio[3] 0.196 0.029 0.141 0.176 0.196 0.215 0.256 1.001 2600
# m.group[1] 10.141 0.145 9.853 10.040 10.140 10.240 10.430 1.001 3000
# m.group[2] 15.021 0.232 14.590 14.870 15.020 15.170 15.500 1.001 3000
# m.group[3] 19.984 0.242 19.510 19.820 19.980 20.150 20.470 1.003 970
# sigma 1.326 0.085 1.179 1.266 1.319 1.378 1.510 1.001 3000
# delta.m[2] 4.881 0.245 4.402 4.719 4.875 5.045 5.353 1.001 3000
# delta.m[3] 4.963 0.299 4.364 4.761 4.969 5.165 5.537 1.001 3000
# shape[1] 59.152 7.140 45.460 54.287 59.025 63.972 73.270 1.001 3000
# shape[2] 129.826 15.976 98.769 119.000 129.650 140.500 161.500 1.001 3000
# shape[3] 229.892 28.998 174.195 210.200 229.550 249.500 288.800 1.001 3000
# rate[1] 5.835 0.714 4.493 5.347 5.825 6.320 7.244 1.001 3000
# rate[2] 8.644 1.066 6.618 7.918 8.631 9.370 10.750 1.001 3000
# rate[3] 11.503 1.440 8.727 10.540 11.480 12.480 14.420 1.001 3000
# deviance 680.517 15.944 655.697 669.200 678.400 689.400 717.600 1.002 1500
母集団のratioの値は0.5, 0.3, 0.2に対し、推定値（中央値）は0.514, 0.289, 0.196。平均は10, 15, 20に対し、10.140, 15.020, 19.980。標準偏差は1.2に対し、1.319。標準偏差がやや大きく推定されたものの、他の推定はかなり良い結果だった。
なお、いずれかの R.hatの値が 1 から大きく離れていたり、n.eff 値が1~2桁に落ちている場合は推定が良くない。n.iter, n.burn, n.thinの値を比例して 5~10倍に増やして再試行したほうがよいだろう。
なお、いずれかの R.hatの値が 1 から大きく離れていたり、n.eff 値が1~2桁に落ちている場合は推定が良くない。n.iter, n.burn, n.thinの値を比例して 5~10倍に増やして再試行したほうがよいだろう。
母集団（実線）と推定値（点線）の確率密度を図示した。標準偏差がやや大きく推定されたため山が低めになったがかなり高い推定精度だ。平均のみをグループ間で変化させたにもかかわらず、ちゃんとバラツキが平均に伴い増加するパターンをモデル化することができた。
# ちなみに、上記のグラフを書くコードはこちら
hist(LengthG, breaks=seq(5*floor(0.2*min(LengthG)), 5*ceiling(0.2*max(LengthG)), by=1.25),
ylim=c(0,40), xlab="Length (mm)", main="", lwd=1.5, col=gray(0.9))
# 棒の区画幅をものぐさするためにコードが複雑になっている…
x <- seq(5,25, 0.1) # x軸を決める
lines(x, N1*dgamma(x, rate=rate1, shape=shape1), lwd=1.5) # 実線のカーブ
lines(x, N2*dgamma(x, rate=rate2, shape=shape2), lwd=1.5)
lines(x, N3*dgamma(x, rate=rate3, shape=shape3), lwd=1.5)
n1 <- round(seg$median$ratio[1]*length(LengthG))
n2 <- round(seg$median$ratio[2]*length(LengthG))
n3 <- round(seg$median$ratio[3]*length(LengthG))
# 点線のカーブ
lines(x, n1*dgamma(x, shape=seg$median$shape[1], rate=seg$median$rate[1]), lty=2, lwd=1.5)
lines(x, n2*dgamma(x, shape=seg$median$shape[2], rate=seg$median$rate[2]), lty=2, lwd=1.5)
lines(x, n3*dgamma(x, shape=seg$median$shape[3], rate=seg$median$rate[3]), lty=2, lwd=1.5)
最後に、ちょっとトリッキーだが、3グループの正規分布からのサンプルを元に正規混合モデル（標準偏差はグループ毎に変化）で推定した結果と、同じサンプルをガンマ混合モデル（標準偏差は一定）で推定した結果を比較してみた。
平均：10, 15, 20、標準偏差：1.0, 1.2, 1.4の3つの正規分布の母集団から、計200サンプルを混合比率：0.5, 0.3, 0.2でサンプリングし、サンプルから母集団のパラメータを推定した。
### 正規混合モデルの推定結果（抜粋）
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
# ratio[1] 0.496 0.037 0.424 0.471 0.495 0.521 0.567 1.001 3000
# ratio[2] 0.272 0.040 0.194 0.246 0.272 0.298 0.349 1.002 1900
# ratio[3] 0.232 0.037 0.167 0.208 0.229 0.254 0.309 1.001 3000
# m.group[1] 10.096 0.093 9.919 10.030 10.090 10.160 10.290 1.001 3000
# m.group[2] 14.545 0.168 14.230 14.430 14.540 14.650 14.890 1.001 3000
# m.group[3] 19.748 0.413 18.840 19.540 19.790 20.010 20.400 1.002 2000
# sigma[1] 0.886 0.072 0.761 0.835 0.882 0.930 1.042 1.001 3000
# sigma[2] 0.950 0.178 0.652 0.829 0.935 1.054 1.341 1.000 3000
# sigma[3] 1.760 0.330 1.246 1.537 1.717 1.934 2.522 1.001 2300
### ガンマ混合モデルの推定結果（抜粋）
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
# ratio[1] 0.506 0.037 0.433 0.481 0.505 0.530 0.579 1.001 3000
# ratio[2] 0.293 0.034 0.229 0.271 0.292 0.315 0.361 1.001 3000
# ratio[3] 0.201 0.029 0.146 0.181 0.200 0.219 0.261 1.001 3000
# m.group[1] 10.186 0.115 9.970 10.110 10.180 10.260 10.420 1.002 1200
# m.group[2] 14.834 0.166 14.510 14.720 14.840 14.940 15.160 1.001 3000
# m.group[3] 20.178 0.192 19.810 20.050 20.180 20.300 20.560 1.001 3000
# sigma 1.114 0.065 0.995 1.067 1.110 1.155 1.248 1.001 3000
比較可能な混合比率、平均で見ると、いずれもいい推定をしているものの、ガンマ混合モデルのほうがより母集団に近い値を推定をした。母集団は正規分布であるのにもかかわらずである。また正規混合モデルの標準偏差の推定はズレが大きく、とくに真ん中のグループが狭められた。なお、DICは前者が707.2、後者が666.3なので、全体としてもガンマ混合モデルのほうが当てはまりがよさそうだ（こんなケースでDICの比較は有効なのか？？）。
DICだけでなく、全体としての母集団との当てはまりの良さを図示して比較し、絵合わせでもチェックしてみた。
実線が母集団の確率密度、黒の点線が正規混合モデルによる推定、赤の点線がガンマ混合モデル。正規混合モデルはデータの重なりにより推定が乱された様子がわかる。第3グループは正規混合モデルのほうがよく当てはまっているようだが、それ以外ではガンマ混合モデルのほうが母集団のラインにより近いものとなった。推定計算にかかる時間はガンマ混合モデルのほうが長いものの、今回のように平均にともなってバラツキが大きくなるようなデータについては、元の分布型にかかわらずより無難な推定を行なってくれることが期待できる。
なお、ここまで単一のヒストグラムを想定してきたが、実際のコホート解析では当然、時系列のデータとなる。発展形として平均や標準偏差を Y[c,t] = alpha + beta*Y[c-1,t-1] のような自己回帰モデルにすると前の時点の情報を利用してもう少し推定がしやすくなるかもしれない（希望的観測？）。だんだん欲が出てきました（c：コホート、t：時間）。
cf.その1で使用したMclustよりも汎用性の高い混合モデルを扱えるパッケージを見つけた。他の分布型も扱えるという大きな利点がある。ただし、推定精度はMclustの方がよいかもしれない。大きく狂うなら、分布型が違ってもMclustの方がよいだろう。そもそもこのパッケージは複数の量的形質を組み合わせた分離を主な用途と想定していると思う。
library(flexmix) # 要インストール
mm <- FLXMRglm(でーた ~ 1, family="gaussian") # ここで分布型を指定できる
mm3 <- flexmix(. ~ ., data=data.frame(でーた), k=3, model=mm, control=list(verbose=1, nrep=5))
# k=3 は分離するグループ数、Mclustと違って一括で比較できないので、検討すべきグループ数全部について行い（k=3, 2, 1）AIC、BIC を比較する必要がある。なお、こちらでは一般的な定義通り、値の小さいモデルがよりよいモデル。
summary(mm3) # 混合比率やAIC、BIC
parameters(mm3) # 平均と標準偏差
2013年8月9日金曜日
混合分布のパラメータ推定：その1（正規混合モデル）
（2024年03月27日追記）より強力なRパッケージ"mixR"が利用可能になっているようです。正規分布だけでなく、lognormal, Gamma, Weibull分布まで対応している模様（追記ここまで）
# 4部構成にしています。こちらのページはより基本形で正規混合モデルを使用しました。
# その1：正規混合モデル（オーソドックスな基本形）
# その2：ガンマ混合モデル（体長の増加に伴い分散も増加するモデル）
# その3：前後関係からユーザがグループ数を仮定して正規混合モデルを実行
# その4：体長と体幅など複数種類のデータで分ける（有効な最終手段）
*******************************************************
（注：mclustでの解析の場合、citation("mclust") で手法を使用の際の引用文献が得られる。mclustは私の作成ではありません）
似通った何かを識別する手段のひとつに統計学的な分離がある。一般的な問題の中にもたくさんあるはず…例えば、外見的な何かの特徴から不良品を識別して排除するなどだろう。沢山のデータを取って比較してみるとそういった傾向が見えてくるのは興味深い。
自分の専門分野では、野外の生物の生活史研究などでよくやる体サイズ頻度のヒストグラムでコホートを分離するとか、近縁種や雌雄間で何らかの量的形質によって両者を識別するといった用途がある。そんな時に判断材料が少なく、1 つの尺度でしか比較できないようなケースでも有効な方法を考えていた。しかし、これまでに見たことのあるやり方は、統計学的にはちょっと強引な手法であったり、自分でアルゴリズムを書かないといけないような敷居が高いものだった。また、提供媒体が古くなっていたりと現代的な解析環境では使いにくくなっているようにも感じたため、Rで実行できる解析方法を紹介することにした。
以下、コホートデータを統計学的に分離し、コホートを正規分布などに当てはめた時のパラメータ推定を試みた。仮定として、サイズ頻度分布上で隣り合うコホート同士が分離されるか否かは、コホート間の平均値が統計学的に差があれば異なるコホートとみなす、コホートの分散は成長の個体差にともなって増加する（分散の差異をコホートの違いの基準には用いない）、とした。
"mclustというR パッケージでも正規混合分布を解ける"という情報提供をいただき、使用法と使い勝手を試行錯誤してみた結果、ある程度の点推定まではこれだけで行けると判断しました。当初は最初からベイズ法で解くつもりでいましたが、ベイズ法の解析についても計算負荷を抑えるためmclust の結果を初期値として計算するように書き直しました。
この手法の有効性を確認する。まず、テストデータとして、体サイズの頻度分布を想定し、予め平均値と標準偏差を指定した 2 つの正規分布を母集団を作成、その母集団から 60, 40 サンプルずつの計 100 のサンプルを抽出し混合した。それぞれ平均 10, 14、標準偏差は 1.2, 1.4 にした。
### 母集団の設定
N <- 100 # サンプル全体のN数
mean1 <- 10
mean2 <- 14
sd1 <- 1.2
sd2 <- 1.4
### 母集団からのサンプリング、N1, N2 ずつ取ってきて混ぜた
ratio <- c(0.6, 0.4) # 6 : 4 の比率でサンプリングした
N1 <- N*ratio[1]
N2 <- N*ratio[2]
set.seed(1) # 乱数のタネを指定
Length2 <- c(rnorm(N1, mean1, sd1), rnorm(N2, mean2, sd2)) # 1, 2を混ぜたサンプル
N <- 100 # サンプル全体のN数
mean1 <- 10
mean2 <- 14
sd1 <- 1.2
sd2 <- 1.4
### 母集団からのサンプリング、N1, N2 ずつ取ってきて混ぜた
ratio <- c(0.6, 0.4) # 6 : 4 の比率でサンプリングした
N1 <- N*ratio[1]
N2 <- N*ratio[2]
set.seed(1) # 乱数のタネを指定
Length2 <- c(rnorm(N1, mean1, sd1), rnorm(N2, mean2, sd2)) # 1, 2を混ぜたサンプル
嫌な感じに混じりました…。
実線のカーブはそれぞれの母集団の正規分布の確率密度。
複数の正規分布が混ざった分布を、正規混合分布（normal mixture distribution）と言う。このサンプルを解析して元の母集団のパラメータを遡り推定し、正しく母集団を再現推定できるかどうかを検証した。パラメータとして、混合比率、各正規分布の平均、標準偏差を求めた。
まず、R パッケージの mclust を使用して大まかな推定をし、その結果をベイズ推定の初期値に使用した。mclustはEMアルゴリズムによるパラメータ推定をする。
（cf. コホート解析に混合正規分布を適用するレビューがあるとの情報提供を頂きました：赤嶺 2005. 混合正規分布とEMアルゴリズム. 水産海洋研究 69: 174–183）
library(mclust) # 要インストール（cf. 追加パッケージのインストール方法）
plot(mclustBIC(Length2))
# グループ分け数のモデルを情報量基準 BIC で比較（図示）
注：mclust では、BIC = 2*対数尤度 になっているようで、一般的な BIC の定義と正負が逆である。このため mclust の結果では BIC 最大のモデルがベストと判断。
E：各グループで分散は共通とするモデル、V：分散は異なるとするモデル。
グループ分け数 2 で最大 → 2 分割するのが妥当である。グループ数の推定は成功した。しかし、いずれのモデルでもグループ分けは2のモデルがベストとされたが、分散が等しいモデルが選ばれてしまうため、冒頭の仮定にしたがい分散が異なるモデルでパラメータを推定した。
cf. 標準偏差が倍くらいの違いがない限り分散が等しいモデルが選択されてしまった。
Mc <- densityMclust(Length2, modelNames="V")
# "V":分散が異なるモデル（cf. Mclustでもいいが、densityMclustの方が図示に便利）
# "V":分散が異なるモデル（cf. Mclustでもいいが、densityMclustの方が図示に便利）
mc <- summary(Mc, parameters=T) # 推定されたパラメータの一覧を見る
mc
# ----------------------------------------------------
# Gaussian finite mixture model fitted by EM algorithm
# ----------------------------------------------------
# Mclust V (univariate, unequal variance) model with 2 components:
# log.likelihood n df BIC ICL
# -213.4482 100 5 -449.9223 -461.0196
# Clustering table:
# 1 2
# 61 39
# Mixing probabilities:
# 1 2
# 0.6119736 0.3880264
# Means:
# 1 2
# 10.18958 14.18402
# Variances: # 標準偏差でなく分散
# 1 2
# 1.165315 1.879066
sqrt(mc$variance) # 標準偏差に直す
# 1 2
# 1.079497 1.370790
推定された確率密度をヒストグラム上に一括で図示できる（見栄えは悪いが推定の確認に便利）
plot(Mc, data=Length2)
#（hist()を用いているようで、breaksの項を加えればヒストグラムの区間を操作可能）
グループ数（G= ）、推定の初期値（prior=）を指定して実行することも可能、何らかの理由で事前情報があるのならば利用していいと思う。
Mc2 <- densityMclust(Length2, modelNames="V", G=3, prior=priorControl(mean=c(15, 21, 32)))
# それでも、相対的に小さいピークは、あるはずなのに認識されないことがある。そんなときはデータの範囲を絞って推定するといい。例えば、25より大きいデータに絞るなら、Length2[Length2 > 25] 平均値、標準偏差ともに母集団のものに近い値を推定できた。点推定でよいならばここまでで十分だろう。標準偏差の推定の粗さも許容するのならば、3グループ以上の場合でもmclustで済ませて構わないと思う。何よりもこの方法は容易で敷居が低いのがよいところだ。
**************************************************************
せっかくなのでmclustの結果を初期値として与え、さらにベイズ法で解析し、これらの値がどれくらいの範囲を取りうるのか、確率分布を求めてみた。
一般的なベイズ法のソフトウェア WinBUGS を用いた単純なコードで推定計算した。…とはいっても、多少の R と WinBUGS の経験がないと使いこなすのは難しいかもしれません。BUGSの推定がうまくいかなかった時に n.iter, n.burn, n.thin の値をどう増やしたらいいかなどは多少の勘が必要です。
コードはこちらの本を参考にした：「マルコフ連鎖モンテカルロ法（豊田秀樹、朝倉書店）」(p. 190~)。
# ベイズモデルの記述。推定計算に必要なデータ、初期値の設定。ディリクレ分布（値の合計が1になる）を2つの集団の混合比率の重み付けに利用することで正規混合モデル化した。なお重み付けをグループの選択に用いるところで、カテゴリカル分布を利用してグループ番号に変換した。
事前分布については、今回は事前情報があるので一般的なベイズ法の場合に比べ範囲を狭めていることに注意
（注：ここでは補助ツールWBUGSを用いて実行するため通常のWinBUGS使用と比べてパラメータの記述は単純化してある）
n.g <- 2 # グループ数の指定
data <- list(n.g=n.g, N=length(Length2), size=Length2, # 上記のデータ Length2 を使用
alpha=mc$pro, # mclustの推定した混合比率を使用
M=mc$mean) # mclustの推定した平均値を使用
par <- list(ratio=mc$pro, tau=1/mc$variance,F, # mclustの推定した混合比率、分散を使用
m.group=mc$mean, # mclustの推定した平均値を使用
group=sample(1:n.g, length(Length2), replace=T),F, # 初期値不明のためランダム
sigma=NA,delta.m=NA, delta.var=NA)
# 初期値が不要なパラメータ（他の数値計算で求まる）
### モデル式
model <- function() {
for (i in 1 : N) { # iループ
size[i] ~ dnorm(m.group[group[i]], tau[group[i]]) # 混合したサイズ分布（グループ毎に変化）
group[i] ~ dcat(ratio[1:n.g]) # グループ選択の事前分布、カテゴリカル分布を使用
} # iループ終わり
ratio[1:n.g] ~ ddirch(alpha[]) # 混合比率の事前分布：ディリクレ分布=ratioの合計が1になる
for (g in 1:n.g) { # gループ（n of group）
m.group[g] ~ dnorm(M[g], 0.5) # 各グループの平均値の事前分布
tau[g] ~ dgamma(1, 0.01) # 各グループの分散の逆数の事前分布
sigma[g] <- 1 / sqrt(tau[g]) # 同、標準偏差の事後分布の計算
} # gループ終わり
for (gg in 2:n.g) { # ggループ
delta.m[gg] <- m.group[gg] - m.group[gg-1] # 各グループの平均値の差の事後分布の計算
delta.var[gg] <- 1/tau[gg] - 1/tau[gg-1] # 各グループの分散値の差の事後分布の計算
} # ggループ終わり
} # model終わり
### 計算の実行
source("WBUGS.R") # cf. WBUGS、R2WinBUGSの補助ツール
seg <- wbugs(data, par, model, n.iter=110000, n.burn=10000, n.thin=100) # 計算実行
# 実行結果
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
# ratio[1] 0.588 0.112 0.311 0.541 0.606 0.657 0.753 1.035 1500
# ratio[2] 0.412 0.112 0.247 0.343 0.394 0.459 0.689 1.042 1300
# m.group[1] 10.194 0.239 9.761 10.050 10.185 10.330 10.680 1.004 2400
# m.group[2] 14.025 0.577 12.530 13.790 14.150 14.410 14.820 1.002 3000
# sigma[1] 1.085 0.207 0.678 0.972 1.079 1.194 1.483 1.028 3000
# sigma[2] 1.488 0.380 0.977 1.222 1.402 1.673 2.383 1.014 2300
# delta.m[2] 3.831 0.509 2.431 3.636 3.932 4.153 4.526 1.004 1300
# delta.var[2] 1.137 1.544 -0.901 0.179 0.789 1.691 5.109 1.001 3000
# deviance 332.039 21.662 308.097 317.400 325.800 339.400 390.100 1.007 2000
#（注：RhatはMCMCのchain間の収束基準、n.effはMCMCの有効サンプル数で自己相関が高いと低下する。ひとつでも Rhat > 1.1 になったり、n.eff < 1000 くらいの場合は推定があやしい。とくにRhatの方はよりシビアで、Rhatの基準が満たされていない場合は、n.iter, n.burn, n.thin を5~10倍に増やして再実行する必要がある）
# 中央値も表示（要は↑の50%の値）
seg$median
$ratio
[1] 0.6058 0.3942
$m.group
[1] 10.19 14.14
$sigma
[1] 1.078 1.414
$delta.m
[1] 3.925
$delta.var
[1] 0.82195
$deviance
[1] 325.8
# 推定値の事後確率分布を図示（MCMCサンプルの頻度分布）
par(mfrow=c(3,2))
hist(seg$sims.list$ratio[,1], main="ratio1")
hist(seg$sims.list$ratio[,2], main="ratio2")
hist(seg$sims.list$m.group[,1], main="mean1")
hist(seg$sims.list$m.group[,2], main="mean2")
hist(seg$sims.list$sigma[,1], main="sd1")
hist(seg$sims.list$sigma[,2], main="sd2")
data <- list(n.g=n.g, N=length(Length2), size=Length2, # 上記のデータ Length2 を使用
alpha=mc$pro, # mclustの推定した混合比率を使用
M=mc$mean) # mclustの推定した平均値を使用
par <- list(ratio=mc$pro, tau=1/mc$variance,F, # mclustの推定した混合比率、分散を使用
m.group=mc$mean, # mclustの推定した平均値を使用
group=sample(1:n.g, length(Length2), replace=T),F, # 初期値不明のためランダム
sigma=NA,delta.m=NA, delta.var=NA)
# 初期値が不要なパラメータ（他の数値計算で求まる）
### モデル式
model <- function() {
for (i in 1 : N) { # iループ
size[i] ~ dnorm(m.group[group[i]], tau[group[i]]) # 混合したサイズ分布（グループ毎に変化）
group[i] ~ dcat(ratio[1:n.g]) # グループ選択の事前分布、カテゴリカル分布を使用
} # iループ終わり
ratio[1:n.g] ~ ddirch(alpha[]) # 混合比率の事前分布：ディリクレ分布=ratioの合計が1になる
for (g in 1:n.g) { # gループ（n of group）
m.group[g] ~ dnorm(M[g], 0.5) # 各グループの平均値の事前分布
tau[g] ~ dgamma(1, 0.01) # 各グループの分散の逆数の事前分布
sigma[g] <- 1 / sqrt(tau[g]) # 同、標準偏差の事後分布の計算
} # gループ終わり
for (gg in 2:n.g) { # ggループ
delta.m[gg] <- m.group[gg] - m.group[gg-1] # 各グループの平均値の差の事後分布の計算
delta.var[gg] <- 1/tau[gg] - 1/tau[gg-1] # 各グループの分散値の差の事後分布の計算
} # ggループ終わり
} # model終わり
### 計算の実行
source("WBUGS.R") # cf. WBUGS、R2WinBUGSの補助ツール
seg <- wbugs(data, par, model, n.iter=110000, n.burn=10000, n.thin=100) # 計算実行
# 実行結果
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
# ratio[1] 0.588 0.112 0.311 0.541 0.606 0.657 0.753 1.035 1500
# ratio[2] 0.412 0.112 0.247 0.343 0.394 0.459 0.689 1.042 1300
# m.group[1] 10.194 0.239 9.761 10.050 10.185 10.330 10.680 1.004 2400
# m.group[2] 14.025 0.577 12.530 13.790 14.150 14.410 14.820 1.002 3000
# sigma[1] 1.085 0.207 0.678 0.972 1.079 1.194 1.483 1.028 3000
# sigma[2] 1.488 0.380 0.977 1.222 1.402 1.673 2.383 1.014 2300
# delta.m[2] 3.831 0.509 2.431 3.636 3.932 4.153 4.526 1.004 1300
# delta.var[2] 1.137 1.544 -0.901 0.179 0.789 1.691 5.109 1.001 3000
# deviance 332.039 21.662 308.097 317.400 325.800 339.400 390.100 1.007 2000
#（注：RhatはMCMCのchain間の収束基準、n.effはMCMCの有効サンプル数で自己相関が高いと低下する。ひとつでも Rhat > 1.1 になったり、n.eff < 1000 くらいの場合は推定があやしい。とくにRhatの方はよりシビアで、Rhatの基準が満たされていない場合は、n.iter, n.burn, n.thin を5~10倍に増やして再実行する必要がある）
# 中央値も表示（要は↑の50%の値）
seg$median
$ratio
[1] 0.6058 0.3942
$m.group
[1] 10.19 14.14
$sigma
[1] 1.078 1.414
$delta.m
[1] 3.925
$delta.var
[1] 0.82195
$deviance
[1] 325.8
# 推定値の事後確率分布を図示（MCMCサンプルの頻度分布）
par(mfrow=c(3,2))
hist(seg$sims.list$ratio[,1], main="ratio1")
hist(seg$sims.list$ratio[,2], main="ratio2")
hist(seg$sims.list$m.group[,1], main="mean1")
hist(seg$sims.list$m.group[,2], main="mean2")
hist(seg$sims.list$sigma[,1], main="sd1")
hist(seg$sims.list$sigma[,2], main="sd2")
全般的に左右非対称のため median を代表値として用います。
実行結果のdelta.m（平均値差）が2.5〜97.5%まですべて同じ符号のため、2つの集団の平均値に差があると言えた。ところが、delta.var（分散差）についてはゼロをまたいでしまいました。その確率をチェック：
dv <- seg$sims.list$delta.var
length(dv[dv<0]) / length(dv) # MCMCサンプル数のうち、負の値だった個数
0.1963333 、ベイズ法でもしっかり0.05を超えていました…。ここは冒頭の仮定を踏襲します！
ratioの真の値は0.6, 0.4に対し、推定値（中央値）は0.6058, 0.3942。mclustの推定より向上した。m.group（平均値）の真の値は10, 14としたが、推定値は10.19, 14.14でmclustと同等のよい推定だった。一方、sigma（標準偏差）の真の値は1.2と1.4としたが、推定値は1.078と1.414、第2グループの推定はmclustより向上したが、前者の方は真の値より小さいままだった。
推定された混合比率、平均値、標準偏差を用いた確率密度を点線で書き加えた。元の母集団の曲線をよく再現することができた。とくに第2グループの方は母集団の曲線にかなり近いものとなった。
第1グループの方は、むしろ推定曲線のほうが当てはまりがいいように見えた。どうやらオーバラップが激しいサンプルを用いたため、第2グループの裾野が第1グループの中心付近にまで入り込み押し上げた結果、真の平均の少し右側がピークになっており、そのピークを見た目通りに真のピークと見なして推定したためだろう。それにより第2グループの裾野は少し第1グループに食われたので、第2グループは少し右にシフトしたのだろう。このあたりが推定の限界のようだ。
なお、任意の値 X がどれくらいの確率で第 1 グループに属するかは、これで求まる。
length(Length2) * 0.6058 * dnorm(X, 10.19, 1.078) # N数 x 正規分布の確率密度関数
ここで length(Length2)*0.6058 は推定された第 1 グループの N 数。X は数字列でもOK、Length2 自身を試してみてもOK（注：X の部分は何かの数字とかに置き換えないと動作しない）。
登録:
コメント (Atom)