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



\begin{document}

\renewcommand{\thesection}{Problem \arabic{section}:}
\renewcommand{\theenumi}{\alph{enumi}}
\addtolength{\belowcaptionskip}{-5mm}

\title{\underline{Problem Set 1 Solutions}}
\author{MAS.622J/1.126J: Pattern Recognition and Analysis}
\date{Originally Due Monday, 15 September 2008}
\maketitle

\section{Why?}
\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.  Limit your answer to a
page.
\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}

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

\section{Probability Warm-Up}
Let $X$ and $Y$ be random variables.  Let $\mu_X \equiv
\mbox{E}[X]$ denote the expected value of $X$ and $\sigma_X^2 \equiv
\mbox{E}[(X - \mu_X)^2] = \mbox{E}[X^2] - \mu_X^2$ denote the variance of $X$.  $a$ and $b$ are constant values. 
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  The conditional expected value ${\mbox E} ( X | Y )$ is a random variable in its own right, whose value depends on the value of Y. Notice that the conditional expected value of X given the event Y = y is a function of y. If we write $E( X | Y = y) = g(y)$ then the random variable $ {\mathop{\mbox E}\nolimits} [X|Y = y] = \sum\limits_x {xP(x|Y = y)} $ is just $g(Y)$.
Show ${\mathop{\mbox E}\nolimits} [X] = {\mathop{\mbox E}\nolimits} [{\mathop{\mbox E}\nolimits} [X|Y]]$ and 
${\mathop{\mbox E}\nolimits} [Y] = {\mathop{\mbox E}\nolimits} [{\mathop{\mbox E}\nolimits} [Y|X]]$.

\item For a real value function $f$ and discrete random variables $X$ and $Y$,   ${\mathop{\mbox E}\nolimits} [f(X,Y)] = \sum\limits_x {\sum\limits_y {f(x,y)P(x,y)} } $. Show ${\mathop{\mbox E}\nolimits} [f(X,Y)] = {\mathop{\mbox E}\nolimits} [{\mathop{\mbox E}\nolimits} [f(X,Y)|Y]] $.

 
\item Let $X_1$ and $X_2$ be independent and identically distributed continuous random variables.  Can $\mbox{Pr}[X_1 \le X_2]$ be calculated?  If so, find its value.  If not, explain.
\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}



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

\item The following is for continuous random variables.  A similar
argument holds for discrete random variables.
\begin{eqnarray*}
\mbox{E}[aX + bY] &=& \int\!\!\int\!(ax+by)\,p(x,y)\,dx\,dy \\
                &=& a\int\!\!\int\!x\,p(x,y)\,dx\,dy + b\int\!\!\int\!y\,p(x,y)\,dx\,dy \\
                &=& a\int\!x\,p(x)\,dx + b\int\!y\,p(y)\,dy \\
                &=& 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 For uncorrelated $\sigma_{xy}$ = 0.  Let $X$ and $Y$ be independent continuous random variables (a
similar argument holds for discrete random variables).  Then,
\begin{eqnarray*}
\mbox{E}[XY] &=& \int\!\!\int\!xy\,p(x,y)\,dx\,dy \\
             &=& \int\!\!\int\!xy\,p(x)\,p(y)\,dx\,dy \\
             &=& \int\!x\,p(x)\,dx \int\!y\,p(y)\,dy \\
             &=& \mbox{E}[X]\,\mbox{E}[Y]
\end{eqnarray*}




\item Let $X$ and $Y$ be discrete random variables such that $X$ takes
on values from $\{0,1\}$ and $Y$ takes on values from $\{-1,0,1\}$.
Let the probability mass function of $X$ be

\begin{eqnarray*}
p_{x}[x=0] &=& 0.5 \\ 
p_{x}[x=1] &=& 0.5
\end{eqnarray*}

\noindent and the probability mass function of $Y$ conditioned on $X$ be

\begin{eqnarray*}
  p_{y|x}[y=-1|x=0] &=& 0.5 \\
  p_{y|x}[y=0|x=0] &=& 0 \\
  p_{y|x}[y=1|x=0] &=& 0.5 \\
  p_{y|x}[y=-1|x=1] &=& 0 \\
  p_{y|x}[y=0|x=1] &=& 1 \\
  p_{y|x}[y=1|x=1] &=& 0.
\end{eqnarray*}

\noindent Given the above, and the fact that $p_{x,y}[x,y] =
p_{y|x}[y|x]\,p_x[x]$, we get
\begin{eqnarray*}
p_{x,y}[x=0,y=-1] &=& 0.25 \\
p_{x,y}[x=0,y=0] &=& 0 \\
p_{x,y}[x=0,y=1] &=& 0.25 \\
p_{x,y}[x=1,y=-1] &=& 0 \\
p_{x,y}[x=1,y=0] &=& 0.5 \\
p_{x,y}[x=1,y=1] &=& 0.
\end{eqnarray*}

\noindent However, the product of the marginals is given by
\begin{eqnarray*}
p_x[x=0]\,p_y[y=-1] &=& 0.125 \\
p_x[x=0]\,p_y[y=0] &=& 0.25 \\
p_x[x=0]\,p_y[y=1] &=&  0.125 \\
p_x[x=1]\,p_y[y=-1] &=&  0.125 \\
p_x[x=1]\,p_y[y=0] &=& 0.25 \\
p_x[x=1]\,p_y[y=1] &=& 0.125.
\end{eqnarray*}

\noindent Thus, we see that $p_{x,y}[x,y] \neq p_x[x]\,p_y[y]$ and $X$
and $Y$ are not independent.  independent.  However, since $XY$ is
identically zero, we also get
\begin{eqnarray*}
\mbox{cov}(X,Y) = \sigma_{XY}^2 &=& \mbox{E}[(X-\mu_X)(Y-\mu_Y)] \\
              &=& \mbox{E}[XY] - \mu_X \mu_Y \\
              &=& \mbox{E}[0] - (0.5)(0) \\
              &=& 0 - 0 \\
              &=& 0.
\end{eqnarray*}
Therefore, $X$ and $Y$ are uncorrelated but not independent.

\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 ${\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 The following is for discrete random variables.  A similar
argument holds for continuous random variables.
$\mbox{E}[X|Y] = \mbox{E}[X|Y = y] $ is a function of $y$, i.e., 
$\mbox{E}[X|Y = y] = \sum\limits_x {xP(x|Y = y)} = g(y)$
\begin{eqnarray*}
 \mbox{E}[\mbox{E}[X|Y]] &=& \mbox{E}[g(y)] \\ 
  &=& \sum\limits_y {g(y)P(y)}  \\ 
  &=& \sum\limits_y {\sum\limits_x {xP(x|Y = y)} P(y)}  \\ 
  &=& \sum\limits_y {\sum\limits_x {xP(x,y)} }  \\ 
  &=& \sum\limits_x {x\sum\limits_y {P(x,y)} }  \\ 
  &=& \sum\limits_x {xP(x)}  \\ 
  &=& \mbox{E}[X] \\ 
\end{eqnarray*}
In like manner, we can prove  $\mbox{E}[\mbox{E}[X|Y]] = \mbox{E}[Y]$

\item The following is for discrete random variables.  A similar
argument holds for continuous random variables.
$\mbox{E}[f(X,Y)|Y] = \mbox{E}[f(X,Y)|Y = y] $ is a function of $y$, i.e., 
$\mbox{E}[f(X,Y)|Y = y] = \sum\limits_x {f(x,y)P(x|Y = y)} = g(y)$
\begin{eqnarray*}
 \mbox{E}[\mbox{E}[f(X,Y)|Y]] &=& \mbox{E}[g(y)] \\ 
  &=& \sum\limits_y {g(y)P(y)}  \\ 
  &=& \sum\limits_y {\sum\limits_x {f(x,y)P(x|Y = y)} P(y)}  \\ 
  &=& \sum\limits_x {\sum\limits_y {f(x,y)P(x,y)} }  \\ 
  &=& \mbox{E}[f(X,Y)] \\ 
\end{eqnarray*}

\item Given that $X_1$ and $X_2$ are continuous random variables, we
know that $\mbox{Pr}[X_1 = x] = 0$ and $\mbox{Pr}[X_2 = x] = 0$ for
any value of $x$.  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 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 = x]$ and $\mbox{Pr}[X_2 = x]$.  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{High-Dimensional Probability}
%Let ${\bf X} = (X_1,X_2,\ldots,X_n)$ be a random vector, where
%the $\{X_i\}$ are independent and identically distributed (i.i.d.)
%continuous random variables with a uniform probability density
%function between 0 and 1:

%$$p(x_i) = \left\{\begin{array}{l}1, \mbox{for $0 \le x_i \le 1$} \\
%                                  0, \mbox{otherwise} 
%		  \end{array}\right.$$

%Each value $\bf x$ of the random vector $\bf X$ can be considered as a
%point in a $n$-dimensional hypercube.  Since the probability density
%function of $\bf X$ is uniform, volume in this $n$-dimensional space
%corresponds directly to probability.  Find an expression for the
%percentage of a $n$-dimensional hypercube's volume located within
%$\epsilon$ of the hypercube's surface.  Plot this percentage as a
%function of $n$ for $1 \le n \le 100000$ and $\epsilon = 0.0001$.
%What do your findings about high-dimensional hypercubes tell you about
%random variables?

%\noindent{\bf Solution}: The volume, $V_n(L)$, of a $n$-dimensional
%hypercube with sides of length $L$ is $$V_n(L) = L^n.$$ Thus, the
%ratio, $R_n(\epsilon)$, of the volumes of a $n$-dimensional hypercube with sides of
%length $1-2\epsilon$ and an $n$-dimensional hypercube with sides of
%length $1$ is 
%$$R_n(\epsilon) = \frac{V_n(1-2\epsilon)}{V_n(1)} =
%\frac{(1-2\epsilon)^n}{1^n} = (1-2\epsilon)^n.$$ Therefore, the
%percentage, $P_n(\epsilon)$, of a $n$-dimensional unit hypercube's
%volume located within $\epsilon$ of the hypercube's surface is 
%$$P_n(\epsilon) = 1 - R_n(\epsilon) = 1 - (1-2\epsilon)^n.$$ See
%Figure \ref{fig:hypercube}.

%\begin{figure}[tb]
%  \centering 
%  \includegraphics[width=\textwidth]{hypercube}
%\caption{The percentage of a $n$-dimensional unit hypercube's volume
%within $\epsilon$ of the surface as a function of $n$.}
%\label{fig:hypercube}
%\end{figure}

%This result shows that randomly choosing a high-dimensional vector in
%which each element is i.i.d. and uniform will with high probability
%result in a point near the edge of the `box' containing all the points
%being considered.  It also hints at the power of considering long
%strings of i.i.d. variables.  See, for example, the asymptotic
%equipartition theorem (AEP), which we may touch on later in class.

%The Python code used to create Figure \ref{fig:hypercube} is

%\lstinputlisting{hypercube.py}

%The Matlab code used to create Figure \ref{fig:hypercube} is

%\lstinputlisting{hypercube.m}

\section{Teatime with Gauss and Bayes}
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=20$, and $\beta = 2.5$.  Plot $p(y)$ and
\hbox{$p(y|x=10.5)$} 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  \\
     &=& \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 Python code produced Figure \ref{fig:GaussBayes}:

%\lstinputlisting{GaussBayes.py}

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

\lstinputlisting{GaussBayes.m}

\begin{figure}[tb]
  \centering 
  \includegraphics[width=\textwidth]{gb.png}
\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}
Let $\Lambda_X = \left[\begin{array}{rr}37 & -15 \\ -15 & 37
\end{array} \right]$.

\begin{enumerate}
\item Verify that $\Lambda_X$ is a valid covariance matrix.
\item Find the eigenvalues and eigenvectors of $\Lambda_X$ by hand.
Show all your work.
\item Write a program to find and verify the eigenvalues and
eigenvectors of $\Lambda_X$.
\item We provide 200 data points sampled from the distribution ${\cal
N}(0,\Lambda_X)$.  Download the dataset from the course website and
plot the data points.  Project the data onto the covariance matrix
eigenvectors and plot the transformed data.  What is the difference
between the two plots?
\end{enumerate}

{\bf Solution}:
\begin{enumerate}

\item The matrix $\Lambda_X$ is a valid covariance matrix if it is
symmetric and positive semi-definite.  Clearly, it is symmetric, since
$\Lambda_X^T = \Lambda_X$.  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 next part of the problem.

\item We can find the eigenvectors and eigenvalues of $\Lambda_X$ by
starting with the definition of an eigenvector.  Namely, an vector ${\bf e}$
is an eigenvector of $\Lambda_X$ if it satisfies
$$\Lambda_X{\bf e} = \lambda{\bf e}$$ for some constant scalar
$\lambda$, which is called the eigenvalue corresponding to ${\bf e}$.
This can be rewritten as
$$(\Lambda_X - \lambda I){\bf e} = {\bf 0}.$$
This is equivalent to 
$$\mbox{det}(\Lambda_X - \lambda I) = 0.$$
Thus, we require that 
$$(37 - \lambda)^2 - 15^2 = 0$$ By inspection, this is true when
$\lambda=52$ and $\lambda=22$, both of which are non-negative, thus
confirming that $\Lambda_X$ is indeed a positive semi-definite matrix.

To find the eigenvectors, we plug the eigenvalues back into the equation above to get
$$(\Lambda_X - 52 I){\bf e} = 
\left[\begin{array}{cc}37-52 & -15 \\ -15 & 37-52\end{array} \right]
\left[\begin{array}{r}a \\ b\end{array} \right] = 
\left[\begin{array}{rr}-15 & -15 \\ -15 & -15\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 = 39$ gives
$$(\Lambda_X - 22 I){\bf e} = 
\left[\begin{array}{cc}37-22 & -15 \\ -15 & 37-22\end{array} \right]
\left[\begin{array}{r}a \\ b\end{array} \right] = 
\left[\begin{array}{rr}15 & -15 \\ -15 & 15\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 $\Lambda_X$:
\lstinputlisting{eigen.m}

\item The following Matlab program generated Figure \ref{fig:covariance}:
\lstinputlisting{covar.m}

\begin{figure}[tb]
  \centering 
  \includegraphics[width=\textwidth]{orig_datadots.png}
  \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}

The second plot in Figure \ref{fig:covariance} shows the data rotated
to align with the eigenvectors of the data's covariance matrix.

\end{enumerate}

\section{Distribution Linearity}

Let $X_1$ and $X_2$ be i.i.d. according to 

$$p(x_i)= \left\{\begin{array}{l}1, \mbox{for $0 \le x_i \le 1$} \\
                                    0, \mbox{otherwise} 
		  \end{array}\right.\mbox{for $i=1,2$}$$

Let $Y = X_1 + X_2$.

\begin{enumerate}
\item Find an expression for $p(y)$.  Plot $p(y)$ for some reasonable
range of $y$.
\item Find an expression for $p(x_1|y)$.  Plot $p(x_1|y)$ as a
function of $x_1$ with $y$ treated as a known parameter for some
reasonable value of $y$ and some reasonable range of $x_1$.
\item Repeat the parts above, this time letting $X_1$ and $X_2$ be
i.i.d. according to ${\cal N}(0,1)$.
\item What was the point of this problem?  Hint: check out the title.
\end{enumerate}

{\bf Solution}: 
\begin{enumerate}
\item From basic probability theory, we know that the probability
density function of the sum of two independent random variables is the
convolution of the two probability density functions.  So,
\begin{eqnarray*}
p_y(y) &=& (p_{x_1}*p_{x_2})(y) \\
       &=& \int_{-\infty}^{\infty}\!p_{x_1}(x)\,p_{x_2}(y-x)\,dx \\
       &=& \int_{0}^{1}\!1\,p_{x_2}(y-x)\,dx \\
       &=& \int_{0}^{1}\!p_{x_2}(y-x)\,dx \\
       &=& \int_{0}^{1}\left\{\begin{array}{ll} 1 & \mbox{for $0 \le y-x \le 1$} \\
                                         0 & \mbox{otherwise} \end{array} \right\} dx \\
       &=& \int_{0}^{1}\left\{\begin{array}{ll} 1 & \mbox{for $y-1 \le x \le y$} \\
                                         0 & \mbox{otherwise} \end{array} \right\} dx \\
       &=& \int_{\mbox{max}(0,y-1)}^{\mbox{min}(1,y)}\!1\,dx \\
       &=& \mbox{max}\{0, \mbox{min}(1,y) - \mbox{max}(0,y-1)\} \\
       &=& \left\{\begin{array}{ll}
             0 & \mbox{for $y \le 0$} \\
             y & \mbox{for $0 \le y \le 1$} \\
             2-y & \mbox{for $1 \le y \le 2$} \\
             0 & \mbox{for $y \ge 2$}
           \end{array}\right.
\end{eqnarray*}

This p.d.f. is shown in Figure \ref{fig:UniformConv}, which was produced using the following Python program:
\lstinputlisting{UniformConv.py}

This p.d.f. is shown in Figure \ref{fig:UniformConv}, which was produced using the following Matlab program:
\lstinputlisting{UniformConv.m}

\begin{figure}[tb]
  \centering 
  \includegraphics[width=\textwidth]{UniformConv}
\caption{The probability density function of the sum of two
independent uniform random variables.}
\label{fig:UniformConv}
\end{figure}

\item Using Bayes' Rule, we have
\begin{eqnarray*}
p_{x_1|y}(x_1|y) &=& \frac{p_{x_1,y}(x_1,y)}{p_y(y)} \\
         &=& \frac{p_{y|x_1}(y|x_1)\,p_{x_1}(x_1)}{p_y(y)}
\end{eqnarray*}
We already know $p_y(y)$ and $p_{x_1}(x_1)$.  Finding
$p_{y|x_1}(y|x_1)$ is a matter of realizing that $y = x_1 + x_2$
implies that, given $x_1$, $y$ is simply $x_2$ offset by a constant.
Thus, $$p_{y|x_1}(y|x_1) = p_{x_2}(y - x_1)$$
and
\begin{eqnarray*}
p_{x_1|y}(x_1|y) &=& \frac{p_{x_2}(y - x_1)\,p_{x_1}(x_1)}{p_y(y)} \\
       &=& \left\{\begin{array}{ll}
             \frac{1}{y} & \mbox{for $0 \le y \le 1$ and $0 \le x_1 \le y$}\\
             \frac{1}{2-y} & \mbox{for $1 \le y \le 2$ and $y-1 \le x_1 \le 1$}\\
             0 & \mbox{otherwise}
           \end{array}\right.
\end{eqnarray*}

See Figure \ref{fig:X1givenY}, which was produced by the following Python program:
\lstinputlisting{X1givenY.py}

See Figure \ref{fig:X1givenY}, which was produced by the following Matlab program:
\lstinputlisting{X1givenY.m}

\begin{figure}[tb]
  \centering 
  \includegraphics[width=\textwidth]{X1givenY}
\caption{The probability density function of $x_1$ given certain
values of $y$, where $y = x_1 + x_2$ and $x_1$ and $x_2$ are
i.i.d. uniform random variables.}
\label{fig:X1givenY}
\end{figure}

\item Repeating the above using normal distributions, we get
\begin{eqnarray*}
p_y(y) &=& (p_{x_1}*p_{x_2})(y) \\
       &=& \int_{-\infty}^{\infty}\!p_{x_1}(x)\,p_{x_2}(y-x)\,dx \\
       &=& \int_{-\infty}^{\infty}
           \left(\frac{1}{\sqrt{2\pi}}e^{\frac{-x^2}{2}}\right)
           \left(\frac{1}{\sqrt{2\pi}}e^{\frac{-(y-x)^2}{2}}\right) \, dx \\
       &=& \int_{-\infty}^{\infty}
           \frac{1}{2\pi}e^{\frac{-(2x^2 - 2xy + y^2)}{2}}\,dx \\
       &=& \int_{-\infty}^{\infty}
           \frac{1}{2\pi}e^{\frac{-(x^2 - xy + \frac{y^2}{4} - \frac{y^2}{4} + \frac{y^2}{2})}{2\frac{1}{2}}}\,dx \\
       &=& \int_{-\infty}^{\infty}
           \frac{1}{2\pi}e^{-\left((x - \frac{y}{2})^2 - \frac{y^2}{4} + \frac{y^2}{2}\right)}\,dx \\
       &=& \frac{1}{\sqrt{4\pi}}e^{-\frac{y^2}{4}}\,
           \int_{-\infty}^{\infty}
           \frac{1}{\sqrt{\pi}}e^{-(x - \frac{y}{2})^2}\,dx \\
       &=& \frac{1}{\sqrt{4\pi}}e^{-\frac{y^2}{4}} \\
       &=& {\cal N}(0,2)
\end{eqnarray*}
Similarly,
\begin{eqnarray*}
p_{x_1|y}(x_1|y) &=& \frac{p_{x_2}(y - x_1)\,p_{x_1}(x_1)}{p_y(y)} \\
                 &=& \frac{\left(\frac{1}{\sqrt{2\pi}}e^{\frac{-(y - x_1)^2}{2}}\right)
                     \left(\frac{1}{\sqrt{2\pi}}e^{\frac{-x_1^2}{2}}\right)}
                     {\left(\frac{1}{\sqrt{4\pi}}e^{\frac{-y^2}{4}}\right)} \\
                 &=& \frac{1}{\sqrt{\pi}}e^{-\frac{y^2 - 4x_1y + 4x_1^2}{4}} \\
                 &=& \frac{1}{\sqrt{\pi}}e^{-(x_1 - \frac{y}{2})^2} \\
                 &=& {\cal N}\left(\frac{y}{2},\frac{1}{2}\right)
\end{eqnarray*}

See Figure \ref{fig:normalX1givenY}, which was produced by the following Python program:
\lstinputlisting{normalX1givenY.py}

See Figure \ref{fig:normalX1givenY}, which was produced by the following Matlab program:
\lstinputlisting{normalX1givenY.m}

\begin{figure}[tb]
  \centering 
  \includegraphics[width=\textwidth]{normalX1givenY}
\caption{The probability density function of $y$ and the probability
density function for $x_1$ given $y=1.6$, where $y = x_1 + x_2$ and
$x_1$ and $x_2$ are i.i.d. ${\cal N}(0,1)$ random variables.}
\label{fig:normalX1givenY}
\end{figure}

\item The point of this problem is to show that probability density
functions are in general not closed under linear combinations of
i.i.d. random variables.  That is, given two i.i.d. random variables
$x_1$ and $x_2$ with distribution of type $A$, the random variable $y
= x_1 + x_2$ does not in general have a distribution of type $A$.
Gaussian (a.k.a normal) distributions are an exception.  In fact,
Gaussians are the only non-trivial family of functions that are both
closed and linear under convolution (and therefore under addition of
i.i.d. random variables): 
$${\cal N}(\mu_1,\sigma^2_1)*{\cal N}(\mu_2,\sigma^2_2) = {\cal N}(\mu_1+\mu_2,\sigma^2_1+\sigma^2_2)$$

\end{enumerate}



\section{Probabilistic Modeling}
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 (variable $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{tabular}{c | c}
$x$ & $y$\\
\hline
1 & 0\\
1 & 1\\
0 & 0\\
1 & 1\\
1 & 0\\
0 & 1\\
0 & 0 \\
\end{tabular}

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{3}{7}$,  $\widehat\theta _2  = \frac{4}{7}$.
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{4}{7}} \right)^8 \left( {\frac{3}{7}} \right)^6  \approx 7.044 \times 10^{ - 5} 
\]
\end{enumerate}


\section{Monty Hall}
To get credit for this problem, you must not only write your own
correct solution, but also write a computer simulation (in either Matlab
or Python) 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
the 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?  Explain your answer.
\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}

{\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)


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 Monty Hall problem:

\lstinputlisting{montysimulation.txt}

Here is a Python program that generates the Monty Hall simulation output above:

\lstinputlisting{montyhall.py}

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

\lstinputlisting{montyhall.m}



\end{itemize}

\end{document}
