Wednesday, December 15, 2010

Bilinen Varyasyonlar ve EM Algoritmasi

Onceki bir yazida 2 boyutlu birkac Gaussian karisimini gosterdik. Bu yazida biri 100 digeri 20 tane olmak uzere iki Gaussian'dan uretilen noktalari nasil modellenebilecegini gosterecegiz. Bu analiz icin genellikle kullanilan yontem EM yani Expectation Maximization algoritmasidir. Uretilen veri icinde "asiri" degerlerin oldugu ama geri kalani normal dagilima sahip bir veriyi modelleyecegiz. Asiri degerler nereden gelebilir? Olcum aletinde bir hata olabilir, arada sirada, yanlislikla asiri degerler kaydedilebilir. EM'e bu asiri degerlerin belli bir alanda odaklanacagini ongorerek o noktalar icin ufak bir kovaryans matrisini verebiliriz, ve o noktalari bir Gaussian altina atarak filtrelenmesini saglayabiliriz.

Not: EM birden fazla Gaussian karisim oranlarini, hem ortalama degerlerini (mu) hem de kovaryans degerlerini veriden hesaplayabilir. Bu ornekte bir kovaryansin bilindigini farz ediyoruz.

Bu algoritmanin neler yaptigina dikkat edelim: Birisi bize bir veri obegi veriyor, biz de diyoruz ki "herhalde bu veri icinde bazi ekstrem noktalar var". Buna gore bir kovaryans ayarliyoruz ve EM'i isletiyoruz, EM pat diye iki Gaussian icin orta noktayi bulup bize veriyor. Kovaryans ta bilmiyor olabilirdik, yine EM isletiyoruz ve bize iki tane Gaussian'i veriden cekip cikartiyor.

Bu onemli bir algoritma. Iteratif olarak isler, detaylar icin literature bakilabilir. Bitis ani log olurluk (likelihood) degerinin artik degismedigi andir.
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats

def gen_data(dim):
sigma_1 = np.eye(dim)*20.0
sigma_2 = np.eye(dim)*2.0

x1 = np.random.multivariate_normal([40.,10.], sigma_1, size=100)
x2 = np.random.multivariate_normal([10.,30.], sigma_2, size=20)
x = np.append(x1, x2, axis=0)

lam_0_1 = np.eye(dim)
lam_0_2 = np.eye(dim)

mu_0_1 = np.array([[1, 1]]).T
mu_0_2 = np.array([[1, 1]]).T
return x.T

def norm_pdf(b,mean,cov):
k = b.shape[0]
part1 = np.exp(-0.5*k*np.log(2*np.pi))
part2 = np.power(np.linalg.det(cov),-0.5)
dev = b-mean
part3 = np.exp(-0.5*np.dot(np.dot(dev.transpose(),np.linalg.inv(cov)),dev))
dmvnorm = part1*part2*part3
return dmvnorm

def EM_known_variance(X, mu, sigma, pi, K, dim, size):
likelihood = -np.Inf
while True:
# E step
for i in range(size):
for k in range(K):
nom = pi[k]*norm_pdf(X[:,i], mu[:,k], sigma[:,:,k])
denom = 0
for j in range(K):
denom += pi[k]*norm_pdf(X[:,i], mu[:,j], sigma[:,:,j])
gamma[k,i] = nom / denom


# M step
N = np.zeros((dim, 1))
for k in range(K):
sum = 0
for i in range(size):
sum += gamma[k, i]
N[k,:] = sum

# new mu
for k in range(K):
sum = 0
for i in range(size):
sum += gamma[k, i]*X[:,i]
mu[:, k] = (1/N[k,:])*sum

# new pi
for k in range(K):
pi[k] = N[k,:] / size

# no estimation of new sigma, we left that one out
# for this example

new_likelihood = 0
for i in range(size):
sum = 0
for k in range(K):
sum += pi[k] * norm_pdf(X[:,i], mu[:,k], sigma[:,:,k])
new_likelihood += np.log(sum)

#print new_likelihood
#print mu[:,0]
#print mu[:,1]

if (np.abs(new_likelihood - likelihood)) < 1e-5: break
likelihood = new_likelihood

print "done"
return mu, pi


K = 2
dim = 2
size = 120
pi = np.ones(dim)*(1.0/K)
X = gen_data(dim)
#mu = np.random.rand(dim, K)*10.0
mu = np.ones((dim, K))*10.0
sigma = np.zeros((dim, dim, K))
gamma = np.zeros((K, size))
sigma[:,:,0] = np.eye(dim)*20.0
sigma[:,:,1] = np.eye(dim)*2.0

mu, pi = EM_known_variance(X, mu, sigma, pi, K, dim, size)
print "EM results"
print "mu", mu
print "pi", pi

No comments: