Rで混合ガウスモデル.1

記事の目的

1次元の混合ガウスモデルのパラメータ推定をEMアルゴリズムで行います。
個人的にはPRMLの解説が分かりやすかったです。
参考1：しくみがわかるベイズ統計と機械学習
参考2 ：[PRML下] (https://www.amazon.co.jp/パターン認識と機械学習-下-ベイズ理論による統計的予測-C-M-ビショップ/dp/4621061240/ref=pd_bxgy_img_2/358-0979251-6826126?_encoding=UTF8&pd_rd_i=4621061240&pd_rd_r=2d3fdac6-bc15-44fc-bbb7-2f912b059365&pd_rd_w=LVErD&pd_rd_wg=QOmz6&pf_rd_p=e64b0a81-ca1b-4802-bd2c-a4b65bccc76e&pf_rd_r=EP4VHHZ9R2230G5C5R8H&psc=1&refRID=EP4VHHZ9R2230G5C5R8H)

目次

1. モデルの説明
2. 使用データ
3. EMアルゴリズム実装
4. 結果

1. モデルの説明

2. 使用データ

10歳の100人分の身長のデータと、20歳の200人分の身長のデータを誤って同時に計測してしまった。

library(dplyr)
set.seed(100)
age10 <- rnorm(100, 140, 5)
age20 <- rnorm(200, 170, 7)
X <- append(age10, age20) 
N <- length(X)
X %>% hist()

3. EMアルゴリズム実装

# (1)
K <- 2
pi <- rep(1/K, K)
u <- rnorm(K, mean(X), sd=10)
sigma <- rep(sd(age)/2, K)

# EMアルゴリズムの計算途中で使用する変数
max_iter <-30
pi.N <- matrix(NA, N, K)
r <- matrix(NA, N, K)
R <- {}

for(s in 1:max_iter){
  # (2) E-step
  for(i in 1:N){
    for(k in 1:K){
      pi.N[i, k] <- pi[k]*dnorm(X[i], u[k], sqrt(sigma[k]))
    }
    r[i, ] <- pi.N[i, ]/sum(pi.N[i,])
  }
  # (3) M-step
  R <- apply(r, 2, sum)
  pi <- R/sum(R)
  for(k in 1:K){
    u[k] <- sum(r[,k]*X)/R[k]
    sigma[k] <- sum(r[, k]*(X-u[k])^2)/R[k]
  }
}

4. 結果

z <- apply(r, 1, which.max)
par(mfrow=c(1,1))
hist(X[z==1], xlim = c(120,200), col = "red")
hist(X[z==2], col = "blue", add=T)