Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Expectation-Maximization Implementation based on the book "Machine Learning" by Tom M. Mitchell
# - Expectation-Maximization Implementation based on the book "Machine Learning" by Tom M. Mitchell
# - Find the mean of One Gaussian; and of a mixture of Two Gaussians
# Copyright (C) 2015 Eric Aislan Antonelo
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
using Gadfly
using Distributions
################## Uma Gaussiana
mu = 5
sigma = 10
n = 1000
srand(123) # Setting the seed
d = Normal(mu,sigma)
points = rand(d, n)
plot(x=1:length(points), y=points)
#### a Estimativa é a média dos pontos
# np.mean(points)
# plt.show()
size(points)
################## Duas Gaussianas
mus = [-10, 30]
sigma = 9
n = 1000
distr = [Normal(mu,sigma) for mu in mus]
x = []
n_each = div(n, length(distr))
for d in distr
x = [x; rand(d, n_each)]
end
x = x[randperm(n)]
#x = np.array([sample(sel_mu) for i in range(n)])
plot(x=1:n, y=x)
##################
#### Usaremos EM (Expectation Maximiation) para estimar as médias (hipóteses) de cada Gaussiana a partir dos dados
function EM_Gaussians(x, n_h, sigma)
## initial guess for hypotheses (mean of Gaussians)
n = length(x)
h = [sample(x) for i in 1:n_h]
h_ = deepcopy(h)
println(h)
println(h_)
println("Estimativa inicial, medias: $h")
# initialize vector for non-observable variables z_i1, z_i2 (which indicate which Gaussian generated the i_th point)
z = zeros(n, length(h))
diff = 10
_inc = 0.00001 # avoid division by zero
iterations = 0
## anonymous functions:
p = ((x,u) -> exp((-1/(2*sigma^2))*(x-u)^2))
expected_z = ((i,j) -> p(x[i],h[j]) / (sum([ p(x[i],h[n]) for n in 1:n_h]) + _inc))
while diff > 0.001
## E-step : estimation step
#- z will hold the expected estimation for the non-observable variables
for i in 1:n
for j in 1:n_h
z[i,j] = expected_z(i, j)
#print i,j,x[i], h[j], z[i][j], sum([ p(x[i],h[n]) for n in range(h.size) ])
end
end
## M-step: find new hypotheses (means of Gaussians) which maximize the likelihood of observing the full data y = [x,z], given the current h
## Soma ponderada dos pontos, peso = probabilidade do ponto pertencer a tal Gaussiana (do passo anterior)
for j in 1:n_h
h_[j] = sum([ z[i,j] * x[i] for i in 1:n ]) / (sum([ z[i,j] for i in 1:n ]) + _inc)
end
println(h_)
# store in diff value to check for loop termination
diff = mean(map(abs, [ h[j] - h_[j] for j in 1:n_h ] ))
# diff = mean(abs([ h[j] - h_[j]) for j in 1:n_h ]))
iterations += 1
#print h, h_
#print diff, iterations
# update current hypothesis
h = deepcopy(h_)
end
(h,z)
end
##################
medias, z = EM_Gaussians(x,length(mus),sigma)
println("Medias das Gaussianas, resultado: $medias")
println("Medias Reais: $mus")
##################
## Plota estimativa final para probabilidades de cada ponto pertencer a qual Gaussiana
j = 1
#plot(z[:,j],'*-')
x_1 = zeros(n)
x_2 = zeros(n)
x_1[:] = NaN
x_2[:] = NaN
for i in 1:n
if z[i,1] > z[i,2]
x_1[i] = x[i]
else
x_2[i] = x[i]
end
end
# preto: medias reais, laranja: medias estimadas por EM
plot(layer(x=1:n, y=[x_1 x_2], Geom.point),
layer(yintercept=medias,Geom.hline(color="orange", size=0.8mm)),
layer(yintercept=mus,Geom.hline(color="black", size=1.5mm)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.