gistfile1.pyw

#!/usr/bin/env python
# -*- coding:utf-8 -*-
# 変分混合ガウス分布を解く。テストデータはOld Faithful
import math
from scipy import *
from scipy.special import polygamma
from scipy import linalg
from random import gauss

def main():
    x = read_file("faithful.txt")

    vgm = VGM( x, 0.01, 0.01, 25, 25, 0, 0, 1, 1, 2, 2)
    vgm.train()
    vgm.plot()

def read_file(filename):
    f = open(filename, "r")
    x = []
    for i in f.readlines():
        #print "i",i
        x.append( map(lambda x: float(x), i[:-1].split() ) )
    return x


class VGM():
    def __init__(self, x, alpha0, alphak, beta0, betak,
                    m0, mk, W0, Wk, nu0, nuk, k = 6):
        # 値の初期化
        # m[k] も D次元ベクトル
        self.D = D = len(x[0])
        self.N = N = len(x)
        self.x = mat(x).T # x[n] = [x1, ... xd]^T
        self.alpha0 = alpha0
        self.beta0 = beta0
        #self.m0 = mat( [gauss(0,1) for i in range(D) ]).T
        self.m0 = mat( [0 for i in range(D) ]).T
        self.W0 = mat([[0.01 if i == j else 0 for i in range(D)] for j in
            range(D)])
        self.nu0 = nu0
        self.K = k

        # 初回は計算出来ないので設定
        # リストの要素は(0..k-1)なので-1する
        self.alpha = [alphak]*k
        self.beta = [betak]*k
        self.m = [mat( [gauss(0,1) for i in range(D) ]).T for j in range(k)]
        self.W = [ self.W0 ]* k
        self.nu = [nuk]*k
        self.gamma = [ [0]*k]*self.N
        self.n = [0] *k


    def train(self):
        for i in range(20):
            self.vb_e_step()
            self.vb_m_step()


   def vb_e_step(self):
        #x[n] はD次元ベクトル
        alpha_hatt = sum([ self.alpha[k] for k in range(0, self.K)])
        rho = [[0]*self.K]*self.N

        for k in range(0, self.K):

            # (10.65)
            E_ln_lambda = sum([ self.psi( (self.nu[k]+1-i)/2. ) for i in\
                range(0,self.D) ]) + self.D*math.log(2.) + \
                math.log( linalg.det(self.W[k]) )

            # (10.66)
            E_ln_pi = self.psi(self.alpha[k]) - self.psi(alpha_hatt)

            sum_rho = [0] * self.N
            for n in range(0,self.N):
                # (10.64)
                E_mu_lambda = self.D/self.beta[k] +\
                    self.nu[k]*((self.x.T[n] - self.m[k].T)*self.W[k]) \
                    *(self.x.T[n].T - self.m[k])

                # (10.46)
                ln_rho = E_ln_pi + E_ln_lambda/2. \
                    -self.D/2.*math.log(2*math.pi) -1./2.*E_mu_lambda

                rho[n][k] = math.exp(ln_rho)
                sum_rho[n] += rho[n][k]

        # (10.49) 負担率
        for k in range(0, self.K):
          for n in range(0, self.N):
            self.gamma[n][k] = rho[n][k]/sum_rho[n]


    def vb_m_step(self):

      for k in range(0, self.K):
        # (10.51) 
        self.n[k] = nk = sum([self.gamma[n][k] for n in range(0, self.N)])

        # (10.52) 
        x_bark = mat([sum([self.gamma[n][k] * self.x.T[n].T[d] for n in range(0,
            self.N)]) / nk for d in range(0,self.D)]).T

        # (10.53)
        sk = sum([self.gamma[n][k]*(self.x.T[n].T - x_bark)*(self.x.T[n] -
            x_bark.T)]) / nk

        # (10.58) 
        self.alpha[k] = self.alpha0 + nk

        # (10.60) 
        self.beta[k] = self.beta0 + nk
        
        # (10.61)
        self.m[k] = (self.beta0*self.m0 + nk*x_bark)/self.beta[k]
        
        # (10.62) 
        Wi = self.W0.I + nk*sk + self.beta0*nk/self.beta[k] \
            * (x_bark - self.m0)*(x_bark.T - self.m0.T)
        self.W[k] = Wi.I
        
        # (10.63)
        self.nu[k] = self.nu0 + nk
        

    @staticmethod
    def psi(x):
        # ディガンマ関数psi(x)
        return polygamma(0,x)
        

    def plot(self):
        import matplotlib
        import pylab
        for n in range(self.N):
            x = self.x.T[n].tolist()[0]
            matplotlib.pyplot.plot( x[0], x[1], 'go')
            
        for k in range(self.K):
            m = self.m[k].T.tolist()[0]
            print m
            if self.n[k] > 0.001:
                matplotlib.pyplot.plot( m[0], m[1], 'ro')
                
        matplotlib.pyplot.show()
        
if __name__ == "__main__":
    main()