#!/usr/bin/gnuplot

# standard normal PDF
phi(x) = (1/sqrt(2*3.141592654)) * exp(-.5*x*x)
set nokey
plot phi(x)
pause -1 "normal density\n"

# likelihood for two observations at a and b for a mixture
# model with one center at 0 and one at x (using mutual
# exclusion so one center can't explain both data points)
a = 2
b = -1.5
f(x) = phi(x-a) * phi(0-b) * .5
g(x) = phi(x-b) * phi(0-a) * .5
plot [-3:4] f(x), g(x), f(x)+g(x)
pause -1 "likelihood\n"

# probability of 2 possible data associations as a function
# of moveable center
p(x) = f(x) / (f(x)+g(x))
plot [-3:4] p(x)
pause -1 "probability of data associations\n"

# the EM bound on likelihood, evaluated so that it is
# tangent at x=0
H(p) = p * log(p) + (1-p) * log(1-p)
q(x,p) = p * log(f(x)) + (1-p) * log(g(x)) - H(p)
p0 = p(0)
plot [-3:4] f(x)+g(x), exp(q(x, p0))
pause -1 "EM bound\n"


# bound in log scale -- note the bound is quadratic
plot [-3:4] [-6:-2.5] log(f(x)+g(x)), q(x, p0) 
pause -1 "EM bound, log scale\n"

# bound in probability scale again
plot [-3:4] f(x)+g(x), exp(q(x, p0))
pause -1 "EM bound\n"

# continue the EM iteration: maximizing previous bound
# wrt x gives x = about .3
p0 = p(.3)
plot [-3:4] f(x)+g(x), exp(q(x, p0))
pause -1 "EM bound iteration\n"
p0 = p(1.1)
plot [-3:4] f(x)+g(x), exp(q(x, p0))
pause -1 "EM bound iteration\n"
p0 = p(1.9)
plot [-3:4] f(x)+g(x), exp(q(x, p0))
pause -1 "EM bound iteration\n"



# same likelihood without mutual exclusion
a = 2
b = -1.5
f(x) = phi(x-a) * phi(0-b) * .25
g(x) = phi(x-b) * phi(0-a) * .25
h(x) = phi(x-a) * phi(x-b) * .25
i(x) = phi(0-a) * phi(0-b) * .25
lik(x) = f(x)+g(x)+h(x)+i(x)
plot [-3:4] f(x), g(x), h(x), i(x), lik(x)
pause -1 "likelihood\n"

# assignment probabilities
p1(x) = phi(x-a) / (phi(x-a)+phi(0-a))
p2(x) = phi(x-b) / (phi(x-b)+phi(0-b))
plot [-3:4] p1(x)
pause -1 "data association probability for point a\n"
plot [-3:4] p2(x)
pause -1 "data association probability for point b\n"

# the EM bound on likelihood, evaluated at x0=0
H(p) = p * log(p) + (1-p) * log(1-p)
q(x,p1,p2) = p1 * (1-p2) * log(f(x)) + (1-p1) * p2 * log(g(x)) + p1 * p2 * log(h(x)) + (1-p1) * (1-p2) * log(i(x)) - H(p1) - H(p2)
x0 = 0
plot [-3:4] lik(x), exp(q(x, p1(x0), p2(x0)))
pause -1 "EM bound\n"

# continue the EM iteration: maximizing previous bound
# wrt x gives x = about .3
x0 = .3
plot [-3:4] lik(x), exp(q(x, p1(x0), p2(x0)))
pause -1 "EM bound iteration\n"
x0 = 1.1
plot [-3:4] lik(x), exp(q(x, p1(x0), p2(x0)))
pause -1 "EM bound iteration\n"
x0 = 1.9
plot [-3:4] lik(x), exp(q(x, p1(x0), p2(x0)))
pause -1 "EM bound iteration\n"