\documentclass{article}
\usepackage{graphicx}
\usepackage{listings}
\begin{document}

\renewcommand{\thesection}{Problem \arabic{section}:}
\renewcommand{\theenumi}{\alph{enumi}}
\addtolength{\belowcaptionskip}{-5mm}
\def\bSigma{\mbox{\boldmath{$\Sigma$}}}


\title{\underline{Problem Set 1}}
\author{MAS 622J/1.126J: Pattern Recognition and Analysis}
\date{Due: 5:00 p.m. on September 20}
\maketitle

\noindent [Note: All instructions to plot data or write a program
should be carried out using Matlab.  In order to maintain a reasonable level of consistency and simplicity
we ask that you do not use other software tools.]
\newline \newline
%If you discussed the problem set with other members of the class please write their names at the top of your solution.
If you collaborated with other members of the class, please write their names at the end of the assignment. Moreover, you will need to write and sign the following statement: ``In preparing my solutions, I did not look at any old homeworks, copy anybody's answers or let them copy mine."


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Why? [5 points]} 
Limit your answer to this problem to a page. 
\begin{enumerate}
\item Describe an application of pattern recognition related to your
research.  What are the features?  What is the decision to be made?
Speculate on how one might solve the problem.  
\item In the same way, describe an application of pattern recognition
you would be interested in pursuing for fun in your life outside of
work.
\end{enumerate}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Solution%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\noindent {\bf Solution}: Refer to examples discussed in lecture.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Probability Warm-Up [20 points]} 
Let $x$ and $y$ be discrete random variables, and $a$ and $b$ are constant values. Let $\mu_x$ denote the expected value of $x$ and $\sigma_x^2$ denote the variance of $x$. 
Use excruciating detail to answer the following:
\begin{enumerate}
\item Show $\mbox{E}[a x + b y] = a \mbox{E}[x] + b \mbox{E}[y]$.
\item Show $\sigma_x^2 = \mbox{E}[x^2] - \mu_x^2$.
\item Show that independent implies uncorrelated.
\item Show that uncorrelated does not imply independent.
\item Let $z = ax + by$.  Show that if $x$ and $y$ are uncorrelated,
then $\sigma_z^2 =  a^2 \sigma_x^2 + b^2 \sigma_y^2$.
\item Let $x_i ~ (i=1,...,n) $ be random variables independently drawn from the same probability distribution with mean $\mu_x$ and variance $\sigma ^2_x$. 
For the sample mean $\overline x  = \frac{1}{n}\sum\limits_{i = 1}^n {x_i } $, show the following:
(i)  $\mbox{E} [\overline x]$ = $\mu_x$.
(ii) $\mbox{Var} [\overline x] $ (variance of the sample mean) = $\sigma_x^2 /  n$.  Note that this is different from the sample variance $s_n^2  = \frac{1}{n}\sum\limits_{i = 1}^n {(x_i  - \overline x )^2 } $.
\item Let $x_1$ and $x_2$ be independent and identically distributed (i.i.d) continuous random variables.  Can $\mbox{Pr}[x_1 \le x_2]$ be calculated?  If so, find its value.  If not, explain. \small{Hint 1: Remember that for a continuous variable $\mbox{Pr}[x_1 = k]$ = 0, for any value of k. Hint 2: Remember the definition of i.i.d. variables.}

\item Let $x_1$ and $x_2$ be independent and identically distributed discrete random variables.  Can $\mbox{Pr}[x_1 \le x_2]$ be calculated?  If so, find its value.  If not, explain.
\end{enumerate}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Solution%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\noindent {\bf Solution}:
\begin{enumerate}

\item The following is for continuous random variables.  A similar
argument holds for continuous random variables.
\begin{eqnarray*}
\mbox{E}[ax + by] &=& \sum_{x\in X, y\in Y} (ax+by)\,p(x,y) \\
                &=& a\sum_{x\in X, y\in Y} x\,p(x,y) + b \sum_{x\in X, y\in Y} y\,p(x,y) \\
                &=& a\sum_{x\in X} x\,p(x) + b \sum_{y \in Y} y\,p(y) \\
                &=& a\mbox{E}[x] + b\mbox{E}[y]
\end{eqnarray*}

\item Making use of the definition of variance and the previous part, we have:
\begin{eqnarray*}
\sigma_x^2 &=& \mbox{E}[(x-\mu_x)^2] \\
           &=& \mbox{E}[x^2-2\mu_xx + \mu_x^2] \\
           &=& \mbox{E}[x^2] - \mbox{E}[2\mu_xx] + \mbox{E}[\mu_x^2] \\
           &=& \mbox{E}[x^2] - 2\mu_x\mbox{E}[x] + \mu_x^2 \\
           &=& \mbox{E}[x^2] - 2\mu_x\mu_x + \mu_x^2 \\
           &=& \mbox{E}[x^2] - 2\mu_x^2 + \mu_x^2 \\
           &=& \mbox{E}[x^2] - \mu_x^2
\end{eqnarray*}

\item In order to check if two discrete random variables $x$ and $y$ are uncorrelated, we have to prove $\sigma_{xy}$ = 0 (the same holds for continuous random variables.)  

From the previous question:
\begin{eqnarray*}
\sigma_{x,y}^2 &=& \mbox{E}[x y] - \mu_x \mu_y 
\end{eqnarray*}
If two variables are independent:
\begin{eqnarray*}
\mbox{E}[xy] &=& \sum_{x \in X, y\in Y} xy\,p(x,y) \\
             &=& \sum_{x \in X, y\in Y} xy\,p(x)\,p(y) \\
             &=& \sum_{x \in X} x\,p(x) \sum_{y \in Y} y\,p(y) \\
             &=& \mbox{E}[x]\,\mbox{E}[y]
\end{eqnarray*}
Finally,
\begin{eqnarray*}
\sigma_{x,y}^2 &=&\mbox{E}[x]\,\mbox{E}[y] - \mu_x \mu_y = 0 
\end{eqnarray*}


\item To prove this, we need to find one case where 1) $p(x,y) \neq p(x) p(y)$ and 2)  $\sigma_{xy}$ = 0 are satisfied. One possible solution is as follows:

Suppose we have the discrete random variables $x$ and $y$, and we observed all possibilities:

\begin{center}
\begin{tabular}{c | c}
$x$ & $y$\\
\hline
1 & 1\\
1 & 1\\
-1 & 1\\
-1 & 1\\
0 & 0\\
0 & 0\\
\end{tabular}
\end{center}

If we look at the case where $x = 1$ and $y = 1$, 1) is satisfied:
\begin{eqnarray*}
p(x=1 , y =1) &=& \frac{2}{6} = \frac{1}{3}\\
p(x=1) p(y =1) &=& \frac{2}{6} \frac{4}{6} = \frac{2}{9}
\end{eqnarray*}

Now it is easy to verify that 2) is also satisfied:
\begin{eqnarray*}
\sigma_{x,y}^2 &=& \mbox{E}[x y] - \mu_x \mu_y  = 0 - 0 \frac{4}{6} = 0
\end{eqnarray*}


\item Given that $z = ax + by$ and that $x$ and $y$ are uncorrelated, we have
\begin{eqnarray*}
\sigma_z^2 &=& \mbox{E}[(z-\mu_z)^2] \\
           &=& \mbox{E}[z^2] - \mu_z^2 \\
           &=& \mbox{E}[(ax+by)^2] - (a\mu_x + b\mu_y)^2 \\
           &=& \mbox{E}[a^2x^2 + 2abxy + b^2y^2] - (a^2\mu_x^2 + 2ab\mu_x\mu_y + b^2\mu_y^2) \\
           &=& a^2\mbox{E}[x^2] + 2ab\mbox{E}[xy] + b^2\mbox{E}[y^2] - a^2\mu_x^2 - 2ab\mu_x\mu_y - b^2\mu_y^2 \\
           &=& a^2(\mbox{E}[x^2] - \mu_x^2) + 2ab(\mbox{E}[xy] - \mu_x\mu_y) + b^2(\mbox{E}[y^2] - \mu_y^2) \\
           &=& a^2\sigma_x^2 + 2ab\sigma_{xy}^2 + b^2\sigma_y^2 \\
           &=& a^2\sigma_x^2 + b^2\sigma_y^2,
\end{eqnarray*}
where only the last equality depends on $x$ and $y$ being uncorrelated.


\item Using the result of (a) and the fact that ${\mathop{\rm E}\nolimits} [x_i ] = \mu _x $, 

\[
{\mathop{\rm E}\nolimits} [\bar x] = {\mathop{\rm E}\nolimits} [\frac{1}{n}\sum\limits_{i = 1}^n {x_i } ] = \frac{1}{n}\sum\limits_{i = 1}^n {{\mathop{\rm E}\nolimits} [x_i } ] = \frac{1}{n}\;n\;\mu _x  = \mu _x 
\]

Also, using the result of (d) and the fact $ {\rm{Var}}[x_i ] = \sigma _x^2 $
\[
{\rm{Var}}[\bar x] = {\rm{Var}}[\frac{1}{n}\sum\limits_{i = 1}^n {x_i } ] = \frac{1}{{n^2 }}\sum\limits_{i = 1}^n {{\rm{Var}}[x_i } ] = \frac{1}{{n^2 }}\;n\;\sigma _x^2  = \sigma _x^2 /n
\]


\item Given that $x_1$ and $x_2$ are continuous random variables, we
know that $\mbox{Pr}[x_1 = k] = 0$ and $\mbox{Pr}[x_2 = k] = 0$ for
any value of $k$.  Thus, $$\mbox{Pr}[x_1 \le x_2] = \mbox{Pr}[x_1 < x_2].$$
Given that $x_1$ and $x_2$ are i.i.d., we know that
replacing $x_1$ with $x_2$ and $x_2$ with $x_1$ will have no effect on
the world.  In particular, we know that 
$$\mbox{Pr}[x_1 < x_2] = \mbox{Pr}[x_2 < x_1].$$
However, since probabilities or the space of possible values must sum to one, we have 
$$\mbox{Pr}[x_1 < x_2] + \mbox{Pr}[x_2 < x_1] = 1.$$
Thus, $$\mbox{Pr}[x_1 \le x_2] = \frac{1}{2}.$$

\item For discrete random variables, unlike the continuous case above,
we need to know the distributions of $x_1$ and $x_2$ in order to find
$\mbox{Pr}[x_1 = k]$ and $\mbox{Pr}[x_2 = k]$.  Thus, the argument we
used above fails.  In general, it is not possible to find
$\mbox{Pr}[x_1 \le x_2]$ without knowledge of the distributions of both
$x_1$ and $x_2$.

\end{enumerate}







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Teatime with Gauss and Bayes [20 points]} 
Let $p(x,y) = \frac{1}{2\pi\alpha\beta}
e^{-\left(\frac{(y-\mu)^2}{2\alpha^2} + \frac{(x-y)^2}{2\beta^2}\right)}$

\begin{enumerate}
\item Find $p(x)$, $p(y)$, $p(x|y)$, and $p(y|x)$.  In addition, give
a brief description of each of these distributions.
\item Let $\mu = 0$, $\alpha=15$, and $\beta = 3$.  Plot $p(y)$ and
\hbox{$p(y|x=9)$} for a reasonable range of $y$.  What is the
difference between these two distributions?
\end{enumerate}


{\bf Solution}:
\begin{enumerate}
\item To find $p(y)$, simply factor $p(x,y)$ and then integrate over $x$:
\begin{eqnarray*}
p(y) &=& \int_{-\infty}^\infty\!p(x,y)\,dx \\ 
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\left(\frac{(y-\mu)^2}{2\alpha^2} + \frac{(x-y)^2}{2\beta^2}\right)}\,dx \\
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\frac{(y-\mu)^2}{2\alpha^2}}\,e^{-\frac{(x-y)^2}{2\beta^2}}\,dx \\
     &=& \frac{1}{\sqrt{2\pi\alpha^2}}e^{-\frac{(y-\mu)^2}{2\alpha^2}}\,\int_{-\infty}^\infty\!\frac{1}{\sqrt{2\pi\beta^2}}e^{-\frac{(x-y)^2}{2\beta^2}}\,dx \\
     &=& \frac{1}{\sqrt{2\pi\alpha^2}}e^{-\frac{(y-\mu)^2}{2\alpha^2}} \\
     &=& {\cal N}(\mu,\alpha^2)
\end{eqnarray*}
The integral goes to $1$ because it is of the form of a probability
distribution integrated over the entire domain.  To find $p(x|y)$,
divide $p(x,y)$ by $p(y)$:
\begin{eqnarray*}
p(x|y) &=& \frac{p(x,y)}{p(y)} \\
       &=& \frac{1}{\sqrt{2\pi\beta^2}}e^{-\frac{(x-y)^2}{2\beta^2}} \\
       &=& {\cal N}(y,\beta^2)
\end{eqnarray*}
Finding $p(x)$ and $p(y|x)$ follows essentially the same procedure,
but the algebra is more involved and requires completing the square in
the exponent.
\begin{eqnarray*}
p(x) &=& \int_{-\infty}^\infty\!p(x,y)\,dy \\ 
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\left(\frac{(y-\mu)^2}{2\alpha^2} + \frac{(x-y)^2}{2\beta^2}\right)}\,dy \\
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\left(\frac{\beta^2(y-\mu)^2 + \alpha^2(x-y)^2}{2\alpha^2\beta^2}\right)}\,dy \\
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\left(\frac{\beta^2y^2 -2\beta^2\mu y + \beta^2\mu^2 + \alpha^2x^2 -2\alpha^2xy +\alpha^2y^2}{2\alpha^2\beta^2}\right)}\,dy \\
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\left(\frac{(\alpha^2 + \beta^2)y^2 -2(\alpha^2x + \beta^2\mu)y + (\beta^2\mu^2 + \alpha^2x^2)}{2\alpha^2\beta^2}\right)}\,dy \\
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\left(\frac{y^2 -2\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}y + \frac{\beta^2\mu^2 + \alpha^2x^2}{\alpha^2 + \beta^2}}{2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}\right)}\,dy \\
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\left(\frac{y^2 -2\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}y + \left(\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2 - \left(\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2 + \frac{\beta^2\mu^2 + \alpha^2x^2}{\alpha^2 + \beta^2}}{2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}\right)}\,dy \\
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}e^{-\left(\frac{\left(y - \frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2 - \left(\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2 + \frac{\beta^2\mu^2 + \alpha^2x^2}{\alpha^2 + \beta^2}}{2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}\right)}\,dy \\
     &=& \int_{-\infty}^\infty\!\frac{1}{2\pi\alpha\beta}
         e^{-\left(\frac{\left(y - \frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2}{2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}} \right)} \, 
         e^{-\left(\frac{\frac{\beta^2\mu^2 + \alpha^2x^2}{\alpha^2 + \beta^2} - \left(\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2}
         {2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}\right)}\,dy \\
     &=& \frac{1}{2\pi\alpha\beta}\sqrt{2\pi\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}
         e^{-\left(\frac{\frac{\beta^2\mu^2 + \alpha^2x^2}{\alpha^2 + \beta^2} - \left(\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2}
         {2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}\right)}
         \int_{-\infty}^\infty\! \frac{1}{\sqrt{2\pi\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}}
         e^{-\left(\frac{\left(y - \frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2}{2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}} \right)}\,dy  \\
	( Note &:&   \int_{-\infty}^\infty\! \frac{1}{\sqrt{2\pi\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}}
         e^{-\left(\frac{\left(y - \frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2}{2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}} \right)}\,dy = \int_{-\infty}^\infty\!p(y | x) \,dy  ) \\
     &=& \frac{1}{\sqrt{2\pi(\alpha^2 + \beta^2)}} 
         e^{-\left(\frac{\frac{\beta^2\mu^2 + \alpha^2x^2}{\alpha^2 + \beta^2} - \left(\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2}
         {2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}\right)} \\
     &=& \frac{1}{\sqrt{2\pi(\alpha^2 + \beta^2)}}
         e^{-\left(\frac{(\alpha^2 + \beta^2)(\beta^2\mu^2 + \alpha^2x^2) - (\alpha^2x + \beta^2\mu)^2}
         {2\alpha^2\beta^2(\alpha^2 + \beta^2)}\right)} \\
     &=& \frac{1}{\sqrt{2\pi(\alpha^2 + \beta^2)}}
         e^{-\left(\frac{\alpha^2\beta^2\mu^2 + \alpha^4x^2 + \beta^4\mu^2 + \alpha^2\beta^2x^2 - \alpha^4x^2 - 2\alpha^2\beta^2\mu x - \beta^4\mu^2}
         {2\alpha^2\beta^2(\alpha^2 + \beta^2)}\right)} \\
     &=& \frac{1}{\sqrt{2\pi(\alpha^2 + \beta^2)}}
         e^{-\left(\frac{\alpha^2\beta^2x^2 - 2\alpha^2\beta^2\mu x + \alpha^2\beta^2\mu^2}
         {2\alpha^2\beta^2(\alpha^2 + \beta^2)}\right)} \\
     &=& \frac{1}{\sqrt{2\pi(\alpha^2 + \beta^2)}}
         e^{-\left(\frac{(x - \mu)^2}
         {2(\alpha^2 + \beta^2)}\right)} \\
     &=& {\cal N}(\mu,\alpha^2 + \beta^2)
\end{eqnarray*}
To find $p(y|x)$ we simply divide $p(x,y)$ by $p(x)$.  In finding
$p(x)$, we already know the form of $p(y|x)$ (see the longest line in
the derivation of $p(x)$ above):

\begin{eqnarray*}
p(y|x) &=& \frac{p(x,y)}{p(x)} \\
       &=& \frac{1}{\sqrt{2\pi\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}}}
         e^{-\left(\frac{\left(y - \frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2}\right)^2}{2\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2}} \right)} \\
       &=& {\cal N}(\frac{\alpha^2x + \beta^2\mu}{\alpha^2 + \beta^2},\frac{\alpha^2\beta^2}{\alpha^2 + \beta^2})
\end{eqnarray*}
Note that all the above distibutions are Gaussian.


\item The following Matlab code produced Figure \ref{fig:GaussBayes}:

\lstinputlisting{GaussBayes.m}
%
%\begin{verbatim}
%m = 0.0
%a = 15.0 
%b = 3 
%x = 9 
%
%y = -100:1:100
%mean = ((a^2)*x + (b^2)*m)/(a^2 + b^2)
%var = ((a*b)^2)/(a^2 + b^2)
%p_y_given_x = (1.0/sqrt(2*pi*var))*exp(-((y-mean).^2)/(2*var))
%var2 = a^2
%p_y = (1.0/sqrt(2*pi*var2))*exp(-((y-m).^2)/(2*var2))
%
%hold off
%plot(y,p_y_given_x,'b')
%hold on
%plot(y,p_y,'r')
%legend('p(y|x)', 'p(y)')
%sy = size(y)
%axis([y(1),y(sy(2)),0,0.2])
%xlabel('y')
%text(-70,0.14, '\mu=0')
%text(-70,0.12, '\alpha=15')
%text(-70,0.1, '\beta=3')
%\end{verbatim}

\begin{figure}[tb]
  \centering 
  \includegraphics[width=\textwidth, bb=0 0 375 375]{gb.pdf}
\caption{The marginal p.d.f. of $y$ and the p.d.f. of $y$ given $x$
for a specific value of $x$.  Notice how knowing $x$ makes your
knowledge of $y$ more certain.}
\label{fig:GaussBayes}
\end{figure}

\end{enumerate}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Covariance Matrix [15 points]} 
Let $\bSigma = \left[\begin{array}{rr} 5 & 4 \\ 4 & 5
\end{array} \right]$

\begin{enumerate}
\item Find the eigenvalues and eigenvectors of $\bSigma$ by hand (include all calculations.) Verify your computations with MATLAB function \textit{eig}.
\item Verify that $\bSigma$ is a valid covariance matrix.
\item We provide 500 data points sampled from the distribution ${\cal
N}([0\, 0],\bSigma)$. Download the dataset from the course website and project the data onto the eigenvectors of the covariance matrix. What is the effect of this projection? Include MATLAB code and plots before and after the projection.
\end{enumerate}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Solution%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Solution}:
\begin{enumerate}


\item We can find the eigenvectors and eigenvalues of $\bSigma$ by
starting with the definition of an eigenvector. Namely, a vector ${\bf e}$
is an eigenvector of $\bSigma$ if it satisfies
$$\bSigma{\bf e} = \lambda{\bf e}$$ for some constant scalar
$\lambda$, which is called the eigenvalue corresponding to ${\bf e}$.
This can be rewritten as
$$(\bSigma - \lambda I){\bf e} = 0$$
This is equivalent to 
$$\mbox{det}(\bSigma - \lambda I) = 0$$
Thus, we require that 
$$(5 - \lambda)^2 - 4^2 = 0$$ By inspection, this is true when
$\lambda=9$ and $\lambda=1$.

To find the eigenvectors, we plug the eigenvalues back into the equation above to get
$$(\bSigma - 9 I){\bf e} = 
\left[\begin{array}{cc}5-9 & 4 \\ 4 & 5-9\end{array} \right]
\left[\begin{array}{r}a \\ b\end{array} \right] = 
\left[\begin{array}{rr}-4 & 4 \\ 4 & -4\end{array} \right]
\left[\begin{array}{r}a \\ b\end{array} \right] = 
\left[\begin{array}{r}0 \\ 0\end{array} \right]$$

which gives $a=b$.  Normalized, this results in the eigenvector
$${\bf e}_1 = \left[\begin{array}{r}\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{array} \right]$$
Similarly, $\lambda = 1$ gives
$$(\bSigma - 1 I){\bf e} = 
\left[\begin{array}{cc}5-1 & 4 \\ 4 & 5-1\end{array} \right]
\left[\begin{array}{r}a \\ b\end{array} \right] = 
\left[\begin{array}{rr}4 & 4 \\ 4 & 4\end{array} \right]
\left[\begin{array}{r}a \\ b\end{array} \right] = 
\left[\begin{array}{r}0 \\ 0\end{array} \right]$$
which gives $a=-b$.  Normalized, this results in the eigenvector
$${\bf e}_1 = \left[\begin{array}{r}\frac{-1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{array} \right]$$

%\item The following Matlab program prints out the eigenvectors and eigenvalues of $\bSigma$:
%%\lstinputlisting{eigen.m}

\item The matrix $\bSigma$ is a valid covariance matrix if it is
symmetric and positive semi-definite.  Clearly, it is symmetric, since
$\bSigma^T = \bSigma$.  One way to prove it is positive
semi-definite is to show that all its eigenvalues are non-negative.
This is indeed the case, as shown in the previous question.


\item Figure \ref{fig:covariance} shows how projecting the data onto the eigenvectors of the data's covariance matrix reduces the correlation between $x$ and $y$. The MATLAB code is as follows:
\lstinputlisting{cov_problem.m}

\begin{figure}[tb]
  \centering 
  \includegraphics[width=\textwidth]{points_projection.pdf}
%  \includegraphics[width=\textwidth]{trans_datadots.png}
\caption{The original data and the data transformed into the
coordinate system defined by the eigenvectors of their covariance
matrix.}
\label{fig:covariance}
\end{figure}

\end{enumerate}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Probabilistic Modeling [20 points]}
Let $x \in \{0,1\}$ denote a person's affective state ($x$ = 0 for ``positive-feeling state", and $x$ = 1 for ``negative-feeling state").  The person feels positive with probability $\theta_1$.  Suppose that an affect-tagging system (or a robot) recognizes her feeling state and reports the observed state, $y$, to you.  But this system is unreliable and obtains the correct result with probability $\theta_2$.

\begin{enumerate}
\item Represent the joint probability distribution $P(x, y | \theta)$ for all $x, y$ (a 2x2 matrix) as a function of the parameters $\theta = (\theta_1, \theta_2)$. 

\item The Maximum Likelihood estimation criterion for the parameter $\theta$ is defined as:
\[
\widehat\theta _{ML}  = \arg \mathop {\max }\limits_\theta  L(t_1 ,...,t_n ;\theta ) = \arg \mathop {\max }\limits_\theta  \prod\limits_{i = 1}^n {p(t_i |\theta )} 
\]
where we have assumed that each data point $t_i$ is drawn independently from the same distribution so that the likelihood of the data is $L(t_1 ,...,t_n ;\theta ) = \prod\limits_{i = 1}^n {p(t_i |\theta )} $.  Likelihood is viewed as a function of the parameters, which depends on the data.  Since the above expression can be technically challenging, we maximize the log-likelihood $ \log L(t_1 ,...,t_n ;\theta )$ instead of likelihood.  Note that any monotonically increasing function (i.e., log function) of the likelihood has the same maxima. Thus,
\[
\widehat\theta _{ML}  = \arg \mathop {\max }\limits_\theta  \log L(t_1 ,...,t_n ;\theta ) = \arg \mathop {\max }\limits_\theta  \sum\limits_{i = 1}^n {\log p(t_i |\theta )} 
\]

Suppose we get the following joint observations $t = (x, y)$. 
\begin{center}
\begin{tabular}{c | c}
$x$ & $y$\\
\hline
1 & 0\\
1 & 1\\
0 & 0\\
1 & 1\\
0 & 0\\
0 & 1\\
0 & 0 \\
1 & 1 \\
\end{tabular}
\end{center}

What are the maximum-likelihood (ML) values of $\theta_1$ and $\theta_2$? ({\it Hint.} Since $P(x,y | \theta) = P(y|x, \theta_2) P(x | \theta_1)$, the estimation of the two parameters can be done separately in the log-likelihood criterion.)

\end{enumerate}



{\bf Solution}: 
\begin{enumerate}
\item The probability mass function (pmf) of $x \in \{0,1\}$ is 
\[
P(x) = \left\{ \begin{array}{l}
 \theta _1 ,\quad \quad \quad x = 0 \\ 
 1-\theta _1 ,\quad \quad x = 1 \\ 
 \end{array} \right\}
\]

The conditional pmf of $y \in \{0,1\}$ given that $x=0$ is 
\[
P(y|x = 0) = \left\{ \begin{array}{l}
 \theta _2 ,\quad \quad y = 0 \\ 
 1 - \theta _2 ,\;\;\;y = 1 \\ 
 \end{array} \right\}
\]

The conditional pmf of $y$ given that $x=1$ is 
\[
P(y|x = 1) = \left\{ \begin{array}{l}
 1 - \theta _2 ,\quad y = 0 \\ 
 \theta _2 ,\;\;\;\quad \;\,y = 1 \\ 
 \end{array} \right\}
\]

Use $P(x,y)=P(y|x)P(x)$ to tabulate the joint pmf of $(x,y)$.
\[
\begin{array}{l}
 P(x,y) = \left( {\begin{array}{*{20}c}
   {P(0,0)} & {P(0,1)}  \\
   {P(1,0)} & {P(1,1)}  \\
\end{array}} \right) = \left( {\begin{array}{*{20}c}
   {\theta _2  \theta _1 } & {(1 - \theta _2 ) \theta _1 }  \\
   {(1 - \theta _2 ) (1 - \theta _1) } & {\theta _2 (1 - \theta _1) }  \\
\end{array}} \right) \\ 
  \\ 
 \end{array}
\]

\item We select $(\theta_1, \theta_2)$ to maximize the log-likelihood of the samples $\{ (x_i, y_i), i = 1, ..., n\}$ which may be expressed as 
\begin{eqnarray*}
 J(\theta _1 ,\theta _2 ) &=& \sum\limits_i {\log P(x_i ,y_i )}  \\ 
  &=& \sum\limits_i {(\log P(y_i |x_i ) + \log P(x_i ))}  \\ 
  &=& \left( {\sum\limits_i {\log P(y_i |x_i )} } \right) + \left( {\sum\limits_i {\log P(x_i )} } \right) \\ 
  &=& J_2 (\theta _2 ) + J_1 (\theta _1 ) \\ 
\end{eqnarray*}
Hence, we choose $\theta_1$ to maximize
\begin{eqnarray*}
 J_1 (\theta _1 ) &=&\sum\limits_i {\log P(x_i )}  \\ 
  &=& N(x = 1)\log (1- \theta _1)  + (n - N(x = 1))\log \theta _1  
 \end{eqnarray*}
where $N(x = 1) = \sum\nolimits_i {x_i }$.  Differentiating w.r.t. $\theta_1$ gives
\[
\frac{{\partial J_1 }}{{\partial \theta _1 }} = \frac{{-N(x = 1)}}{{1-\theta _1 }} + \frac{{n - N(x = 1)}}{{ \theta _1 }}
\]
We set this derivative to zero and solve for $\theta_1$ to obtain 
\[
\widehat\theta _1  =1- \frac{{N(x = 1)}}{n}
\]
Similarly, we choose $\theta_2$ to maximize
\begin{eqnarray*}
 J_2 (\theta _2 ) &=& \sum\limits_i {\log P(y_i |x_i )}  \\ 
  &=& N(x = y)\log \theta _2  + (n - N(x = y))\log (1 - \theta _2 ) 
\end{eqnarray*}
where $N(x = y) = \sum\nolimits_i {(x_i y_i  + (1 - x_i )(1 - y_i ))} $.  Differentiating $J_2$ w.r.t. $\theta_2$, setting to zero and solving for $\theta_2$ gives
 \[
\widehat\theta _2  = \frac{{N(x = y)}}{n}
\]
For the example data, $ \widehat\theta _1  = \frac{4}{8}$,  $\widehat\theta _2  = \frac{6}{8}$.
Thus,
\[
\begin{array}{l}
 {\widehat P(x,y)} = \left( {\begin{array}{*{20}c}
   {\widehat\theta _2  \widehat\theta _1 } & {(1 - \widehat\theta _2 ) \widehat\theta _1 }  \\
   {(1 - \widehat\theta _2 ) (1 - \widehat\theta _1) } & {\widehat\theta _2 (1 - \widehat\theta _1) }  \\
\end{array}} \right) \\ 
  \\ 
 \end{array}
\]
The maximum likelihood of the data under this model is
\[
\prod\limits_i {\widehat P(x_i, y_i )}  = \left( {\frac{6}{8}} \right)^6  \left( {\frac{4}{8}} \right)^8 \left( {\frac{2}{8}} \right)^2 \approx 4.34\times 10^{ - 5} 
\]
\end{enumerate}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Ring Problem [20 points]}
To get credit for this problem, you must not only write your own
correct solution, but also write a computer simulation of the process of playing this game:

\vspace{3mm}

Suppose I hide the ring of power in one of three identical boxes while
you weren't looking.  The other two boxes remain empty.  After hiding
he ring of power, I ask you to guess which box it's in.  I know which
box it's in and, after you've made your guess, I deliberately open the
lid of an empty box, which is one of the two boxes you did not choose.
Thus, the ring of power is either in the box you chose or the
remaining closed box you did not choose.  Once you have made your
initial choice and I've revealed to you an empty box, I then give you
the opportunity to change your mind -- you can either stick with your
original choice, or choose the unopened box.  You get to keep the
contents of whichever box you finally decide upon.

\vspace{3mm}

\begin{itemize}
\item What choice should you make in order to maximize your chances of receiving
the ring of power?  Justify your answer using Bayes' rule.
\item Write a simulation.  There are two choices in this game for the
contestant in this game: (1) choice of box, (2) choice of whether or not to
switch.  In your simulation, first let the host choose a random box to place the ring of power.  Show a trace of your program's output for a single game play,
as well as a cumulative probability of winning for 1000 rounds of the two
policies (1) to choose a random box and then switch and (2) to choose a 
random box and not switch.
\end{itemize}

%\section{Python Prime Practice}
%[{\bf Extra Credit}] \\
%A billboard found on Massachusetts Avenue near MIT states:
%
%$$\left\{ \begin{array}{l}
%          \mbox{first 10-digit prime found} \\
%          \mbox{in consecutive digits of $e$} 
%	  \end{array}
%  \right\}\mbox{\bf.com}$$
%
%Use this problem as an excuse to hone your Python skills.  Yes, the
%problem and solution are both well-documented on the Internet.
%However, the point of this exercise is for you to understand how to
%solve it using Python.  Your work should be in the form of a Python
%program, a copy of that program's output, and whatever else you wish
%to include.  Your work should show that you understand Python well
%enough to solve the problem, not that you are capable of copying and
%pasting the results of an Internet search.  That said, showing that
%you understand what you found through an Internet search is
%acceptable.

{\bf Solution}:

\begin{itemize}
\item Always switch your answer to the box you didn't
choose the first time.  This reason is as follows.  You have a 1/3 chance of
initially picking the correct box.  That is, there is a 2/3 chance the
correct answer is one of the other two boxes.  Learning which of the
two other boxes is empty does not change these probabilities; your
initial choice still has a 1/3 chance of being correct.  That is,
there is a 2/3 chance the remaining box is the correct answer.
Therefore you should change your choice.


%More formally, 
%
%event\_right\_first\_choice = the event that your first choice is right
%
%event\_wrong\_first\_choice = the event that your first choice is wrong
%
%event\_right\_when\_change = the event that you get the ring when changing your initial choice 
%
%event\_an\_empty\_box\_opened = the event that an empty box is opened after your first choice  
%
%First, P(event\_right\_first\_choice) = 1/3
%
%Second,  P(event\_right\_when\_change $|$ event\_an\_empty\_box\_opened) 
%
%%= P(event\_right\_when\_change, (event\_right\_first\_choice or event\_wrong\_first\_choice) $|$ event\_an\_empty\_box\_opened) 
%
%= P(event\_right\_when\_change, event\_right\_first\_choice $|$  event\_an\_empty\_box\_opened) + P(event\_right\_when\_change, event\_wrong\_first\_choice $|$ event\_an\_empty\_box\_opened)
%
%%= P(event\_right\_when\_change $|$  event\_right\_first\_choice, event\_an\_empty\_box\_opened) P(event\_right\_first\_choice $|$ event\_an\_empty\_box\_opened) + P(event\_right\_when%\_change $|$  event\_wrong\_first\_choice, event\_an\_empty\_box\_opened) P(event\_wrong\_first\_choice $|$ event\_an\_empty\_box\_opened)		   
%
%= P(event\_right\_when\_change $|$  event\_right\_first\_choice, event\_an\_empty\_box\_opened) P(event\_right\_first\_choice) + P(event\_right\_when\_change $|$  event\_wrong\_first\_choice, event\_an\_empty\_box\_opened) P(event\_wrong\_first\_choice)		   
%
%= 0*1/3 + 1*2/3 = 2/3

%Note that P(event\_right\_first\_choice $|$  event\_an\_empty\_box\_opened) = P(event\_right\_first\_choice), 
%P(event\_wrong\_first\_choice $|$  event\_an\_empty\_box\_opened) = P(event\_wrong\_first\_choice), because your first choice is independent of the event "event\_an\_empty\_box\_opened".    

%Thus,  P(event\_right\_when\_change $|$ event\_an\_empty\_box\_opened) $>$ P(event\_right\_first\_choice)

Using Bayes:
\newline \newline
R = Box with the Ring. \newline
O = Box Opened. \newline 
S = Box Selected. \newline
\begin{eqnarray}
P(R | O, S) =  \frac{P(O | R, S) P(R | S)}{P(O | S)}
\end{eqnarray}
Where:
\begin{eqnarray}
P(O | S) =  \Sigma^{3}_{R=1} P(O, R | S) = \Sigma^{3}_{R=1} P(O | R, S) P(R | S)
\end{eqnarray}
\begin{eqnarray}
P(R | O, S) =  \frac{1*\frac{1}{3}}{\frac{1}{2}*\frac{1}{3} + 1 * \frac{1}{3} + 0 * \frac{1}{3}} = \frac{2}{3}
\end{eqnarray}

Another way to understand the problem is to extend it to 100 boxes,
only one of which has the ring of power.  After you make your initial
choice, I then open 98 of the 99 remaining boxes and show you that
they are empty.  Clearly, with very high probability the ring of power
resides in the one remaining box you did not initially choose.

\item 

Here is a sample simulation output for the Ring problem:

\lstinputlisting{montysimulation.txt}

Here is a Matlab program that simulates the Ring simulation output above:

\lstinputlisting{montyhall.m}



\end{itemize}



\end{document}
