Matplotlibで混合ガウス分布の確率密度関数をプロットする

PythonとNumpy/Scipyの練習．前回はNumpyを使って混合ガウス分布のEMアルゴリズムを実装で混合ガウス分布について取り扱ったので，今回は混合ガウス分布についての数式をおさらいしつつ，確率密度関数をプロットしようと思う．

混合ガウス分布の概要

混合ガウス分布は，複数のガウス分布の線形結合で表すことができる．

p(\bold{x})=\sum_{k=1}^{K}\pi_k\mathcal{N}(\bold{x}|\bold{\mu}_k,\bold{\Sigma}_k)

ここで，\pi_kを混合係数，\mathcal{N}(\bold{x}|\bold{\mu}_k,\bold{\Sigma}_k)を混合要素と呼ぶ．混合系数は確率の条件\sum_{k=1}^{K}\pi_k = 1 および0 \le \pi_k < 1 を満たす．

また，混合要素\pi_k = p(k)をk番目の混合要素が選ばれる事前確率とし，\mathcal{N}(\bold{x}|\bold{\mu}_k,\bold{\Sigma}_k) をkが与えられた時のx の条件付き密度とすると，x の周辺分布は

p(\bold{x})=\sum_{k=1}^{K}p(k)p(\bold{x}|k)

で表すことができる．このk の選択が，EMを用いたGMMでの隠れ変数に対応している．

一次元の混合ガウス分布

1
2
3
4
5
6
7
8
9
10
11
12
13

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

x = np.linspace(-5, 5,200)
pi_k = np.array([0.3,0.7])
norm1 = mlab.normpdf(x, -1, 1)
norm2 = mlab.normpdf(x, 2, 1)

plt.plot(x,pi_k[0]*norm1,color="blue")
plt.plot(x,pi_k[1]*norm2,color="blue")
plt.plot(x,pi_k[0]*norm1+pi_k[1]*norm2, color="red")
plt.show()

gmm_1.png

今回は正規分布の確率密度関数にmatplotlibのmlab.normpdfを使ったが，scipy.statsのnorm.pdfを使うこともできる．

1
2
3
4

import numpy as np
rom scipy import stats
norm1 = stats.norm.pdf(x, loc=-1, scale=1)
norm2 = stats.norm.pdf(x, loc=2, scale=1)

二次元の混合ガウス分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

pi_k = np.array([0.3,0.7])
x = np.linspace(-3.0, 3.0,200)
y = np.linspace(-2.0, 2.0,200)
X, Y = np.meshgrid(x, y)
Z1 = mlab.bivariate_normal(X, Y, 0.5, 1.0, -1, -1)
Z2 = mlab.bivariate_normal(X, Y, 1.0, 1.0, 1, 1)
Z = pi/k[0]*Z1 + pi_k[1]*Z2

CS = plt.contour(X, Y, Z)
plt.clabel(CS, inline=1, fontsize=10)
plt.show()

gmm_2.png

参考