Statistical modeling and classification definitions

Last updated 9/8/2003.

General definitions

Annoying warnings

Off[General::spell];
Off[General::spell1];
Off[Part::partd];

Load packages

Required by all modeling techniques

<< Statistics`MultinormalDistribution`;
<< Graphics`Colors`;
<< Graphics`Graphics`;

Histogram modeling

<< Statistics`DataManipulation`;
<< DiscreteMath`ComputationalGeometry`;

Gaussian modeling

<< Graphics`InequalityGraphics`;
<< Graphics`ImplicitPlot`;
<< Graphics`ContourPlot3D`;
<< Graphics`Animation`;

Random variable generation

Uniform random variable

ru[range_] := Random[Real, {-range,  range}];

ru2[range_] := 
        Random[Real, {range[[1]], range[[2]]}];

Two-dimensional uniform random variable

r2[range_] :=  {ru[range],  ru[range]};

Two-dimensional Gaussian random variable

rg[μ_,  Σ_] := 
        Random[MultinormalDistribution[μ,  Σ]];

Bounded two-dimensional Gaussian random variable

rgBounded[{μ_, Σ_}, 
        pr_] := 
Module[{xt, yt,  flag, xmin, xmax, ymin, 
          ymax}, 
{xmin, xmax, ymin, ymax} = 
          Flatten[pr]; 
{xt, yt}  =  
          rg[μ, Σ]; 
flag = -1; 
While[
          flag  <  0, 
If[
            xmin  <  xt  <  xmax  &&  ymin  <  yt  <  ymax,  
            flag  =  1, {xt, yt} = 
              rg[μ,  Σ]]
]; 
{xt, yt}
]

General helper functions

Distance measure between two vectors

dist[a_,  b_] :=  Sqrt[(a - b) . (a - b)];

Convergence function for iterative algorithms

conv[param1_,  
          param2_] := 
Module[{df}, 
df  =  Flatten[param1]  -  
              Flatten[param2]; 
Sqrt[(Plus  @@  (df . 
                      df))/Length[df]]
];

Visualization definitions

Short-cuts

nS =  DisplayFunction -> Identity;
yS  =  
        DisplayFunction -> $DisplayFunction;
imSize  =  ImageSize -> 400;

Two-dimensional data plotting

PlotData2D[x_, pr_, style_:  Red, 
          opt___] := 
Module[{xmin, xmax, ymin, ymax, 
            glist}, 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr] + {- .01,  .01, - .01,  .01}; 
glist  =  
            If[Length[style]  ==  0, {style, Point  /@  x},  
              Join[style, Point  /@  x]]; 
Show[
            Graphics[glist], AspectRatio -> Automatic,  
            PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
            Frame -> True, opt]
];

Histogram modeling definitions

Uniform quantization

UniformQuantization[xrange_,  yrange_,  
          n_] := 
Module[{dx, 
            dy}, 
dx  =  (xrange[[2]] - 
                  xrange[[1]])/(n); 
dy  =  (yrange[[2]] - yrange[[1]])/(n); 
Partition[
              Flatten[Table[{x, y}, {x, xrange[[1]] + dx/2, 
                    xrange[[2]], dx +  .00001}, 
{y,  
                    yrange[[1]] + dy/2, yrange[[2]], dy +  .00001}]], 
              2]  //  N
];

Vector quantization

Classification of data into current centroids

VQClass[
          x_,  μs_] := 
Module[{tmp}, 
tmp = dist[
                  x,  #] &  /@  μs; 
Position[
                tmp,  Min[tmp]]  //  Flatten  //  
            First
];

Splitting of centroids

VQSplitCentroid[μ_, ϵ_:  
            0.0001] := 
{μ  +  ϵ,  μ  -  ϵ};

VQSplitAllCentroids[μs_,  ϵ_:  
            0.0001] := 
Module[{tmp}, 
tmp = 
            VQSplitCentroid  /@  μs; 
Partition[
            Flatten[tmp], Length[First[μs]]]];

One step of VQ algorithm

VQOne[
          data_,  μs_] := 
Module[{k, 
            cl}, 
k = 
            Length[μs]; 
cl  =  VQClass[#,  μs] &  /@  
              data; 
Mean  /@  (Part[data, 
                    Flatten[Position[cl, #]]] &  /@  
                Range[k])
];

One iteration of Linde-Buzo-Grey (LBG) VQ algorithm

VQOneLBG[data_, μs_,  
          tst_:  .01] := 
NestWhile[
          VQOne[data, #] &, μs,  (conv[#1, #2] > tst) &,  
          2];

Full LBG VQ algorithm

VQ[data_, m_:  1, 
          tst_:  .01] := 
Module[{μs}, 
μs  =  {Mean[
                data]}; 
Do[
μs  =  
              VQSplitAllCentroids[μs]; 
μs  =  
              VQOneLBG[data, μs, 
                tst], {m}
]; 
μs
];

Histogram computation

Count number of points in each cluster

CountFrequencies[
          data_,  μs_] := 
Module[{k, cl,  freq,  
            f2}, 
k = 
            Length[μs]; 
cl  =  VQClass[#,  μs] &  /@  data; 
freq  =  
            Frequencies[cl]; 
f2  =  
            Table[0, {Length[μs]}]; 
(f2[[#[[2]]]]  =  #[[1]]) &  /@  
            freq; 
f2/Length[data]  //  
            N
];

Compute histogram matrix

Histogram[data_, μs_, {xrange_, yrange_}, 
          n_:  20] := 
Module[{freq, dx, 
            dy}, 
freq  =  
            CountFrequencies[
              data, μs]; 
dx  =  (xrange[[2]] - 
                  xrange[[1]])/
              n; 
dy  =  (yrange[[2]] - 
                  yrange[[1]])/n; 
Table[
              freq[[VQClass[{x, y}  //  N, μs]]], {x, 
                xrange[[1]],  xrange[[2]], dx}, {y, yrange[[1]], 
                yrange[[2]], dy}]  //  
            Transpose
];

Histogram classification

Compute two-class classification

HistogramClass[hList_] := 
Map[
          If[#[[1]]  <  #[[2]], 0, 1] &, 
          Transpose  /@  Transpose[hList], {2}];

Draw two-class classification

DrawHistogramClass[hList_,  {c1_, c2_}, 
          opts___] := 
Module[{cl, 
            drw}, 
cl = 
            HistogramClass[hList]; 
drw  = 
            Map[If[#  ==  1,  List  @@  c1, List  @@  c2] &,  
              cl, {2}]; 
Show[
            Graphics[Raster[drw,  ColorFunction -> RGBColor]], 
            AspectRatio -> Automatic, 
            opts]
];

Histogram visualization

Visualize output of Histogram[ ]

DrawHistogram[h_, pr_, 
          opts___] := 
Module[{meshRange}, 
meshRange = 
            Partition[Flatten[pr] + {- .01,  .01, - .01,  .01}, 
              2]; 
ListPlot3D[h, 
            MeshRange -> meshRange, PlotRange -> All, 
            Mesh -> False, opts]
];

Visualize partitioning of 2d space

DrawPartitioning[μs_,  pr_, 
          opts___] := 
Module[{vert, val,  
            plotRange}, 
plotRange = 
            Partition[Flatten[pr] + {- .02,  .02, - .02,  .02}, 
              2]; 
{vert, val} = 
            VoronoiDiagram[
              vqCentroids]; 
DiagramPlot[μs, vert, val, 
            LabelPoints -> False,  PlotRange -> pr,  
            Frame -> True, opts]
];

Gaussian modeling definitions

One-dimensional Gaussian modeling

Gaussian PDF

gaussian1D[data_,  x_] := 
        PDF[NormalDistribution[Mean[data], StandardDeviation[data]], 
            x]  //  Simplify;

Log of Gaussian PDF

logGaussian1D[data_, x_] := 
        Log[gaussian1D[data, x]]  //  PowerExpand;

One-dimensional Gaussian

pdf1D[{μ_, σ_},  
          x_:  x] := 
PDF[
            NormalDistribution[μ,  σ],  x]  //  
          Simplify;

Log-likelihood function for two-dimensional Gaussian

logpdf1D[{μ_, σ_},  
          x_:  x] := 
Log[
            PDF[NormalDistribution[μ,  σ],  x]  //  
              Simplify]  //  PowerExpand;

One-dimenisonal Gaussian visualization

One-dimensional  pdf plot

PlotGaussianPDF1D[pdf_, pr_, x_, opts___] := 
 
        Plot[pdf, {x, pr[[1]], pr[[2]]},  PlotRange -> All, 
          Frame -> True, opts];

Plot two one-dimensional pdfs

ShowTwoGaussianPDFs1D[pdf1_,  pdf2_,  pr_,  x_,  
          opts___] := 
Module[{g1, g2, g3, 
            g4}, 
g1 = 
            PlotGaussianPDF1D[pdf1, pr, x, nS,  
              PlotStyle -> {Red,  Thickness[ .01]}, 
              opts]; 
g2 = 
            PlotGaussianPDF1D[pdf2, pr, x, nS,  
              PlotStyle -> {Blue,  Thickness[ .01]}, 
              opts]; 
Show[g1, 
            g2]
];

Two-dimensional Gaussian modeling

Gaussian PDF

gaussian[data_,  {x_,  y_}] := 
        PDF[MultinormalDistribution[Mean[data], CovarianceMatrix[data]], {x, 
              y}]  //  Simplify;

Log of Gaussian PDF

logGaussian[data_, {x_, y_}] := 
        Log[gaussian[data, {x, y}]]  //  PowerExpand;

Equal scalar variances in both dimensions

gaussianScalar[
          data_,  {x_,  
            y_}] := 
Module[{μ, Σ}, 
μ  =  
            Mean[data]; 
Σ  =  (Plus  @@  ((# . #) &  /@  data))/Length[data]*
              IdentityMatrix[2]; 
PDF[
              MultinormalDistribution[μ,  Σ], {x, 
                y}]  //  Simplify];

Log of Gaussian PDF with scalar variance

logGaussianScalar[
          data_, {x_, y_}] := Log[gaussianScalar[data, {x, y}]]  //  
            PowerExpand  //  Chop;

Diagonal covariance matrix

gaussianDiagonal[
          data_,  {x_,  
            y_}] := 
Module[{μ, Σ}, 
μ  =  
            Mean[data]; 
Σ  =  
            DiagonalMatrix[((# . #) &  /@  Transpose[data])/
                Length[data]]; 
PDF[
              MultinormalDistribution[μ,  Σ], {x, 
                y}]  //  Simplify];

Log of Gaussian PDF with diagonal covariance matrix

logGaussianDiagonal[
          data_, {x_, y_}] := Log[gaussianDiagonal[data, {x, y}]]  //  
            PowerExpand  //  Chop;

Log-likelihood function for two-dimensional Gaussian

logpdf[μ_, Σ_,  x_:  x, 
          y_:  y] := 
Log[
            PDF[MultinormalDistribution[μ,  Σ],  {x, 
                  y}]  //  Simplify]  //  PowerExpand;

Two-dimenisonal Gaussian visualization

Three-dimensional pdf plot

PlotGaussianPDF[pdf_, pr_, x_,  y_, opts___] := 
 
        Plot3D[pdf, {x, pr[[1, 1]], pr[[1, 2]]}, {y, pr[[2, 1]], 
            pr[[2, 2]]},  PlotPoints -> 40,  Mesh -> False,  
          PlotRange -> All, opts];

Contour pdf plot

ContourPlotGaussianPDF[pdf_, pr_, x_,  y_, 
          opts___] := 
 
        ContourPlot[
          pdf, {x, pr[[1, 1]], pr[[1, 2]]}, {y, pr[[2, 1]], 
            pr[[2, 2]]},  PlotPoints -> 40,  
          PlotRange -> All, opts];

Plot two classes in graphics array

ShowTwoGaussianPDFs[pdf1_,  pdf2_,  pr_,  x_,  y_,  
          opts___] := 
Module[{g1, g2, g3, 
            g4}, 
g1 = 
            PlotGaussianPDF[pdf1, pr, x, y, Axes -> False, 
              nS]; 
g2 = 
            ContourPlotGaussianPDF[pdf1, pr, x, y, Axes -> False, 
              FrameTicks -> None, nS]; 
g3 = 
            PlotGaussianPDF[pdf2, pr, x, y, Axes -> False, 
              nS]; 
g4 = 
            ContourPlotGaussianPDF[pdf2, pr, x, y, Axes -> False, 
              FrameTicks -> None, nS]; 
Show[
            GraphicsArray[{{g1, g3}, {g2, g4}}], 
            opts]
];

Plot one classification region for two classes

PlotClassificationRegion[pdf1_, pdf2_, pr_, x_, y_, flag_, color_, 
          opts___] := 
Module[{xmin, xmax, ymin, ymax, 
            plotRange, lpdf1, 
            lpdf2}, 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr]; 
plotRange = 
            Partition[Flatten[pr] + {- .01,  .01, - .01,  .01}, 
              2]; 
lpdf1  =  
            Log[pdf1]  //  PowerExpand; 
lpdf2  =  
            Log[pdf2]  //  PowerExpand; 
InequalityPlot[
              flag*lpdf1 > flag*lpdf2,  {x,  xmin, xmax}, {y, ymin, 
                ymax},  PlotRange -> plotRange, opts]  /.  
            Hue[1/2] -> color
];

Plot both classification regions for two classes

PlotTwoClassificationRegions[pdf1_, pdf2_, pr_, x_, y_, colors_, 
          opts___] := 
Module[{g1, 
            g2}, 
g1  =  
            PlotClassificationRegion[pdf1, pdf2, pr, x, y, 1, 
              colors[[1]], nS]; 
g2  =  
            PlotClassificationRegion[pdf1, pdf2, pr, x, y, -1, 
              colors[[2]], nS]; 
Show[g1, g2, 
            opts]
];

Plot classification for two Gaussian pdfs with base graphics overlay

ShowTwoClassificationRegions[pdf1_, pdf2_, pr_, x_, y_,  colors_, 
          gbase_, opts___] := 
Module[{g1}, 
g1 = 
            PlotTwoClassificationRegions[pdf1, pdf2, pr, x, y, colors, 
              nS]; 
If[gbase  ==  False, 
            Show[g1, opts], , 
Show[g1, gbase, 
              opts]]
];

Plot classification for two Gaussian pdfs from Gaussian parameters

ShowDecisionBoundary[μ1_,  Σ1_,  μ2_,  Σ2_,  pr_, cRegion_,  cPDF_, 
          opts___] := 
Module[{lp1,  lp2,  g1, g2, g3, 
            g4, g5, gf, x1, x2, y1, y2, x, y,  
            plotRange}, 
{{x1, x2}, {y1, y2}} = 
            pr; 
plotRange = 
            Partition[Flatten[pr] + {- .01,  .01, - .01,  .01}, 
              2]; 
lp1  = 
            logpdf[μ1, Σ1, x, 
              y]; 
lp2  =  
            logpdf[μ2, Σ2, x, y]; 
g1 =  
            InequalityPlot[lp1 > lp2,  {x,  x1, x2}, {y, y1, y2},  
                nS]  /.  Hue[1/2] -> cRegion[[1]]; 
          g2 =  InequalityPlot[lp1 < lp2,  {x,  x1, x2}, {y, y1, y2},  
                nS]  /.  Hue[1/2] -> 
                cRegion[[2]]; 
g3 = 
            Graphics[{PointSize[ .02], Point[μ1], 
                Point[μ2]}]; 
g4 =  
            ImplicitPlot[
              lp1  ==  (lp1  /.  {x -> μ1[[1]], 
                        y -> μ1[[2]]}) - 0.5,  {x, x1, 
                x2}, {y, y1, y2}, nS,  
              PlotStyle -> {Thickness[ .015], 
                  cPDF[[1]]}]; 
g5 =  
            ImplicitPlot[
              lp2  ==  (lp2  /.  {x -> μ2[[1]], 
                        y -> μ2[[2]]}) - 0.5,  {x, x1, 
                x2}, {y, y1, y2}, nS,  
              PlotStyle -> {Thickness[ .015], 
                  cPDF[[2]]}]; 

gf  =  Show[g1, g2,  g3, g4, g5, yS,  PlotRange -> plotRange,  
              Frame -> True,  Axes -> False, 
              opts]
];

Three-dimensional decision regions

Log-likelihood (3d)

logpdf3D[μ_, Σ_,  x_:  x, y_:  y,  
          z_:  z] := 
Log[
            PDF[MultinormalDistribution[μ,  Σ],  {x, y, 
                  z}]  //  Simplify]  //  PowerExpand;

Three-dimensional decision surface

ShowDecisionSurface[μ1_,  Σ1_,  μ2_,  Σ2_,  
          range_: {{-1.02, 1.02}, {-1.02, 1.02}, {-1.02, 
                1.02}}] := 
Module[{lp1,  lp2,  g1, gf, 
            x1, x2, y1, y2, z1, z2, x, y, 
            z}, 
{{x1, x2}, {y1, y2}, {z1, z2}} = 
            range; 
lp1  = 
            logpdf3D[μ1, Σ1, x, y, 
              z]; 
lp2  = 
            logpdf3D[μ2, Σ2, x, y, 
              z]; 
g1 =  
            ContourPlot3D[
              lp1 - lp2,  {x,  x1, x2}, {y, y1, y2},  {z, z1, z2},  
              DisplayFunction -> Identity]; 
gf  =  
            Show[g1, yS,  PlotRange -> range, ImageSize -> 288,  
              Axes -> False]];

Gaussian mixture modeling (Expectation-Maximization) definitions

EM initialization function

Scalar covariance initialization

\(\(EMInitializationScalar[x_, \ k_] := 
        Module[{n, \ rl, \ μ, \ σ, \ priors, \ sm, \ del2}, 
	
          n\  = \ Length[x]; 
	
          rl\  = \ Table[{Random[Integer, {1, n}]}, \ {k}]; 
	μ\  = \ 
            Partition[Flatten[\(Take[x, \ #] &\)\  /@ \ rl], \ 2]; 
	
          sm\  = \ 
            Table[\(dist[#, \ μ[\([i]\)]] &\)\  /@ \ μ, \ {i, \ 1, \ 
                k}]; 
	
          del2\  = \ 
            Table[{i, i}\ \ , \ {i, \ 1, 
                k}]; 
	σ\  = \ \(\((#*#)\) &\)\  /@ \ \((Min\  /@ \ 
                  Delete[sm, \ del2])\); 
	
          priors\  = \ Table[1/k, \ {k}]\  // \ N; 
	{μ, \ σ, \ 
            priors}
];\)\)

Full covariance initialization

EMInitializationFull[x_, 
          k_] := 
Module[{μ, σ, priors, 
            dim}, 
{μ, σ, priors}  =  
            EMInitializationScalar[x,  k]; 
dim  =  
            Length[First[μ]]; 
{μ, (# IdentityMatrix[dim]) &  /@  σ, 
            priors}
];

Mathematica-native EM mixture modeling functions

Σ conversion routine to full matrix

convertSigma[ {μ_, Σ_}] := 
Module[{type, dim, s}, 
type  =  
            Length[Dimensions[Σ  //  First]] + 
              1; 
dim  =  
            Length[μ[[1]]]; 
s  =  
            Switch[type, 
1,  (# IdentityMatrix[
                        dim]) &  /@  Σ, 
2,  DiagonalMatrix[#] &  /@  Σ, 
3,  Σ]; 
s];

Log likelihood of single point (requires 1d data to be list)

LogLikelihood[
          x_,  μ_, Σ_] := 
Module[{type, 
            dim, s}, 
type  =  
            Length[Dimensions[Σ]] + 
              1; 
dim  =  
            Length[x]; 
s  =  
            Switch[type, 
1,  Σ IdentityMatrix[dim], 
2,  
              DiagonalMatrix[Σ], 
3,  Σ]; 
Log[
            PDF[MultinormalDistribution[μ,  s], 
              x]]
];

Mixture likelihood of single point

MixtureLikelihoodNative[x_, {μ_,  Σ_,  priors_}] := 
        Module[{type, dim, s}, 
type  =  
            Length[Dimensions[Σ  //  First]] + 
              1; 
dim  =  
            Length[x]; 
s  =  
            Switch[type, 
1,  (# IdentityMatrix[
                        dim]) &  /@  Σ, 
2,  DiagonalMatrix[#] &  /@  Σ, 
3,  Σ]; 
Plus  @@ 
            Table[priors[[i]] 
                PDF[MultinormalDistribution[μ[[i]],  s[[i]]], 
                  x], {i, 1, 
                Length[μ]}]
];

Mixture log likelihood of data set

LogLikelihoodMixtureDataSet[
          x_, {μ_,  Σ_,  
            priors_}] := 
Module[{type, dim, s,  f, 
            v}, 
type  =  
            Length[Dimensions[Σ  //  First]] + 
              1; 
dim  =  
            Length[x[[1]]]; 
s  =  
            Switch[type, 
1,  (# IdentityMatrix[
                        dim]) &  /@  Σ, 
2,  DiagonalMatrix[#] &  /@  Σ, 
3,  Σ]; 
          f[v_] = Log[
              Plus  @@ Table[
                  priors[[i]] 
                    PDF[MultinormalDistribution[μ[[i]],  
                        s[[i]]], Table[v[[j]], {j, 1, dim}]], {i, 1, 
                    Length[μ]}]]; 
Plus  @@  (f  /@  x)
];

P((Subscript[ω, i]) | θ, (Subscript[x, j]))

ClassProbabilitiesDataSet[
          x_, {μ_,  Σ_,  
            priors_}] := 
Module[{type, dim, s,  f, 
            v}, 
type  =  
            Length[Dimensions[Σ  //  First]] + 
              1; 
dim  =  
            Length[x[[1]]]; 
s  =  
            Switch[type, 
1,  (# IdentityMatrix[
                        dim]) &  /@  Σ, 
2,  DiagonalMatrix[#] &  /@  Σ, 
3,  Σ]; 
          f[v_] = Table[
              priors[[i]] 
                PDF[MultinormalDistribution[μ[[i]],  s[[i]]], 
                  Table[v[[j]], {j, 1, dim}]], {i, 1, 
                Length[μ]}]; 
((#/(Plus @@ #)) &  /@  (f  /@  x))  //  
            Chop
];

One step EM helper function

diff[x_,  μ_] :=  (# - μ) &  /@  x;

One step of Expectation-Maximization algorithm (scalar, full covariances only)

EMOneStep[
          x_,  {μ_,  Σ_,  
            priors_}] := 
Module[{ptbl, 
            ptmp,  μnew,  Σnew,  priorsnew,  tmp2,  
            type, dim, n}, 
n  =  
            Length[x]; 
dim = 
            Length[First[μ]]; 
type  =  
            Length[Dimensions[Σ  //  First]] + 
              1; 
ptbl  =  
            ClassProbabilitiesDataSet[
              x,  {μ,  Σ,  
                priors}]; 
ptmp =  ((Plus  @@  #) &  /@  Transpose[ptbl]); 
priorsnew  =  
            ptmp/n; 
μnew  =  (Transpose[ptbl] . 
                  x)/ptmp; 
Switch[
            type, 
1, 
tmp2  =  
              Map[(# . #) &, (diff[
                          x,  #] &  /@  μnew)  //  
                  Transpose, {2}]; Σnew  =  ((Plus  @@  #) &  /@  ((ptbl tmp2)  //  Transpose))/ptmp/
                dim, 
3, 
tmp2  =  
              Map[Outer[
                    Times, #, #] &, (diff[
                          x,  #] &  /@  μnew)  //  
                  Transpose, {2}]; Σnew  =  ((Plus  @@  #) &  /@  ((ptbl tmp2)  //  Transpose))/
                ptmp]; 
{μnew, Σnew, 
            priorsnew}
];

Full EM algorithm

EM[x_,  {μ_,  Σ_,  priors_},  
          tst_] := 
NestWhileList[
          EMOneStep[x, #] &,  {μ,  Σ,  
            priors},  (conv[#1, #2] > tst) &,  2];

Mixture likelihood over a range of inputs

MixtureLikelihoodMatrix[range_,  
          n_, {μ_,  Σ_,  
            priors_}] := 
Module[{type, dim, s,  f, v, xd, 
            xmin, xmax, ymin, ymax, dx, 
            dy}, 
{xmin, xmax, ymin, ymax}  =  
            Flatten[range]; 
dx  =  (xmax  -  xmin)/
              n; 
dy  =  (ymax  -  ymin)/
              n; 
type  =  
            Length[Dimensions[Σ  //  First]] + 
              1; 
dim  =  
            Length[x[[1]]]; 
s  =  
            Switch[type, 
1,  (# IdentityMatrix[
                        dim]) &  /@  Σ, 
2,  DiagonalMatrix[#] &  /@  Σ, 
3,  Σ]; 
          f[v_] = Plus  @@ 
              Table[priors[[i]] 
                  PDF[MultinormalDistribution[μ[[i]],  s[[i]]], 
                    Table[v[[j]], {j, 1, dim}]], {i, 1, Length[μ]}]; 
          xd =  Table[{x, y}, {x, xmin, xmax + dx/2, dx}, {y, ymin, 
                ymax + dy/2, dy}]; 
Map[f, xd, {2}]  //  
            Transpose];

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

RandomMixtureModelOne[{μ_,  Σ_,  priors_}] := 
        Module[{type, dim, s, class}, 
type  =  
            Length[Dimensions[Σ]] + 
              1; 
dim  =  
            Length[First[μ]]; 
s  =  
            Switch[type, 
1,  Σ IdentityMatrix[dim], 
2,  
              DiagonalMatrix[Σ], 
3,  Σ]; 
class = 
            PickOne[priors]; 
Random[
            MultinormalDistribution[μ[[class]],  
              s[[class]]]]
];

RandomMixtureModel[{μ_,  Σ_,  priors_},  
          n_] := 
Table[
          RandomMixtureModelOne[{μ, Σ, 
              priors}], {n}];

Two-dimensional Gaussian mixture model visualization

Mixture means without overlay in 2d

PlotMeans[{μ_, Σ_, priors_},  pr_,  
          color_:  Black,  
          opt___] := 
Module[{lp, s, x, y, xmin, xmax, 
            ymin, ymax, g1, 
            g2}, 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr]; 
s  =   convertSigma[{μ, Σ}]; 
g1 = 
            Table[
lp  = 
                logpdf[μ[[i]], s[[i]], x, y]; 
              ImplicitPlot[
                lp  ==  (lp  /.  {x -> μ[[i, 1]], 
                          y -> μ[[i, 2]]}) - 0.5,  {x, xmin, 
                  xmax}, {y, ymin, ymax}, nS,  
                PlotStyle -> {Thickness[ .005], color},  
                PlotRange ->  All], {i, 1, 
                Length[μ]}]; 
g2  =  
            Show[Graphics[{PointSize[ .03], Point  /@  μ}], 
              nS]; 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr] + {- .01,  .01, - .01,  .01}; 
Show[g1, g2,  
            PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
            Frame -> True, AspectRatio -> Automatic, 
            opt]
];

Mixture model contour plot (with shading)

PlotContourPDF[{μ_, Σ_, priors_},  pr_, n_:  50,  
          opt___] := 
Module[{xmin, xmax, ymin, ymax, 
            g1}, 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr]  //  N; 
ListContourPlot[
            MixtureLikelihoodMatrix[{xmin, xmax, ymin, ymax}, 
              n,  {μ,  Σ,  priors}], 
            MeshRange -> {{xmin, xmax}, {ymin, ymax}},  PlotRange -> All, 
            opt]
];

Mixture model contour plot (no shading)

PlotContourPDFOverlay[{μ_, Σ_, priors_},  pr_, 
          n_:  50,  
          opt___] := 
Module[{xmin, xmax, ymin, ymax, 
            g1}, 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr]  //  N; 
ListContourPlot[
            MixtureLikelihoodMatrix[{xmin, xmax, ymin, ymax}, 
              n,  {μ,  Σ,  priors}],  
            ContourShading -> False, 
            MeshRange -> {{xmin, xmax}, {ymin, ymax}},  PlotRange -> All, 
            opt]
];

Mixture model 3d surface plot

PlotPDF[{μ_, Σ_, priors_},  pr_, n_:  50,  
          opt___] := 
Module[{xmin, xmax, ymin, ymax, 
            g1}, 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr]  //  N; 
ListPlot3D[
            MixtureLikelihoodMatrix[{xmin, xmax, ymin, ymax}, 
              n,  {μ,  Σ,  priors}], 
            MeshRange -> {{xmin, xmax}, {ymin, ymax}},  PlotRange -> All, 
            Mesh -> False, opt]
];

Log-likelihood trajectory

PlotLogLikelihood[em_,  xd_,  
          opt___] := 
Module[{logp,  traj,  g1, 
            g2}, 
logp  =  LogLikelihoodMixtureDataSet[
                  xd, #] &  /@  em; 
traj  = 
            Transpose[{Table[i,  {i, 0, Length[logp] - 1}], 
                logp}]; 
g1  =  
            ListPlot[traj,   Axes -> False,  Frame -> True, 
              PlotRange -> All, PlotJoined -> True, 
              DisplayFunction  ->  Identity]; 
          g2  =  Show[
              Graphics[{Red, PointSize[ .02],  
                      Point[#]} &  /@   traj],  
              DisplayFunction  ->  Identity,  Frame -> True, 
              PlotRange -> All, opt]; 
Show[ g1,  
            g2,  yS,  AspectRatio  -> 1/GoldenRatio,  
            GridLines -> Automatic,  opt]
];

Log-likelihood trajectory (varying number of classes)

PlotLogLikelihood2[em_,   xd_, 
          opt___] := 
Module[{logp,  traj,  g1, g2,  
            nc}, 
nc  =  Length[First[#]] & /@  
              em; 
logp  =  LogLikelihoodMixtureDataSet[
                  xd, #] &  /@  em; 
traj  = 
            Transpose[{nc, logp}]; 
g1  =  
            ListPlot[traj,   Axes -> False,  Frame -> True, 
              PlotRange -> All, PlotJoined -> True, 
              DisplayFunction  ->  Identity]; 
          g2  =  Show[
              Graphics[{Red, PointSize[ .02],  
                      Point[#]} &  /@   traj],  
              DisplayFunction  ->  Identity,  Frame -> True, 
              PlotRange -> All, opt]; 
Show[ g1,  
            g2,  yS,  AspectRatio  -> 1/GoldenRatio,  
            GridLines -> Automatic,  opt]
];

Plot multiple results for mixture model

EMAllPlot[{μ_, Σ_, priors_},  gin_, 
          pr_: {{0, 1}, {0, 1}}] := 
Module[{g1, g2, g4, 
            g5}, 
g2  =  
            PlotMeans[{μ, Σ, priors},  pr,  Blue,  
              nS]; 
g1 = 
            Show[gin, g2,  nS]; 
g4  =  
            PlotContourPDF[{μ, Σ, priors},  pr, 50, nS,  
              FrameTicks -> None]; 
g5  =  
            PlotPDF[{μ, Σ, priors},  pr, 50, nS,  
              Axes -> False]; 
Show[
            GraphicsArray[{{g1, gin}, {g4, g5}}], yS,  
            ImageSize -> 400]
];