混合分布のパラメータ推定：その３（グループ数をユーザ側で指定する場合のコード…）

（このページは URL を2013ベントス・プランクトン学会大会の要旨に掲載したものからリンクしています）

# ４部構成にしてあります、ここではグループ数の推定がうまくいかない場合の対処について。
# より基本形の例は、その１、２の方でまとめてあります。
# その１：正規混合モデル（オーソドックスな基本形）
# その２：ガンマ混合モデル（体長の増加に伴い分散も増加するモデル）
# その３：前後関係からユーザがグループ数を仮定して正規混合モデルを実行
# その４：体長と体幅など複数種類のデータで分ける（有効な最終手段）

*******************************************************

（その１からの続きです）
体サイズ頻度分布のヒストグラムでよくあることですが、高齢のコホートが他と比べてかなり少ない場合には、解析すると隣のグループに含められてしまうことが多々ある（想定するグループ数と異なると感じる場合…それは主観的かもしれないが）。


もし、前の時点でのグループ数などの事前情報によって、そこに確かにコホートが存在するはずという仮定ができるのならば、グループ数は既知であるとして解析してもよいだろう。完全に仮定無しでの解析はなかなか難しい…混合分布問題は完全解決には至っていないと思う。


その２でやり方を書いたような、グループ数指定（G=◯の項を加える）をしても、途中にもう一つグループを作られて、肝心の最高齢集団が無視されたりとなかなか言うことを聞いてくれない場合がある。

その２のガンマ混合モデルならば解決できるだろうが、ベイズを使うのは敷居が高いだろうし計算時間もかかる。ここでは mclust（正規混合モデル）を用いた別の解決方法を紹介する。

mclust （正規混合モデル）を使用して、グループ数をユーザ側で決めた場合の混合分布モデリングを行う。一番若い齢（体長などの値の小さい）のコホートから順に推定し、推定し終わったら除いて残りのデータで解析するというステップ・バイ・ステップ法で処理することによって、集団サイズの小さい最高齢のコホートも検出することができる。

### 母集団に用いるサンプルデータ

N <- 200 # サンプル全体のN数
mean1 <- 10
mean2 <- 15
mean3 <- 20
sd1 <- 1.0
sd2 <- 1.2
sd3 <- 1.4

### 母集団からのサンプリング、N1, N2, N3 ずつ取ってきて混ぜた
ratio <- c(0.5, 0.3, 0.2) # 5:3:2 の比率でサンプリングした
 N1 <- N*ratio[1]
 N2 <- N*ratio[2]
 N3 <- N*ratio[3]

set.seed(1) # 乱数のタネを指定
Length3 <- c(rnorm(N1, mean1, sd1), rnorm(N2, mean2, sd2), rnorm(N3, mean3, sd3)) # 1, 2, 3を混ぜたサンプル

### mclust による正規混合モデリングを最若齢〜最高齢のコホートへと順に１つずつ推定するプログラム（いったん読み込ませたら、Rを終了するまで何度でも再利用できる）

step.mclust <- function(data, n.group, var="V", priors=NULL) {
require(mclust)
x.min <- 5*floor(0.2*min(data))
x.max <- 5*ceiling(0.2*max(data))
x <- seq(x.min, x.max, by=(x.max - x.min)/100)
summ <- numeric(0)
for (g in 1:n.group-1) {
mc <- summary(Mclust(data, G=n.group - g, modelNames=var), parameters=T)
n <- round(mc$pro[1]*length(data)) # グループ１の推定個体数
n2 <- round(mc$pro[2]*length(data)) # グループ２の推定個体数
m <- mc$mean[1] # グループ１の推定平均値
m2 <- mc$mean[2]
s <- sqrt(mc$variance[1]) # グループ１の推定標準偏差
s2 <- sqrt(mc$variance[2])
# １の中では確率密度が小さく、かつ２の方が確率密度が大きいものは２である
# １の中では確率密度が小さいが、２よりも確率密度が大きいものも１である
d <- dnorm(data, m, s) # dataがグループ１に属する確率密度
d2 <- dnorm(data, m2, s2) # 同、グループ２に…
dd <- d > d2 # グループ２よりもグループ１に属する確率が高い個体
p.rank <- rank(1/d) # 確率密度の高い順
gr <- data[dd ==T & p.rank <= n] # グループ１に属する可能性が高い個体（でもまだグループ１が残ってる）
# 面倒でも二段階抽出をしないとグループ１の左の裾が落ちてしまう
gr2 <- setdiff(data, gr) # dataからgrを除いた残り（この中にグループ1が残っている）
r.remain <- rank(1/dnorm(gr2, m, s)) # 残り物の中で、グループ１に属する確率が高いものの順
remain <- gr2[r.remain <= n - length(gr)] # グループ１の残党
gr <- c(gr, remain) # remainをグループ１に追加
data <- setdiff(gr2, remain) # グループ２からremainを除外したものを次のループのdataにする
summ <- rbind(summ, data.frame(n, m, s)) }
rownames(summ) <- 1:n.group
return(summ) }

### 使用法
step.mclust(Length3, 3) # データ、グループ数を指定するだけ

### 出力結果（n：各グループの推定 N数、m：同 平均値、s：同 標準偏差）

#      n        m         s
# 1  99 10.09472 0.8822265
# 2  60 14.68063 1.4391196
# 3  41 19.91919 2.3426576

# 最後のグループの標準偏差の推定が良くない。
# また、もっと最後のグループの集団サイズを減らすとけっこう推定は悪くなる…。とくに集団サイズの推定がよくない。
そういうケースにもきっちりと対応したい場合には、その２のガンマ混合モデルでがんばるなど試行錯誤するしかないだろう。