#!/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"