EMアルゴリズムを用いた混合ガウス分布モデルのクラスタリング処理が実装できない

Question

現在以下のサイトに載っているEMアルゴリズムを用いた混合ガウス分布モデルのプログラムを使ってRGB座標の画素値をクラスタリングしようと考えています。(全体コードのを使っています)
https://tips-memo.com/python-emalgorithm-gmm#i-17

import numpy as np
import csv
from numpy import linalg as la
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import sys

txt_dir = sys.argv[1]
csv_name = sys.argv[2]
dat_name = sys.argv[3]

with open(txt_dir) as f:
    reader = csv.reader(f)
    l = [row for row in reader]
    for i in range(len(l)):
      for j in range(len(l[i])):
        l[i][j] = float(l[i][j])
l = np.array(l)

def calc(x, mu, sigma_inv, sigma_det):
  D = x.shape[0]
  exp = -0.5*(x - mu).T@sigma_inv.T@(x - mu)
  denomin = np.sqrt(sigma_det)*(np.sqrt(2*np.pi)**D)
  return np.exp(exp)/denomin

def gauss(X, mu, sigma):
  output = np.array([])
  eps = np.spacing(1)
  Eps = eps*np.eye(sigma.shape[0])
  sigma_inv = la.inv(sigma)
  sigma_det = la.det(sigma)
  N = X.shape[0]
  for i in range(N):
    output = np.append(output, calc(X[i], mu, sigma_inv, sigma_det))
  return output

def mix_gauss(X, Mu, Sigma, Pi):
  k = len(Mu)
  output = np.array([Pi[i]*gauss(X, Mu[i], Sigma[i]) for i in range(k)])
  return output, np.sum(output, 0)[None,:]

def setInitial(X, k):
  D = X.shape[1]
  Mu = np.random.randn(k, D)
  Sigma = np.array([np.eye(D) for i in range(k)])
  Pi = np.array([1/k for i in range(k)])
  return Mu, Sigma, Pi

def log_likelihood(X, Mu, Sigma, Pi):
  K = Mu.shape[0]
  D = X.shape[1]
  N = X.shape[0]
  _, out_sum = mix_gauss(X, Mu, Sigma, Pi)
  logs = np.array([np.log(out_sum[0][n]) for n in range(N)])
  return np.sum(logs)

def EM(X, k, Mu, Sigma, Pi, thr):
  K = Mu.shape[0]
  D = X.shape[1]
  N = X.shape[0]
  log_list = np.array([])
  log_list = np.append(log_list, log_likelihood(X, Mu, Sigma, Pi))
  count = 0
  while True:
    out_com, out_sum = mix_gauss(X, Mu, Sigma, Pi)
    gamma = out_com / out_sum
    N_k = np.sum(gamma, 1)[:,None]
    Mu = (gamma @ X) / N_k

    sigma_list = np.zeros((N, K, D, D))
    for k in range(K):
      for n in range(N):
        sigma_com = gamma[k][n]*(X[n] - Mu[k])[:,None]@((X[n] - Mu[k]))[None,:]
        sigma_list[n][k] = sigma_com
    Sigma = np.sum(sigma_list, 0) / N_k[:,None]
    Pi = N_k/N
    
    log_list = np.append(log_list, log_likelihood(X, Mu, Sigma, Pi))
    if np.abs(log_list[count] - log_list[count+1]) < thr:
      return count+1, log_list, Mu, Sigma, Pi, gamma
    else:
      print("Previous log-likelihood gap:" + str(np.abs(log_list[count] - log_list[count+1])))
      count += 1

thr = 0.01
k = 4
Mu, Sigma, Pi = setInitial(l, k)
n_iter, log_list, Mu, Sigma, Pi, gamma = EM(l, k, Mu, Sigma, Pi, thr)
print("Iteration:"+str(n_iter))
print("log-likelihood:"+str(log_list))

params = np.array([Mu.ravel(), Sigma.ravel(), Pi.ravel()])
index = np.argmax(gamma, 0)

cm = plt.get_cmap("tab10")
fig = plt.figure()
ax = Axes3D(fig)
N = l.shape[0]

for n in range(N):
  ax.plot([l[n][0]], [l[n][1]], [l[n][2]],  "o", color=cm(index[n]), ms=1.5)
ax.view_init(elev=30, azim=45)
plt.show()

with open(csv_name, 'w') as file:
    writer = csv.writer(file, lineterminator='\n')
    writer.writerows(gamma.T)

with open(dat_name, 'w') as file:
    writer = csv.writer(file, lineterminator='\n\n')
    writer.writerows(params)

ここに書かれている10000個の3次元データがどのような形で保存されているかわからず。試しに右から各点のRGB座標をランダムに置きました。(全部で9つの点です。)

しかし以下のように実装しても以下のようにPrevious log-likelihood gap:nanがずっと出てしまいURLに載っているようなクラスタリング処理が上手くいきません。(ちなみにnew1.csvは先ほどの9つのランダム座標でh.csvとp.csvはそれぞれ負担率とパラメータを入れるための空のcsvファイルにしています。)

C:\pleiades\workspace\GMM>python program.py new1.csv h.csv p.csv
program.py:55: RuntimeWarning: divide by zero encountered in log
  logs = np.array([np.log(out_sum[0][n]) for n in range(N)])
program.py:67: RuntimeWarning: invalid value encountered in true_divide
  gamma = out_com / out_sum
Previous log-likC:\pleiades\python\3\lib\site-packages\numpy\linalg\linalg.py:2159: RuntimeWarning: invalid value encountered in det
  r = _umath_linalg.det(a, signature=signature)elihood gap:nan
Previous log-likelihood gap:nan
Previous log-likelihood gap:nan
Previous log-likelihood gap:nan

Comments

[meg_]

> しかし以下のように実装しても以下のようにエラーが出てしまいます。
メッセージはWarningですが、掲載されているものとは別にエラーが出ているのでしょうか？そうであればそのエラーメッセージを追記してください。

[mainichiramen]

すみませんエラーではなく上手く実装されずPrevious log-likelihood gap:nanがずっと出てしまいます。
書き直します

Answer 1

out_sumについて2種類のRuntimeWarningが出ています。
program.py:55: RuntimeWarning: divide by zero encountered in log
logs = np.array([np.log(out_sum[0][n]) for n in range(N)])
program.py:67: RuntimeWarning: invalid value encountered in true_divide
gamma = out_com / out_sum
logの引数であることから、本来はout_sumは0より大きい正の値であるはずなのに、そうなってはいないのではないでしょうか？log(0)だとnanが返されるので、Previous log-likelihood gapの値と関係ありそうです。
関数log-likelihoodの中のout_sumの値やout_sumを計算しているmix_gaussの引数の値などが正しく設定されていないかもしれません。print文を使ってそれらの値を出力し、想定通りの値になっているか確認してはどうでしょう？いまのところ、どのような値が出力されますか？