(* EXAMPLE OF LIKELIHOOD CALCULATION FOR A MIXTURE DISTRIBUTION *)

(* Probability Density Function for univariate normal with mean mu and variance var *)
NormalPdf[x_,mu_,var_] := 1/Sqrt[2*Pi*var]*Exp[-(x-mu)^2/(2*var)]
Plot[NormalPdf[x,0,1],{x,-3,3}];

(* Likelihood of data x arising from a mixture of two normals with means \
m[[1]], m[[2]], both with unit variance, and equal mixing proportions *)
Likelihood[x_,m_] := 
  Product[Sum[NormalPdf[x[[i]],m[[j]],1]/2, {j,1,2}],{i,1,Length[x]}]

(* Plot of likelihood for 9 points (x) in 3 clusters.  
   Small power of Likelihood[...] just to scale it for better visualization. *)
z=List[-.2,0,.2];
x=Join [z-10,z,z+12]
margin = .5*(Max[x]-Min[x]);
Plot3D[Likelihood[x,{m1,m2}]^.02,
    {m1,Min[x]-margin,Max[x]+margin},
    {m2,Min[x]-margin,Max[x]+margin}, 
    PlotPoints\[Rule]50,
    PlotRange\[Rule]{0,.17}];

(* Save plot as .pdf *)
Display["~/Desktop/likelihood.pdf",%,"PDF"];

(* Find 4 modes: *)
Maximize[{Likelihood[x,{m1,m2}]^.02,m1<m2 && m2>11},{m1,m2}]
Maximize[{Likelihood[x,{m1,m2}]^.02,m1<m2 && m2<08},{m1,m2}]
Maximize[{Likelihood[x,{m1,m2}]^.02,m1>m2 && m1<06},{m1,m2}]
Maximize[{Likelihood[x,{m1,m2}]^.02,m1>m2 && m1>11},{m1,m2}]