GMMとEMアルゴリズム | Python + Numpy
GMM混合ガウス分布MixtureGaussianModel
混合ガウス分布と多変量ガウス分布は違うものだよ。EMは実装が容易なので、手を動かすとすぐに理解できます。参考資料は自分で探すこと。
実装例
#! coding:utf-8 """ gmm.py 混合ガウス分布のためのEMアルゴリズム http://aidiary.hatenablog.com/entry/20100521/1274445325 """ import numpy as np from pylab import * # 混合ガウス分布の数(固定) K = 2 def scale(X): """ データ行列Xを属性ごとに標準化したデータを返す :param X: :return: """ # 属性の数(=列の数) col = X.shape[1] # 属性ごとに平均値と標準偏差を計算 mu = np.mean(X, axis=0) sigma = np.std(X, axis=0) # 属性ごとにデータを標準化 for i in range(col): X[:, i] = (X[:, i] - mu[i]) / sigma[i] # TODO:Pysonicに改善あり return X def gaussian(x, mean, cov): """ 多変量ガウス分布 :param x: 観測データ :param mean: 平均 :param cov: 分散 :return: """ A = 1. / ((2. * np.pi) ** (x.size / 2.0)) B = 1. / (np.linalg.det(cov) ** 0.5) C = -0.5 * np.dot(np.dot(x - mean, np.linalg.inv(cov)), x - mean) return A * B * np.exp(C) def likelihood(X, mean, cov, pi): """ 対数裕度関数 :param X: 観測データ :param mean: 平均 :param cov: 分散 :param pi: 混合計数 :return: 尤度 """ sum = 0. for n in range(len(X)): temp = 0.0 for k in range(K): temp += pi[k] * gaussian(X[n], mean[k], cov[k]) sum += np.log(temp) return sum def main(): # 訓練データをロード data = np.genfromtxt("faithful.txt") X = data[:, :] # データの標準化 x = scale(X) # データ数 N = len(x) # --------------------------------------------- # 訓練データから混合ガウス分布のパラメータを # EMアルゴリズムで推定する # --------------------------------------------- # 平均、分散、混合計数を初期化 mean = np.random.rand(K, 2) cov = np.zeros((K, 2, 2)) for k in range(K): cov[k] = [[1., 0.], [0., 1.]] pi = np.random.rand(K) # 負担率の空配列を用意 gamma = np.zeros((N, K)) # 対数裕度の初期値を計算 like = likelihood(x, mean, cov, pi) print 'like: ', like turn = 0 MAX_ITERATE = 100 THR_VALUE = 0.01 for l in range(MAX_ITERATE): print 'like[%d] = %0.1f' % (l, like) # E-step : 現在のパラメータを使って、負担率を計算 for n in range(N): # 分母はkによらないので最初に1回だけ計算 denominator = 0.0 for j in range(K): denominator += pi[j] * gaussian(x[n], mean[j], cov[j]) # 各kについて負担率を計算 for k in range(K): gamma[n, k] = pi[k] * gaussian(x[n], mean[k], cov[k]) / denominator # M-step : 現在の負担率を使って、パラメータを再計算 for k in range(K): # Nkを計算 Nk = 0.0 for n in range(N): Nk += gamma[n, k] # 平均を再計算 mean[k] = np.array([0., 0.]) for n in range(N): mean[k] += gamma[n, k] * x[n] mean[k] /= Nk # 共分散を再計算 cov[k] = np.array([[0., 0.], [0., 0.]]) for n in range(N): # おそらくx[n,:] - mean[k] (1xK)行列 temp = x[n] - mean[k] # おそらく cov[k] += gamma[n,k] * np.dot(tmep.T, temp) cov[k] += gamma[n, k] * np.matrix(temp).reshape(2, 1) * np.matrix(temp).reshape(1, 2) cov[k] /= Nk # 混合計数を再計算 pi[k] = Nk / N # 収束判定 new_like = likelihood(x, mean, cov, pi) diff = new_like - like if diff < THR_VALUE: break like = new_like # # ガウス分布の平均を描画 for k in range(K): scatter(mean[k, 0], mean[k, 1], c='r', marker='o') # 等高線を描画 xlist = np.linspace(-2.5, 2.5, 50) ylist = np.linspace(-2.5, 2.5, 50) gx, gy = np.meshgrid(xlist, ylist) for k in range(K): gz = bivariate_normal(gx, gy, np.sqrt(cov[k, 0, 0]), np.sqrt(cov[k, 1, 1]), mean[k, 0], mean[k, 1], cov[k, 0, 1]) cs = contour(gx, gy,gz, 3, colors='k', linewidths=1) # 訓練データを描画 plot(x[:, 0], x[:, 1], 'gx') grid(True) savefig('gmm.png') show() pass if __name__ == "__main__": main()
結果
使用しているデータセット
faithful.txt
その他データセットが欲しい方は
http://naoyat.hatenablog.jp/entry/2011/12/29/064419
ゼロから作るDeep Learning ―Pythonで学ぶディープラーニングの理論と実装
- 作者: 斎藤康毅
- 出版社/メーカー: オライリージャパン
- 発売日: 2016/09/24
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (6件) を見る
退屈なことはPythonにやらせよう ―ノンプログラマーにもできる自動化処理プログラミング
- 作者: Al Sweigart,相川愛三
- 出版社/メーカー: オライリージャパン
- 発売日: 2017/01/26
- メディア: 単行本(ソフトカバー)
- この商品を含むブログを見る
【Amazon.co.jp限定 】 ASUS ゲーミングモニター 23型フルHDディスプレイ (応答速度1ms / HDMI×2,D-sub×1 / スピーカー内蔵 / 3年保証) VX238H-P
- 出版社/メーカー: Asustek
- 発売日: 2014/05/28
- メディア: Personal Computers
- この商品を含むブログを見る
- 作者: 裕時悠示
- 出版社/メーカー: SBクリエイティブ
- 発売日: 2016/11/17
- メディア: Kindle版
- この商品を含むブログを見る
- 作者: 丸戸史明,深崎暮人
- 出版社/メーカー: KADOKAWA
- 発売日: 2016/11/19
- メディア: 文庫
- この商品を含むブログ (2件) を見る
- 作者: 近江のこ
- 出版社/メーカー: 小学館
- 発売日: 2016/11/18
- メディア: Kindle版
- この商品を含むブログを見る