１次元混合正規分布の等高線の描き方

書籍、「続・わかりやすいパターン認識」のｐ182の図9-3の「混合正規分布のパラメータ推定」
の図でμ1とμ2の等高線を描きたいのですが、この図のμ1とμ2の関係式はどのようにして導かれる
のでしょうか？ｐ181の図9.3１の説明では、式（9.17）により計算される対数尤度logp(x;θ)の
等高線が描かれているとあります。

この書籍勉強されている方で、分かる方、御教示でがえればと思います。

No.3 回答者： kamiyasiro 回答日時： 2017/03/22 23:14

#2です。

他のみなさんのために、テキストの記述の概略を示します。

１次元混合分布を考えます。
N( 3,1^2)に従う300個と
N(-1,1^2)に従う200個の混合分布です。
今、μ1，μ2を未知数として、500個の実データから決めたいのです。
そのために、μ1，μ2を振りながら、対数尤度が最大になるところを探します。

このシーンに付随して、尤度の等高線が示されており、
それがご質問者の言うｐ182の図9.3です。
(#2のスクリプトが描く図です)

この図は紙面に垂直にy軸（対数尤度）があり、
その大きさが等高線になっています。
横軸は、μ1，μ2の水準を振っています。

ここからは、ご質問に対する回答です。

yの値は、μ1，μ2の関数ではありません。
関係式は?と書かれているので、誤解されていると思います。
yの値は、500個の乱数発生させたデータに依存して変化します。

対数尤度は、データ500個の同時確率の値です。
それは、式(9.17)のΣlog p(x;θ)ですが、この式のp(x;θ)は、
式(9.4)のΣpi p(x;θ)です。piは、この場面では、
μ1のクラスが0.6、μ2のクラスが0.4です。

上の２つのサンメイションΣの違いはお分かりですよね。
上のは、サンプル500個の和
下のは、2クラスですので、2個の和です。

まず、個々のサンプルxの確率値p(x;θ)を求めます。
θは、分布のパラメータだから、μとσです。
Rでは、dnorm(x,μ,σ)で計算できます。
この値は、正規分布曲線の高さに相当します。

式(9.4)を実際に式にすると、２つのクラスの和になります。
dnorm(x,mu1,1) * Pr1 + dnorm(x,mu2,1) * Pr2

この対数の和が対数尤度で、式(9.17)です。実際に式にすると、
sum(log(dnorm(x,mu1,1) * Pr1 + dnorm(x,mu2,1) * Pr2))

これをμ1，μ2の水準を振りながら計算したいので、muを配列に置き換え、
sum(log(dnorm(x,z[i,1],1) * Pr1 + dnorm(x,z[i,2],1) * Pr2))
これをforループで回します。

あとは、各格子点のデータを等高線としてプロットすればよいのです。

#2のスクリプトでは、ついでに尤度が最大値を取るときの
μ1，μ2を求めています。そのときの対数尤度の値も出力しています。
おおかた183ページの上から3行目の記述に合いますので、
是非確認して下さい。

この回答へのお礼

等高線は解析関数（y = f(x)）の形）では求められないということですよね。
Rは少しかじった程度ですが、本があるのでこれを元にPythonに落としてみます。
ありがとうございました。

お礼日時：2017/03/27 12:54

No.2 ベストアンサー 回答者： kamiyasiro 回答日時： 2017/03/22 22:12

企業で統計を推進する立場にある者です。

まず、等高線を描くRスクリプトを投稿します。
機械学習を学ばれている方なら、Rくらいは使えますよね。

次の投稿で、対数尤度の部分の解説をしたいと思います。

～～～～～～～～～～～～～～～～～～～～～～

rm(list=ls())
par(ask=T)

# 1次元混合モデルの尤度の等高線を描く
# μ1,μ2は未知数としてパラメータ扱いする

n <- 500
Pr1 <- 0.6
Pr2 <- 0.4

# 500個のデータを乱数生成しヒストグラムを描く

x <- c(rnorm(n * Pr1,3,1),rnorm(n * Pr2,-1,1))
bk <- seq((-5-1/6),(7+1/6),by=1/3) # 図9.2の区切り方
hist(x,breaks=bk,xlim=c(-4,6))

# 等高線図を作るためのグリッドを生成する

mu1 <- mu2 <- seq(-4,6,by=0.2)
z <- expand.grid(mu1,mu2)
len <- nrow(z)

# 500個のデータの対数尤度を、格子点毎に計算する

y <- NULL
for(i in 1:len){
l <- sum(log(dnorm(x,z[i,1],1) * Pr1 + dnorm(x,z[i,2],1) * Pr2))
y <- append(y,l)
}

# 格子点データを等高線図，外観図として描く

y <- matrix(y,ncol=sqrt(len))
contour(mu1,mu2,y,nlevels=50,drawlabels=FALSE)
persp(mu1,mu2,y,theta=-20,phi=45,expand=0.5,col="lightblue",shade=0.75)

# 最尤値

index <- which(y == max(y))
z[index,]
max(y) # そのときの対数尤度

No.1 回答者： rabbit_cat 回答日時： 2017/03/18 21:58

さすがに、特定の本を持っていることを前提に、具体的な内容を全く書かない質問っていうのは、ちょっとどうかと思います。

（まあ、有名な本ではあるんで、その本は私自身は持ってはいますが。。）