Needs["Statistics`Master`"]

EM Algorithm (Expectation - Maximization)

McLachlan, G. J. and Krishnan, T.: The EM Algorithm and Extensions. Wiley, New York 1997

The principle: There are some "missing" data; solve the likelihood equation indirectly, iteratively.

Y the random vector corresponding to the observed data y having density g(y,Ψ), Ψ  vector of unknown parameters. The observed data y are viewed as being incomplete and is regarded as an observable function of the so called complete-data. Let x  denote the vector containing the augement or so called complete-data, and let z denote the vector containing the additional data, reffered to as as the  unobservable or missing data.
Let (Subscript[g, c])(x;Ψ) denote the pdf of X corresponding to the complete-data vector x -> complete-data log likelihood function that could be formed for Ψ  if  x  were fully observable is
	
	log (Subscript[L, c])(Ψ) = log (Subscript[g, c])(x;Ψ) .
	
Formally, we have two sample spaces, X and Y, and a many -to-one mapping from X to Y. Instead of observing the complete-data vector x in X, we observe the incomplete-data vector y=y(x) in Y ->

	g(y,Ψ) =(∫_(y = y(x)))  (Subscript[g, c])(x; Ψ)dx	



A crucial point: To find a formulation of the estimation problem with missing data or, alternatively, to find a proper space X.

2 steps of EM algorithm:

given the estimates from the k-th iteration (the initial estimates should be chosen carefully)
E - step:
Calculate  Q(Ψ;(Ψ^((k))))=(Subscript[E, (Ψ^((k)))]){log (Subscript[L, c])(Ψ) | y}
M -step:
Choose (Ψ^((k + 1))) = arg max {Q(Ψ;(Ψ^((k)))), Ψ  from the corresponding parametric space}

Repeate the E- and M- steps until  |(Ψ^((k + 1))) - (Ψ^((k)))|  becomes sufficiently small.

Example (Randomly censored exponentially distributed data)
W ~ Exp(μ), y = ((((Subscript[y, 1]^T), ...,(Subscript[y, n]^T)))^T), where (Subscript[y, j]) = ((((Subscript[c, j]),(Subscript[δ, j])))^T), (Subscript[δ, j]) = 0 or 1 according as the observation is censored or uncensored at time (Subscript[c, j]).

McDonald, J. B.: Probability Distributions for Financial Models, Chapter 14 In: Rao, C. R., Maddala, G. S. (eds.): Statistical Methods in Finance, Handbook of Statistics 14, Elsevier, Amsterdam 1996 (useless)

fr[x_, μ_, σ_] := 
      Exp[-(x - μ)^2/σ^2]
mixture1[x_, p1_, p2_, 
        p3_, μ1_, μ2_, μ3_, σ_] := 
      p1*fr[x, μ1, σ] + p2*fr[x, μ2, σ] + 
        p3*fr[x, μ3, σ]

mixture1[x, p1, p2, 
      p3, μ1, μ2, μ3, σ]

Plot[
        Table[mixture1[x, 0.2, 0.6, 0.2, -j/2. , 0. , (5 - j)/2. , 
              1.0], {j, 0, 5, 2}] // Evaluate, {x, -3, 3}];

{m  =  0,  s  =  1,  eps  =  3. ,  p1  =  0.2,  
        p3  =  0.2};
\( (*\ Number\ of\ observations\ *) \ \(pocet\  = \ 
        100;\)
 (*\(n\  = \ pocet;\)*) \)
p2  =  1  -  p1  -  p3;
m1  =  m  -  eps;
m2  =  m;
m3  =  m  +  eps;
pole1  =  
        RandomArray[NormalDistribution[m1,  s],  
          IntegerPart[p1*pocet]];
pole2  =  
        RandomArray[NormalDistribution[m2,  s],  
          IntegerPart[p2*pocet]];
pole3  =  
        RandomArray[NormalDistribution[m3,  s],  
          IntegerPart[p3*pocet]];
xx  =  Join[pole1,  pole2,  pole3];
{Mean[xx], Variance[xx], s + 2  p1*eps^2}
\(mixture[x_] := 
      p1*fr[x, m1, s] + 
        p2*fr[x, m2, s] + \(\(p3*fr[x, m3, s]\)\(
\)\(		\)\)\)

Plot[mixture[x], {x, -5, 5}];

\(
Do[n = Length[xx]; 

      p = Map[Function[x, 
              Map[fr[x, #, σ] &, μ]/
                Apply[Plus, α*Map[fr[x, #, σ] &, μ]]], xx] // 
          N; 
p = Map[\((#*α)\) &, p]; 
α = 
        Apply[Plus, p]; 
μ = 
        xx . p/α; 
α = α/n; 
σ = 
        Sqrt[Apply[Plus, 
              Apply[Plus, p*Transpose[Map[\((xx - #)\)^2 &, μ]]]]/n]; 

      Print["μ= ", μ]; 
Print[	"α= ", α]; 

      Print[	"σ2= ", σ^2], {20}]\)
\(Print["ϵ=  ", ϵ = μ[\([2]\)] - μ[\([1]\)], "   ", μ[\([3]\)] - μ[\([2]\)]]\)
\(Print["Estimated variance= ", σ^2 + α . μ^2]
\)