Two-parameter mixture-model example

This notebook explores a two-parameter mixture-model example and assumes that the definitions in the notebook "statistical_modeling.nb" have already been evaluated.

Local definitions

Random numbers from Gaussian mixture model

Uniform random number in range [0, 1]

RR  :=  Random[Real,  {0. ,  1. }];

Random number generation helper functions

TransformProbabilityVector[
          prob_] := 
Module[{s  =  0,  
            len}, 
len  =  
            Length[prob] - 1; 
Join[
            Table[s  =  s  +  prob[[i]],  {i,  1,  
                len}],  {1}]
];

PickOne[ 
          priors_] := 
Module[{r}, 
r  =  RR; (Position[(r  <=  #) &  /@  
                    TransformProbabilityVector[priors], True]  //  
                Flatten)  //  First
];

Random number generation from Gaussian mixture model (1D)

RandomMixtureModelOne1D[{μ_,  σ_,  priors_}] := 
        Module[{class}, 
class = 
            PickOne[priors]; 
{class, 
            Random[NormalDistribution[μ[[class]],  σ[[class]]]]}
];

RandomMixtureModel1D[{μ_,  σ_,  priors_},  
          n_] := 
Table[
          RandomMixtureModelOne1D[{μ, σ, 
              priors}], {n}];

1D Gaussian mixture modeling

1D Mixture model PDF

mixturepdf1D[{μ_,  σ_,  priors_}, 
          x_] := 
Plus  @@ 
          Table[PDF[NormalDistribution[μ[[i]], σ[[i]]], 
                x] priors[[i]], {i, 1, Length[μ]}];

1D log mixture model PDF

mixturelogpdf1D[{μ_,  σ_,  priors_}, 
          x_] := 
Log[
          mixturepdf1D[{μ, σ, priors}, x]];

Find best estimate for μ1 and μ2

likelihoodFunction1D[data_, 
          params_] := 
Plus  @@  (mixturelogpdf1D[
                  params  //  N, #[[2]]] &  /@  
              data);

estimateParams1D[
          data_,  {μ_, σ_, 
            priors_}] := 
Module[{mlf, μ1, μ2,  
            params,  
            sol}, 
params  =  {{μ1, μ2}, σ, priors}; 
mlf  =  
            likelihoodFunction1D[data, params]; 
sol  =  
            FindMinimum[-mlf, {μ1,  μ[[1]]}, {μ2, μ[[2]]}]; 
{{{μ1, μ2}  /.  
                sol[[2]] , σ, 
              priors}, -sol[[1]]}
];

Special plotting functions

Plot 1D labeled data

PlotDataLabeled1D[data_,  c1_,  c2_,  
          opts___] := 
Module[{colors}, 
colors  =  {c1, c2}; 
Show[
            Graphics[{colors[[#[[1]]]], PointSize[ .02], 
                    Point[{#[[2]], 1}]} &  /@  data], 
            Axes -> {True, False}, Frame -> False, 
            AspectRatio -> Automatic,  PlotRange -> All,  
            opts]
];

Plot log-likelihood as a function of μ1 and μ2

PlotLogLikelihoodLocal[
          data_, {{μ1_, μ2_}, σ_, priors_}, 
          opts___] := 
Module[{lf}, 
lf  =  Plus  @@  (mixturelogpdf1D[{{μ1, μ2}, σ, 
                          priors}  //  N, #[[2]]] &  /@  
                  data); 
ContourPlot[
            lf, {μ1, -6, 6}, {μ2, -6, 6},  
            PlotPoints -> 30, opts]
];

PlotLogLikelihoodLocal2[
          data_, {{μ1_, μ2_}, σ_, priors_}, 
          opts___] := 
Module[{g1, g2, est1, est2, c1, 
            c2}, 
g1  =  
            PlotLogLikelihoodLocal[
              data, {{μ1, μ2}, σ, priors}, 
              nS]; 
est1  =  
            estimateParams1D[
              data,  {{-6, 6}, σ, 
                priors}]; 
est2  =  
            estimateParams1D[
              data,  {{6, -6}, σ, priors}]; 
If[
            Last[est1]  <  Last[est2], c1 = Blue; c2 = Red, c1 = Red; 
            c2 = Blue]; 
g2 = 
            Show[Graphics[{{c1, PointSize[ .02], 
                    Point[First[est1][[1]]]}, 
{c2, 
                    PointSize[ .02], Point[First[est2][[1]]]}}], 
              Frame -> True, AspectRatio -> Automatic,  
              nS]; 
Show[g1, g2, yS,  
            PlotRange -> All,  
            opts]; 
{First[est1][[1]], First[
                est2][[1]]}
];

PlotLogLikelihoodLocal3[
          data_, {{μ1_, μ2_}, σ_, priors_}, 
          opts___] := 
Module[{lf}, 
lf  =  Plus  @@  (mixturelogpdf1D[{{μ1, μ2}, σ, 
                          priors}  //  N, #[[2]]] &  /@  
                  data); 
Plot3D[
            lf, {μ1, -6, 6}, {μ2, -6, 6},  
            PlotPoints -> 30, opts]
];

Plot optimization algorithm trajectory

PlotTrajectory[pts_,  color_,  
          opts___] := 
Module[{g1, 
            g2}, 
g1  =  
            Graphics[{PointSize[ .015], color, 
                Point  /@  pts}]; 
g2  =  
            Graphics[{color, Line[pts]}]; 
Show[g1, g2, 
            opts]];

Mixture of two Gaussians (unknown (Subscript[μ, 1])and (Subscript[μ, 2]))

Mixture model definitions

σ  =  {1, 1};
Pω  =  {1/3, 2/3};
μ  =  {-2, 2};
μ1 =. ; μ2 =. ;
μv  =  {μ1,  μ2};

params1  =  {μ,  σ, 
          Pω};
params2  =  {μv, σ, Pω};

Visualize mixture model

pdf  =  
        Plot[mixturepdf1D[params1, x], {x, -5, 5},  
          Frame -> True,  Axes -> False,  
          PlotStyle -> {Thickness[ .01], Green}];

Generate labeled data from mixture model

n  =  100;
data  =  RandomMixtureModel1D[params1, n];

Visualize generated data

PlotDataLabeled1D[data, Red,  Blue];

Plot log likelihood of data as a function of μ1 and μ2

{est1, est2} = PlotLogLikelihoodLocal2[data, params2]

PlotLogLikelihoodLocal3[data, params2];

likelihoodFunction1D[data,  
      Join[{est1},  Drop[params2, 1]]]

likelihoodFunction1D[data,  
      Join[{est2},  Drop[params2, 1]]]

Visualize resulting pdfs

g1  =  
        Plot[mixturepdf1D[{est1, σ, Pω}, x], {x, -5, 5},  
          Frame -> True,  Axes -> False,  
          PlotStyle -> {Thickness[ .01], Red},  
          Evaluate[nS]];
g2  =   Plot[
          mixturepdf1D[{est2, σ, Pω}, x], {x, -5, 5},  
          Frame -> True,  Axes -> False,  
          PlotStyle -> {Thickness[ .01], Blue},  
          Evaluate[nS]];
Show[pdf, g1, g2,  yS];