3
$\begingroup$

I would like to use the expectation maximization algorithm within Mathematica. The issue is that I am very basic user and progress slowly with understanding functional programming.. I have found this "featured example" on Wolfram web page, where an algorithm for a mixture of 2 normal distributions is demonstrated. Now I would like to understand what is going on in order to change it later to be applicable to my problem (My problem is to estimate a mixture of 3 distributions, given the real data) Could please anyone explain line by line (or more coarsely) what is going on in this code? (And how would one go to use the real dataset (excel file)) Thanks a lot!

{Manipulate[
 Module[{d1, d2, sim, updates},
  SeedRandom[sr];
  d1 = NormalDistribution[μ1, σ1];
  d2 = NormalDistribution[μ2, σ2];
  sim = With[{musAndsigmas = 
      RandomChoice[{p1, 
         1 - p1} -> {{μ1, σ1}, {μ2, σ2}}, 
       numsamples]},
    musAndsigmas[[All, 2]] RandomReal[NormalDistribution[], 
       numsamples] + musAndsigmas[[All, 1]]];
  updates = 
   NestList[em[#, sim] &, {{-10, 5}, {5, 5}, {.2, .8}}, numiterations];
  Grid[{{Show[{Histogram[sim, Automatic, "PDF", 
        PlotRange -> {{-30, 30}, {0, Automatic}}], 
       Plot[{updates[[-1, 3, 1]]*
          PDF[Apply[NormalDistribution, updates[[-1, 1]]], x], 
         updates[[-1, 3, 2]]*
          PDF[Apply[NormalDistribution, updates[[-1, 2]]], 
           x]}, {x, -30, 30}, 
        PlotStyle -> {{Red, Thick}, {Blue, Thick}}, 
        PerformanceGoal -> "Speed", PlotRange -> All]}, 
      ImageSize -> {400, 275}, 
      AxesOrigin -> {-30, 
        0}]}, {Grid[{{Panel[
         TableForm[updates[[-1, 1 ;; 2]], 

          TableHeadings -> {updates[[-1, 3]], {"μ", 
             "σ"}}]]}}]}}]
  ],
 Style["first distribution", Bold],
 {{μ1, 8}, -20, 20, ImageSize -> Tiny, Appearance -> "Labeled"},
 {{σ1, 5}, 1, 10, 1, ImageSize -> Tiny, 
  Appearance -> "Labeled"},
 Delimiter,
 Style["second distribution", Bold],
 {{μ2, -12}, -20, 20, ImageSize -> Tiny, Appearance -> "Labeled"},
 {{σ2, 3}, 1, 10, 1, ImageSize -> Tiny, 
  Appearance -> "Labeled"},
 Delimiter,
 Style["probability", Bold],
 {{p1, .3, "p"}, .1, .9, ImageSize -> Tiny, Appearance -> "Labeled"},
 Delimiter,
 Style["number of samples", Bold],
 {{numsamples, 100, ""}, 50, 200, 10, ImageSize -> Tiny, 
  Appearance -> "Labeled"},
 Style["number of iterations\nof EM algorithm", Bold],
 {{numiterations, 3, ""}, 1, 20, 1, ImageSize -> Tiny, 
  Appearance -> "Labeled"},
 Delimiter,
 {{sr, 326, "random seed"}, 1, 9999, 1, ImageSize -> Tiny, 
  Appearance -> "Labeled"},
 ControlPlacement -> Left,
 Initialization :> (
   problists[θ_, y_] := Block[{probs, totalprobs},
     probs = 
      Table[θ[[3, i]]*
        Map[PDF[Apply[NormalDistribution, θ[[i]]], #] &, 
         y], {i, 2}];
     totalprobs = Total[probs];
     Map[#/totalprobs &, probs]];
   pi[j_, p_] := Mean[p[[j]]];
   emu[sim_, j_, p_] := Total[sim*p[[j]]]/Total[p[[j]]];
   emstd[sim_, j_, u_List, p_] := Total[(sim - u[[j]])^2*p[[j]]]/
    Total[p[[j]]];
   em[params_, sim_] := Module[{theprobs, mus, vars, probs},
     theprobs = problists[params, sim];
     mus = Map[emu[sim, #, theprobs] &, {1, 2}];
     vars = 
      Map[emstd[sim, #, params[[1 ;; 2, 1]], theprobs] &, {1, 2}];
     probs = Map[pi[#, theprobs] &, {1, 2}];
     Append[Transpose[{mus, vars^.5}], probs]];
   )]}
$\endgroup$
  • $\begingroup$ For reference, this is the Wolfram Demo Project you mean, I presume? demonstrations.wolfram.com/… $\endgroup$ – dr.blochwave Sep 5 '14 at 13:41
  • $\begingroup$ yes, exactly this one.... I'd really like to understand how all is done. $\endgroup$ – french cactus Sep 5 '14 at 13:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.