Mixture-of-Gaussians example #2

This notebook illustrates Gaussian mixture modeling for some 2d synthetic data, and assumes that the definitions in the notebook "statistical_modeling.nb" have already been evaluated.

Global user definitions

Colors representing two classes

Data colors

data1Color = Red;
data2Color = Blue;

style1  =  {data1Color, 
          PointSize[ .02]};
style2  =  {data2Color, PointSize[ .02]};

Classification region colors

class1Color  =  RGBColor[1, 1,  .5];
class2Color  =  RGBColor[ .5, 1, 1];

classColors  =  {class1Color,  class2Color};

Default plot range

pr  =  {{0, 1}, {0, 1}};

No frame ticks

nF = FrameTicks -> None;
FS = Sequence[yS,  nF];

Local definitions

Special data generation functions

\(\(withinRectangle[p_, \ r_] := 
	
        Module[{rx1, \ ry1, \ rx2, \ ry2, \ x, 
            y}, 
		{rx1, \ ry1, \ rx2, \ ry2}\  = \ r; 
		{x, y}\  = \ 
            p; 
		If[
            rx1\  < \ x\  < \ rx2\  && \ 
              ry1\  < \ y\  < \ ry2, 
			True, \ False]];\)\)

\(withinRegion2[p_] := 
	
      If[
			\ withinRectangle[p, \ {0.2, \ 0.6, \ 0.8, \ 0.8}]\  || \ 
          withinRectangle[p, \ {0.2, \ 0.2, \ 0.8, \ 0.4}]\  || \ 
          withinRectangle[p, \ {0.8, \ 0.0, \ 1.0, \ 1.0}], 
		True, \ 
        False]\)

\(\(region1Data[n_] := 
	
        Module[{count, \ xd, p}, 
		count\  = \ 1; 
		xd\  = \ 
            Table[0, \ {n}]; 
		While[
            count\  <= \ 
              n, 
			p\  = \ {Random[Real, \ {0. , 1. }], \ 
                Random[Real, \ {0. , 1. }]}; 
			If[
              withinRegion2[p] == 
                False, 
				xd[\([\(count++\)]\)]\  = \ p]]; 
	
          xd];\)\)

\(\(region2Data[
          n_] := 
		Module[{count, \ xd, p}, 
		count\  = \ 
            1; 
		xd\  = \ Table[0, \ {n}]; 
		While[
            count\  <= \ 
              n, 
			p\  = \ {Random[Real, \ {0. , 1. }], \ 
                Random[Real, \ {0. , 1. }]}; 
			If[
              withinRegion2[p] == True, 
				xd[\([\(count++\)]\)]\  = \ 
                p]]; 
	xd];\)\)

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

\(\(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];\)\)

\(\(GenerateFromMixturePDF[{mu_, \ sigma_, \ priors_}] := 
	
        Module[{k, \ d, \ class}, 
		k\  = \ Length[priors]; 
		d\  = \ 
            Length[mu[\([1]\)]]; 
	\ \ class\  = \ 
            PickOne[priors]; 
	\ \ Random[
            MultinormalDistribution[mu[\([class]\)], \ 
              sigma[\([class]\)]\ IdentityMatrix[d]]]\ 
		];\)\)

specialClassPlot  := 
        Show[Graphics[{{class1Color, 
                Rectangle[{0, 0}, {1, 1}]},  {class2Color,  
                Rectangle[{0.8, 0}, {1, 1}],  
                Rectangle[{0.2, 0.6}, {0.8, 0.8}],  
                Rectangle[{0.2, 0.2},  {0.8, 0.4}]}}],  
          AspectRatio -> Automatic,  Frame -> True,  nF];

Visualize interlocking regions in 2d

Saved data definition

poly  = {Partition[ {0, 0, 8, 0, 8, 2, 2, 2, 2, 4, 8, 4, 8, 6, 2, 
                6, 2, 8, 8, 8, 8, 10, 0, 10}, 2]/10. };

Draw polygonal region

drawPolygonalRegion[polyL_,  pr_,  color_:  LightPink, 
          opt___] := 
Module[{xmin, xmax, ymin, 
            ymax}, 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr] + {- .01,  .01, - .01,  .01}; 
Show[Graphics[{color, Polygon  /@  polyL}], 
            PlotRange -> {{xmin, xmax}, {ymin, ymax}},  
            Frame -> True, AspectRatio -> Automatic, 
            opt]
];

drawAntiPolygonalRegion[polyL_,  pr_,  color_:  LightPink, 
          opt___] := 
Module[{xmin, xmax, ymin, 
            ymax}, 
{xmin, xmax, ymin, ymax} = 
            Flatten[pr] + {- .01,  .01, - .01,  .01}; 
Show[
            Graphics[{{color, 
                  Rectangle[{pr[[1, 1]], 
                      pr[[2, 1]]}, {pr[[1, 2]], 
                      pr[[2, 2]]}]}, {RGBColor[1, 1, 1], 
                  Polygon  /@  polyL}}], 
            PlotRange -> {{xmin, xmax}, {ymin, ymax}},  
            Frame -> True, AspectRatio -> Automatic, 
            opt]
];

Plot multiple results for polygonal region

PolygonalAllPlot[
          polyL_,  {μ_, Σ_, 
            priors_}] := 
Module[{pr, g1, g2, g3, g4, g5, 
            ga, gb}, 
pr  =  {{0, 1}, {0, 
                1}}; 
g1  =   drawPolygonalRegion[
              polyL,  pr, LightPink,  nF, nS]; 
g2  =  
            PlotMeans[{μ, Σ, priors},  pr,  Blue,  nF, 
              nS]; 
g3  =  
            PlotContourPDFOverlay[{μ, Σ, priors},  pr, 50, 
              nF, nS]; 
g4  =  
            PlotContourPDF[{μ, Σ, priors},  pr, 50, nF, 
              nS]; 
g5  =  
            PlotPDF[{μ, Σ, priors},  pr, 50, 
              Axes -> False, nS]; 
ga  =  
            Show[g1, g2]; 
gb  =  
            Show[g1, g3]; 
Show[
            GraphicsArray[{{ga, gb}, {g4, g5}}], yS,  
            imSize]
];

Plot multiple results for polygonal region

AntiPolygonalAllPlot[
          polyL_,  {μ_, Σ_, 
            priors_}] := 
Module[{pr, g1, g2, g3, g4, g5, 
            ga, gb}, 
pr  =  {{0, 1}, {0, 
                1}}; 
g1  =   drawAntiPolygonalRegion[
              polyL,  pr, LightPink,  nF, nS]; 
g2  =  
            PlotMeans[{μ, Σ, priors},  pr,  Blue,  
              nS]; 
g3  =  
            PlotContourPDFOverlay[{μ, Σ, priors},  pr, 50, 
              nF, nS]; 
g4  =  
            PlotContourPDF[{μ, Σ, priors},  pr, 50, nF, 
              nS]; 
g5  =  
            PlotPDF[{μ, Σ, priors},  pr, 50, 
              Axes -> False, nS]; 
ga  =  
            Show[g1, g2]; 
gb  =  
            Show[g1, g3]; 
Show[
            GraphicsArray[{{ga, gb}, {g4, g5}}], yS, 
            imSize]
];

Data generation and visualization

Idealized data classes

idealPlot  =  specialClassPlot;

Generate data from above distributions

Number of points in each class

n1  =  1000;
n2  =  1000;

Generate data

data1  =  region1Data[n1];
data2  =  region2Data[n2];

Visualize data

gd1  =  
        PlotData2D[data1, pr, style1, nS];
gd2  =  
        PlotData2D[data2, pr, style2, nS];
gd  =  Show[gd1, gd2,  FS];

EM experiment #1: red class, scalar covariances

Number of Gaussians, assumes scalar covariances

Choose number of Gaussians in mixture model

k = 9;

Assumes scalar covariances

initScalar  =  EMInitializationScalar[data1, k];

Compute EM algorithm

Timing[em  =  EM[data1, initScalar, 0.001];]

Visualization

Initial model

PolygonalAllPlot[poly, em // First];

Final model

PolygonalAllPlot[poly, em // Last];

Plot log-likelihood of data given the model as a function of EM iteration

PlotLogLikelihood[em, data1];

Save solution

scalar1  =  em  //  Last;

EM experiment #2: red class, full covariances

Number of Gaussians, assumes full covariance matrices

Choose number of Gaussians in mixture model

k = 4;

Assumes full covariance matrices

initFull  =  EMInitializationFull[data1, k];

Compute EM algorithm

Timing[em  =  EM[data1, initFull, 0.001];]

Visualization

Initial model

PolygonalAllPlot[poly, em // First];

Final model

PolygonalAllPlot[poly, em  //  Last];

Plot log-likelihood of data given the model as a function of EM iteration

PlotLogLikelihood[em, data1];

Save solution

full1  =  em  //  Last;

EM experiment #3: blue class, scalar covariances

Number of Gaussians, assumes scalar covariances

Choose number of Gaussians in mixture model

k = 9;

Assumes scalar covariances

initScalar  =  EMInitializationScalar[data2, k];

Compute EM algorithm

Timing[em  =  EM[data2, initScalar, 0.001];]

Visualization

Initial model

AntiPolygonalAllPlot[poly, em // First];

Final model

AntiPolygonalAllPlot[poly, em // Last];

Plot log-likelihood of data given the model as a function of EM iteration

PlotLogLikelihood[em, data2];

Save solution

scalar2  =  em  //  Last;

EM experiment #4: blue class, full covariances

Number of Gaussians, assumes full covariance matrices

Choose number of Gaussians in mixture model

k = 3;

Assumes full covariance matrices

initFull  =  EMInitializationFull[data2, k];

Compute EM algorithm

Timing[em  =  EM[data2, initFull, 0.001];]

Visualization

Initial model

AntiPolygonalAllPlot[poly, em // First];

Final model

AntiPolygonalAllPlot[poly, em  //  Last];

Plot log-likelihood of data given the model as a function of EM iteration

PlotLogLikelihood[em, data2];

Save solution

full2  =  em  //  Last;

Classification

Scalar covariance models

Which models

model1  =  scalar1;
model2  =  scalar2;

Generate classification data

pts  =  150;
reg1  = 
        MixtureLikelihoodMatrix[{0. , 1. , 0. , 1. }, pts, 
          model1];
reg2  = 
        MixtureLikelihoodMatrix[{0. , 1. , 0. , 1. }, pts, 
          model2];

class  = 
        Partition[If[#[[1]]  <  #[[2]],  1, 
                0] &  /@  ({Flatten[reg1],  Flatten[reg2]}  //  
                Transpose), pts + 1];
classPlot  = 
        ListDensityPlot[class,  
          ColorFunction -> (If[# == 1, RGBColor[0, 0, #],  
                  RGBColor[1, 0, #]] &),  Mesh -> False,  
          MeshRange -> {{0, 1}, {0, 1}}, nF, nS];

g1  =  
        drawPolygonalRegion[poly, {{0, 1}, {0, 1}},  Red,  nF, 
          nS];
g2  =  
        drawAntiPolygonalRegion[poly, {{0, 1}, {0, 1}},  Blue,  nF, 
          nS];
idealPlot  =  Show[g2, g1];

Visualize result

Show[GraphicsArray[{idealPlot,  classPlot}], yS, 
        imSize];

Full covariance models

Which models

model1  =  full1;
model2  =  full2;

Generate classification data

pts  =  150;
reg1  = 
        MixtureLikelihoodMatrix[{0. , 1. , 0. , 1. }, pts, 
          model1];
reg2  = 
        MixtureLikelihoodMatrix[{0. , 1. , 0. , 1. }, pts, 
          model2];

class  = 
        Partition[If[#[[1]]  <  #[[2]],  1, 
                0] &  /@  ({Flatten[reg1],  Flatten[reg2]}  //  
                Transpose), pts + 1];
classPlot  = 
        ListDensityPlot[class,  
          ColorFunction -> (If[# == 1, RGBColor[0, 0, #],  
                  RGBColor[1, 0, #]] &),  Mesh -> False,  
          MeshRange -> {{0, 1}, {0, 1}}, nF, nS];

g1  =  
        drawPolygonalRegion[poly, {{0, 1}, {0, 1}},  Red,  nF, 
          nS];
g2  =  
        drawAntiPolygonalRegion[poly, {{0, 1}, {0, 1}},  Blue, nF, 
          nS];
idealPlot  =  Show[g2, g1];

Visualize result

Show[GraphicsArray[{idealPlot,  classPlot}], yS, 
        imSize];

Animation (slow)

Red class

PolygonalAllPlot[poly, #] &  /@  em;

Blue class

AntiPolygonalAllPlot[poly, #] &  /@  em;