Homework 1 solutions

This notebook  assumes that the definitions in the notebook "statistical_modeling.nb" have already been evaluated.

General

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

Problem 1

The pdf of a d length vector is a joint pdf of the corresponding d random variables. If the variables are independent, then their joint distribution is the product of the marginals.

pX[x1_,  x2_] := KX;
pY[y1_, y2_]  :=  
        KY (E^(((y1  +  2 y2)/4)));
px1[x1_]  :=  Kx1;
px2[x2_]  :=  Kx2;
py1[y1_]  :=  
        Ky1 (E^((y1/4)));
py2[y2_]  := Ky2 (E^(((2 y2)/4)));

Part a) Find (Subscript[K, X]) and (Subscript[K, Y]).

For(Subscript[K, X]),the volume of the function must be 1,so we can say:

KX =. ;
S1  =  Solve[(3 - 2) (3 - 0) KX  ==  1,  KX]  //  
            Flatten  // N;
KX  =  KX  /.  S1

Because the variables are independent, the integration can be separated as follows:
\(TraditionalForm\`∫\_0\%\(Log[2]\)\(∫\_0\%\(Log[2]\)p\_Y[y\_1, \ y\_2] dY\)\) = \(TraditionalForm\`\(K\_Y\) \(∫\_0\%\(Log[2]\)\(e\^\(y\_1\/4\)\) dy\_1 \(∫\_0\%\(Log[2]\)\(e\^\(\(2\ y\_2\)\/4\)\) dy\_2\)\)\)= 1

KY =. ;
S2 = Solve[KY ((∫_0^(Log[2]))(E^((y1/4))) \[DifferentialD]y1) ((∫_0^(Log[2]))(E^(((2 y2)/4))) \[DifferentialD]y2) == 1, KY] // Flatten// N 
KY  =  KY  /.  S2

Part b) Find the mean vectors (no they aren't bullies) (Subscript[μ, X])and (Subscript[μ, Y]).

Again, because the variables are independent, we can solve for each of their means individually. First, we need to find (Subscript[K, (Subscript[x, 1])]), (Subscript[K, (Subscript[x, 2])]), (Subscript[K, (Subscript[y, 1])]), and (Subscript[K, (Subscript[y, 2])]), where (Subscript[K, X]) =  (Subscript[K, (Subscript[x, 1])]) (Subscript[K, (Subscript[x, 2])])and (Subscript[K, Y]) =  (Subscript[K, (Subscript[y, 1])]) (Subscript[K, (Subscript[y, 2])]). Then we solve for the means.

Independent variable method of finding μ.

Kx1 = (1/1)
Kx2 =  (1/3)  // N
Ky1 =  
      (1/((∫_0^(Log[2]))(E^((y1/4))) \[DifferentialD]y1))  //  N
Ky2 =  
      (1/((∫_0^(Log[2]))(E^(((2  y2)/4))) \[DifferentialD]y2 ))  //  N

Just to prove that this worked notice that (Subscript[K, X]) =  (Subscript[K, (Subscript[x, 1])]) (Subscript[K, (Subscript[x, 2])])and (Subscript[K, Y]) =  (Subscript[K, (Subscript[y, 1])]) (Subscript[K, (Subscript[y, 2])]) holds true.

Kx1 Kx2

Ky1 Ky2

μx1  =  
      Kx1 (∫_2^3) x1 \[DifferentialD]x1  // 
        N
μx2  =  
      Kx2 (∫_0^3) x2 \[DifferentialD]x2  // 
        N
μy1  =  
      Ky1 (∫_0^(Log[2]))y1 (E^((y1/4))) \[DifferentialD]y1  // N
μy2  =  
      Ky2 (∫_0^(Log[2]))y2 (E^(((2 y2)/4))) \[DifferentialD]y2  // N

μx  =  {μx1,  μx2}
μy  =  {μy1,  μy2}

Finding μ from the definition

μx1  =  (∫_2^3)(∫_0^3) 
            x1 pX[x1, x2] \[DifferentialD]x2 \[DifferentialD]x1

μx2  =  (∫_2^3)(∫_0^3) 
            x2 pX[x1, x2] \[DifferentialD]x2 \[DifferentialD]x1

μy1  =  (∫_0^(Log[2]))(∫_0^(Log[2]))y1 pY[y1, y2] \[DifferentialD]y2 \[DifferentialD]y1

μy2  =  (∫_0^(Log[2]))(∫_0^(Log[2]))y2 pY[y1, y2] \[DifferentialD]y2 \[DifferentialD]y1

Part c) Find the covariance matrices (Subscript[Σ, X]) and (Subscript[Σ, Y]).

Independent variable method

Because the variables are independent, all the off-diagonal elements are zero. The diagonal elements are as follows:

Σx1x1  =  (∫_2^3)(((x1 - μx1))^2) px1[
              x1] \[DifferentialD]x1  // N
Σx2x2  =  (∫_0^3)(((x2  -  μx2))^2) px2[x2] \[DifferentialD]x2  // N
Σy1y1  =  (∫_0^(Log[2]))(((y1  -  μy1))^2) py1[y1] \[DifferentialD]y1  // N
Σy2y2  =  (∫_0^(Log[2]))(((y2  -  μy2))^2) py2[y2] \[DifferentialD]y2  // N

Σx  =  {{Σx1x1,  
            0},  {0,  Σx2x2}};
Σx  //  MatrixForm
Σy  =  {{Σy1y1,  
            0},  {0,  Σy2y2}};
Σy  // MatrixForm

Definition method

Sx1x1  =  (∫_2^3)(∫_0^3)(((x1  -  μx1))^2) pX[x1, x2] \[DifferentialD]x2 \[DifferentialD]x1

Sx2x2  =  (∫_2^3)(∫_0^3)(((x2  -  μx2))^2) pX[x1, x2] \[DifferentialD]x2 \[DifferentialD]x1

Sx1x2  =  (∫_2^3)(∫_0^3)(x1  -  μx1) (x2  -  μx2) 
              pX[x1, x2] \[DifferentialD]x2 \[DifferentialD]x1  //  
        Chop

Sx2x1  =  Sx1x2;

Sy1y1  =  (∫_0^(Log[2]))(∫_0^(Log[2]))(((y1  -  μy1))^2) pY[y1, 
              y2] \[DifferentialD]y2 \[DifferentialD]y1

Sy2y2  =  (∫_0^(Log[2]))(∫_0^(Log[2]))(((y2  -  μy2))^2) pY[y1, 
              y2] \[DifferentialD]y2 \[DifferentialD]y1

Sy1y2  =  (∫_0^(Log[2]))(∫_0^(Log[2]))(y1  -  μy1) (y2  -  μy2) pY[y1, 
                y2] \[DifferentialD]y2 \[DifferentialD]y1  //  
        Chop

Sy2y1  =  Sy1y2;

ΣX = (({{Sx1x1, Sx1x2}, {Sx2x1, Sx2x2}}));
ΣY = (({{Sy1y1, Sy1y2}, {Sy2y1, Sy2y2}}));
ΣX  // 
      MatrixForm
ΣY  // 
      MatrixForm

Problem 2

Part a) Which is the most likely covariance matrix?

A covariance matrix is always symmetric and positive semidefinite. Symmetric means that element (Subscript[m, ij])  =  (Subscript[m, ji]). Positive semidefinite means that the determinant of a covariance matrix can never be negative.  Note that the off diagonal elements can be positive or negative. To analyze a covariance matrix, beside these two general properties, one must also consider its problem specific properties.

The symmetric requirement eliminates covariance 2 (general). The determinant of covariance 6 is negative, as shown below (general). Altitude and pressure are inversely related, so the off diagonal elements must be negative (problem specific). This eliminates covariances 1, 5, and 8. Covariance 4 can also be eliminated because its off diagonals are small, indicating little relationship between the two variables (problem specific). Covariance 9 can be eliminated because the huge variance in altitude is not followed by the huge variance in pressure (problem specific).  Also, covariance 9's correlation coefficient, (Subscript[ρ, 9]), is low, indicating that the two variables are not very related (problem specific). There are not specific enough details in the problem statement to determine whether covariance 3 or 7 is more correct. Covariances 3 and 7 are equally likely from real data. So the answer is covariances 3 and 7.

Σ6 = (({{1.232, -3.307}, {-3.307, 2.791}}));
Det[Σ6]

ρ9  =  -((Sqrt[0.867])/((Sqrt[1232.5]) (Sqrt[0.91])))

Part b) Find the eigenvalues and eigenvectors for the solution to Part a).

Σ3 = (({{1.232, -0.867}, {-0.876, 3.307}}));
Σ7 = (({{1.232, -0.812}, {-0.812, 2.791}}));
Eigenvalues[Σ3]  // 
      TraditionalForm
Eigenvectors[Σ3]  // 
      TraditionalForm
Eigenvalues[Σ7]  // 
      TraditionalForm
Eigenvectors[Σ7]  // 
      TraditionalForm

Problem 3

Note: You must run statistical_modeling.nb prior to evaluating this section of code.

Plot the contour plot and Bayes decision boundary 2-D normal distributions with the following parameters:

μ1p3 = {0, 1};
μ2p3 = {0, -1};

Default plot range

pr = {{-1, 1}, {-1, 1}};

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};

Part a)

Σ1ap3 = (({{2, 0}, {0, 1}}));
Σ2ap3 = (({{1, 0}, {0, 1}}));

prLocal  =  {{-5, 5},  {-5,  5}};

Log likelihood of first Gaussian

lp1  =  
        logpdf[μ1p3,  Σ1ap3]; 
lp1  //  TraditionalForm

Log likelihood of second Gaussian

lp2  =  
        logpdf[μ2p3,  Σ2ap3];
lp2  //  TraditionalForm

Decision boundary

Simplify[lp1  -  lp2  ==  0]   //  
      TraditionalForm

View Gaussians

ShowTwoGaussianPDFs[Exp[lp1  //  N],  Exp[lp2  //  N],  
        prLocal,  x, y, yS,  imSize];

View decision boundary

ShowDecisionBoundary[μ1p3, Σ1ap3, μ2p3, Σ2ap3, prLocal, classColors, {data1Color, data2Color}, 
        yS];

Part b)

Σ1bp3 = (({{1, 1/2}, {1/2, 1}}));
Σ2bp3 = (({{1, -1/2}, {-1/2, 1}}));

prLocal  =  {{-5, 5},  {-5,  5}};

Log likelihood of first Gaussian

lp1 = 
        logpdf[μ1p3, Σ1bp3] // Rationalize;

Log likelihood of second Gaussian

lp2 =  
        logpdf[μ2p3, Σ2bp3] // Rationalize;

Decision boundary

Simplify[lp1  -  lp2 ==  0]   //  
      TraditionalForm

View Gaussians

ShowTwoGaussianPDFs[Exp[lp1  //  N],  Exp[lp2  //  N],  
        prLocal,  x, y, yS,  imSize];

View decision boundary

Mathematica is a computational tool and cannot deal with singularities. In solving our decision boundary equation, we have a product of two terms which must equal zero. This yields an infinite set of solutions. As a prudent user, you must pay attention for these cases. If you used rational covariance matrices, the ShowDecisionBoundary function works. If you used decimal representation then you needed to write your own function as I did below (largely copied from Michael's ShowDecisionBoundary). Remember that with computers: garbage in -> garbage out. From the boundary equation, the boundary consists of two lines: x = -2, and y = 0.

ShowSpecialCaseBoundary[μ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[(x > -2 &&  
                      y > 0)  ||  (x  < -2 && y < 0),  {x,  
                  x1, x2}, {y, y1, y2},  nS]  /.  
              Hue[1/2] -> cRegion[[1]]; 
          g2 =  InequalityPlot[(x < -2 && 
                      y > 0)  ||  (x > -2  &&  y < 0),  {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, -5, 
                5}, {y, -5, 5},  nS, 
              PlotStyle -> {Thickness[ .015], 
                  cPDF[[2]]}]; 
gf  =  
            Show[g1, g2,  g3, g4, g5, yS,  PlotRange -> plotRange,  
              Frame -> True,  Axes -> False, 
              opts];
];

ShowDecisionBoundary[μ1p3, Σ1bp3, μ2p3, Σ2bp3, prLocal, classColors, {data1Color, data2Color}, 
        yS];

ShowSpecialCaseBoundary[μ1p3, Σ1bp3, μ2p3, Σ2bp3, prLocal, classColors, {data1Color, data2Color}, 
        yS];

Problem 4

Compute the 3-D ML estimates of the mean and covariance for the gator and smilie pictures. Note that there is a difference between the unbiased covariance (which is typically used in statistics) and the ML estimate of the covariance. As the number of samples increases, the ML estimate of the covariance matrix does not converge to the true parameter value (See Duda and Hart textbook pp. 89). Therefore, if you are using a function in a math program, be sure you know what formula it is using. For instance, cov function in Matlab by default computes the unbiased estimate as shown in this excerpt from its help file:
COV Covariance matrix. 
COV(X), if X is a vector, returns the variance. For matrices, where each row is an observation, and each column a variable, COV(X) is the covariance matrix. DIAG(COV(X)) is a vector of variances for each column, and SQRT(DIAG(COV(X))) is a vector of standard deviations. COV(X,Y), where X and Y are vectors of equal length, is equivalent to COV([X(:) Y(:)]). 
COV(X) or COV(X,Y) normalizes by (N-1) where N is the number of observations. This makes COV(X) the best unbiased estimate of the covariance matrix if the observations are from a normal distribution. 
COV(X,1) or COV(X,Y,1) normalizes by N and produces the second moment matrix of the observations about their mean. COV(X,Y,0) is the same as COV(X,Y) and COV(X,0) is the same as COV(X). The mean is removed from each column before calculating the result. 
See also CORRCOEF, STD, MEAN.

Gator Image

\(\(path\  = \ "/Users/nechyba/www/eel6825/assignments/hw1/homework1/";\)\)

gator  =  
        Import[path <> "gator.ppm",  "PPM"];
gatorjpg  =  
        Import[path <> "gator.jpg", "JPG"];
Show[gator];

gatorRGB  =  
        Partition[Flatten[gator[[1, 1]]], 
          3];
gatorRGBjpg  =  
        Partition[Flatten[gatorjpg[[1, 1]]], 3];

Dimensions[gator[[1, 1]]]
Dimensions[gator[[1, 1]]]
n  =  Part[Dimensions[gatorRGB],  1]
μgator  =  (1/n) (∑_(i = 1)^n) 
              gatorRGB[[i]]  //  N
μgator = Mean[gatorRGB] // N
μgatorjpg  =  Mean[gatorRGBjpg]  // N

ΣgatorML  =  (1/n) (∑_(i = 1)^n) 
                Transpose[{(gatorRGB[[i]]  -  μgator)}] . {(gatorRGB[[i]]  -  μgator)} // 
        MatrixForm
ΣgatorML  =  CovarianceMatrixMLE[gatorRGB] // N // 
        MatrixForm
ΣgatorMLjpg  =  CovarianceMatrixMLE[gatorRGBjpg]  // 
          N  // MatrixForm

ΣgatorUNBIASED  =  (1/(n - 
                  1)) (∑_(i = 1)^n) 
                  Transpose[{(gatorRGB[[i]]  -  μgator)}] . {(gatorRGB[[i]]  -  μgator)} // N // 
        MatrixForm
ΣgatorUNBIASED = CovarianceMatrix[gatorRGB]  // N // 
        MatrixForm

Smilie Image

smilie  =  
        Import[path <> "smilie.ppm",  "PPM"];
smiliejpg  =  
        Import[path <> "smilie.jpg", "JPG"];
Show[smilie];

smilieRGB  =  
        Partition[Flatten[smilie[[1, 1]]], 
          3];
smilieRGBjpg  =  
        Partition[Flatten[smiliejpg[[1, 1]]],  3];

Dimensions[smilie[[1, 1]]]
Dimensions[smiliejpg[[1, 1]]]
n  =  Part[Dimensions[smilieRGB],  1]
μsmilie  =  (1/n) (∑_(i = 1)^n) 
              smilieRGB[[i]]  //  N
μsmilie  =  Mean[smilieRGB]  // N
μsmiliejpg  =  Mean[smilieRGBjpg]  // N

ΣsmilieML  =  (1/n) (∑_(i = 1)^n) 
                  Transpose[{(smilieRGB[[i]]  -  μsmilie)}] . {(smilieRGB[[i]]  -  μsmilie)}  // N  // 
        MatrixForm
ΣsmilieML  =  CovarianceMatrixMLE[smilieRGB]  // 
          N  // MatrixForm
ΣsmilieMLjpg  =  CovarianceMatrixMLE[smilieRGBjpg]  // 
          N  // MatrixForm

ΣsmilieUNBIASED  =  (1/(n - 
                  1)) (∑_(i = 1)^n) 
                  Transpose[{(smilieRGB[[i]]  -  μsmilie)}] . {(smilieRGB[[i]]  -  μsmilie)}  // N  // 
        MatrixForm
ΣsmilieUNBIASED  =  CovarianceMatrix[smilieRGB]  // 
          N  // MatrixForm