Problem Set 1 Solutions

MAS.622J/1.126J: Pattern Recognition and Analysis

Originally Due Monday, 15 September 2008

Problem 1: Why?

a. 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.

b. In the same way, describe an application of pattern recognition you would be interested in pursuing for fun in your life outside of work.

Solution: Refer to examples discussed in lecture.

Problem 2: Probability Warm-Up

Let X and Y be random variables. Let μ_{X} ≡ E[X] denote the expected value of X and σ_{X}^{2} ≡ E[(X − μ_{X})^{2}] = E[X^{2}] − μ_{X}^{2} denote the variance of X. a and b are constant values. Use excruciating detail to answer the following:

a. Show E[aX + bY] = aE[X] + bE[Y].

b. Show that independent implies uncorrelated.

c. Show that uncorrelated does not imply independent.

d. Let Z = aX + bY. Show that if X and Y are uncorrelated, then σ_{Z}^{2} = a^{2}σ_{X}^{2} + b^{2}σ_{Y}^{2} .

e. Let X_{i} (i = 1, ..., n) be random variables independently drawn from the same probability distribution with mean μ_{X} and variance σ_{X}^{2}. For the

sample mean \overline{X} = \frac{1}{n} ∑_{i=1}^{n} X_{i}, show the following: (i) E[\overline{X}] = μ_{X}. (ii)

Var[\overline{X}] (variance of the sample mean) = σ_{X}^{2}/n. Note that this is different

from the sample variance s_{n}^{2} = \frac{1}{n} ∑_{i=1}^{n} (X_{i} − \overline{X})^{2}.

f. The conditional expected value 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 E[X|Y = y] = ∑_{x} xP (x|Y = y) is just g(Y). Show E[X] = E[E[X|Y]] and E[Y] =

E[E[Y |X]].

g. For a real value function f and discrete random variables X and Y, E[f(X, Y)] = ∑_{x} ∑_{y} f(x, y)P (x, y). Show E[f(X, Y)] = E[E[f(X, Y)|Y]].

h. Let X_{1} and X_{2} be independent and identically distributed continuous random variables. Can Pr[X_{1} ≤ X_{2}] be calculated? If so, find its value. If not, explain.

i. Let X_{1} and X_{2} be independent and identically distributed discrete random variables. Can Pr[X_{1} ≤ X_{2}] be calculated? If so, find its value. If not, explain.

Solution:

a. The following is for continuous random variables. A similar argument holds for discrete random variables.

E[aX + bY] = ∫ ∫ (ax + by) p(x, y) dx dy
= a ∫ ∫ x p(x, y) dx dy + b ∫ ∫ y p(x, y) dx dy
= a ∫ x p(x) dx + b ∫ y p(y) dy
= aE[X] + bE[Y]

b. For uncorrelated σ_{xy} = 0. Let X and Y be independent continuous random variables (a similar argument holds for discrete random variables). Then,

E[XY] = ∫ ∫ xy p(x, y) dx dy
= ∫ ∫ xy p(x) p(y) dx dy
= ∫ x p(x) dx ∫ y p(y) dy
= E[X] E[Y]

c. 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

p_{x}[x = 0] = 0.5
p_{x}[x = 1] = 0.5

and the probability mass function of Y conditioned on X be

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.

Given the above, and the fact that p_{x,y}[x, y] = p_{y|x}[y|x] p_{x}[x], we get

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.

However, the product of the marginals is given by

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.

Thus, we see that p_{x,y}[x, y] ≠ p_{x}[x] p_{y}[y] and X and Y are not independent. independent. However, since XY is identically zero, we also get

cov(X, Y) = σ_{XY}^{2} = E[(X − μ_{X})(Y − μ_{Y})]
= E[XY] − μ_{X}μ_{Y}
= E[0] − (0.5)(0)
= 0 − 0
= 0.

Therefore, X and Y are uncorrelated but not independent.

d. Given that Z = aX + bY and that X and Y are uncorrelated, we have

σ_{Z}^{2} = E[(Z − μ_{Z})^{2}]

= E[Z^{2}] − μ_{Z}^{2}
= E[(aX + bY)^{2}] − (aμ_{X} + bμ_{Y})^{2}
= E[a^{2}X^{2} + 2abXY + b^{2}Y^{2}] − (a^{2}μ_{X}^{2} + 2abμ_{X}μ_{Y} + b^{2}μ_{Y}^{2} )
= a^{2}E[X^{2}] + 2abE[XY] + b^{2}E[Y^{2}] − a^{2}μ_{X}^{2} − 2abμ_{X}μ_{Y} − b^{2}μ_{Y}^{2}
= a^{2}(E[X^{2}] − μ_{X}^{2}) + 2ab(E[XY] − μ_{X}μ_{Y}) + b^{2}(E[Y^{2}] − μ_{Y}^{2} )
= a^{2}σ_{X}^{2} + 2abσ_{XY}^{2} + b^{2}σ_{Y}^{2}
= a^{2}σ_{X}^{2} + b^{2}σ_{Y}^{2} ,

where only the last equality depends on X and Y being uncorrelated.

e. Using the result of (a) and the fact E[X_{i}] = μ_{X},

E[X̄] = E[\frac{1}{n} ∑_{i=1}^{n} X_{i}] = \frac{1}{n} ∑_{i=1}^{n} E[X_{i}] = \frac{1}{n} n μ_{X} = μ_{X}

Also, using the result of (d) and the fact Var[X_{i}] = σ_{X}^{2}

Var[X̄] = Var[\frac{1}{n} ∑_{i=1}^{n} X_{i}] = \frac{1}{n^{2}} ∑_{i=1}^{n} Var[X_{i}] = \frac{1}{n^{2}} n σ_{X}^{2} = σ_{X}^{2}/n

f. The following is for discrete random variables. A similar argument holds for continuous random variables. E[X|Y] = E[X|Y = y] is a function of

y, i.e., E[X|Y = y] = ∑_{x} xP (x|Y = y) = g(y)

E[E[X|Y]] = E[g(y)]
= ∑_{y} g(y)P (y)
= ∑_{y} ∑_{x} xP (x|Y = y)P (y)
= ∑_{y} ∑_{x} xP (x, y)
= ∑_{x} x ∑_{y} P (x, y)
= ∑_{x} xP (x)
= E[X]

In like manner, we can prove E[E[X|Y]] = E[Y]

g. The following is for discrete random variables. A similar argument holds for continuous random variables. E[f(X, Y)|Y] = E[f(X, Y)|Y = y] is a function of y, i.e., E[f(X, Y)|Y = y] = ∑_{x} f(x, y)P (x|Y = y) = g(y)

E[E[f(X, Y)|Y]] = E[g(y)]
= ∑_{y} g(y)P (y)
= ∑_{y} ∑_{x} f(x, y)P (x|Y = y)P (y)
= ∑_{x} ∑_{y} f(x, y)P (x, y)
= E[f(X, Y)]

h. Given that X_{1} and X_{2} are continuous random variables, we know that Pr[X_{1} = x] = 0 and Pr[X_{2} = x] = 0 for any value of x. Thus,

Pr[X_{1} ≤ X_{2}] = 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

Pr[X_{1} < X_{2}] = Pr[X_{2} < X_{1}].

However, since probabilities must sum to one, we have

Pr[X_{1} < X_{2}] + Pr[X_{2} < X_{1}] = 1.

Thus,

Pr[X_{1} ≤ X_{2}] = \frac{1}{2}.

i. 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 Pr[X_{1} = x] and Pr[X_{2} = x]. Thus, the argument we used above fails. In general, it is not possible to find Pr[X_{1} ≤ X_{2}] without knowledge of the distributions of both X_{1} and X_{2}.

Problem 3: Teatime with Gauss and Bayes

Let p(x, y) = \frac{1}{2παβ} e^{−(\frac{(y−μ)^{2}}{2α^{2}} + \frac{(x−y)^{2}}{2β^{2}} )}.

a. Find p(x), p(y), p(x|y), and p(y|x). In addition, give a brief description of each of these distributions.

b. Let μ = 0, α = 20, and β = 2.5. Plot p(y) and p(y|x = 10.5) for a reasonable range of y. What is the difference between these two distributions?

Solution:

a. To find p(y), simply factor p(x, y) and then integrate over x:

p(y) = ∫_{−∞}^{∞}p(x, y) dx

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(\frac{(y−μ)^{2}}{2α^{2}} + \frac{(x−y)^{2}}{2β^{2}} )} dx
= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{− \frac{(y−μ)^{2}}{2α^{2}}} e^{− \frac{(x−y)^{2}}{2β^{2}}} dx
= \frac{1}{\sqrt{2πα^{2}}} e^{− \frac{(y−μ)^{2}}{2α^{2}}} ∫_{−∞}^{∞} \frac{1}{\sqrt{2πβ^{2}}} e^{− \frac{(x−y)^{2}}{2β^{2}}} dx
= \frac{1}{\sqrt{2πα^{2}}} e^{− \frac{(y−μ)^{2}}{2α^{2}}}
= N (μ, α^{2})

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):

p(x|y) = \frac{p(x, y)}{p(y)}
= \frac{1}{\sqrt{2πβ^{2}}} e^{− \frac{(x−y)^{2}}{2β^{2}}}
= N (y, β^{2})

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.

p(x) = ∫_{−∞}^{∞}p(x, y) dy

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(\frac{(y−μ)^{2}}{2α^{2}} + \frac{(x−y)^{2}}{2β^{2}} )} dy

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(\frac{β^{2}(y−μ)^{2}+α^{2}(x−y)^{2}}{2α^{2}β^{2}} )} dy

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(\frac{β^{2}y^{2}−2β^{2}μy+β^{2}μ^{2}+α^{2}x^{2}−2α^{2}xy+α^{2}y^{2}}{2α^{2}β^{2}} )} dy

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(\frac{(α^{2}+β^{2})y^{2}−2(α^{2}x+β^{2}μ)y+(β^{2}μ^{2}+α^{2}x^{2})}{2α^{2}β^{2}})} dy

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(\frac{y^{2}−2 \frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}} y+ \frac{β^{2}μ^{2}+α^{2}x^{2}}{α^{2}+β^{2}}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}} )} dy

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(^{\frac{y^{2}−2 \frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}} y+(\frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}−(\frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}+ \frac{β^{2}μ^{2}+α^{2}x^{2}}{α^{2}+β^{2}}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} )} dy

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(^{\frac{(y− \frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}−(\frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}+ \frac{β^{2}μ^{2}+α^{2}x^{2}}{α^{2}+β^{2}}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} )} dy

= ∫_{−∞}^{∞} \frac{1}{2παβ} e^{−(^{\frac{(y− \frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} )} e^{−(^{\frac{\frac{β^{2}μ^{2}+α^{2}x^{2}}{α^{2}+β^{2}} −(\frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} )} dy

= \frac{1}{2παβ} \sqrt{2π \frac{α^{2}β^{2}}{α^{2} + β^{2}}} e^{−(^{\frac{\frac{β^{2}μ^{2}+α^{2}x^{2}}{α^{2}+β^{2}} −(\frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} )} ∫_{−∞}^{∞} \frac{1}{\sqrt{2π \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} e^{−(^{\frac{(y− \frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} )} dy

= \frac{1}{\sqrt{2π(α^{2} + β^{2})}}e^{−(^{\frac{\frac{β^{2}μ^{2}+α^{2}x^{2}}{α^{2}+β^{2}} −(\frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} )}

= \frac{1}{\sqrt{2π(α^{2} + β^{2})}}e^{−(\frac{(α^{2}+β^{2})(β^{2}μ^{2}+α^{2}x^{2})−(α^{2}x+β^{2}μ)^{2}}{2α^{2}β^{2}(α^{2}+β^{2})} )}

= \frac{1}{\sqrt{2π(α^{2} + β^{2})}}e^{−(\frac{α^{2}β^{2}μ^{2}+α^{4}x^{2}+β^{4}μ^{2}+α^{2}β^{2}x^{2}−α^{4}x^{2}−2α^{2}β^{2}μx−β^{4}μ^{2}}{2α^{2}β^{2}(α^{2}+β^{2})} )}

= \frac{1}{\sqrt{2π(α^{2} + β^{2})}}e^{−(\frac{α^{2}β^{2}x^{2}−2α^{2}β^{2}μx+α^{2}β^{2}μ^{2}}{2α^{2}β^{2}(α^{2}+β^{2})} )}

= \frac{1}{\sqrt{2π(α^{2} + β^{2})}}e^{−(\frac{(x−μ)^{2}}{2(α^{2}+β^{2})})}
= N (μ, α^{2} + β^{2})

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):

p(y|x) = \frac{p(x, y)}{p(x)}

= \frac{1}{\sqrt{2π \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} e^{−(^{\frac{(y− \frac{α^{2}x+β^{2}μ}{α^{2}+β^{2}})^{2}}{2 \frac{α^{2}β^{2}}{α^{2}+β^{2}}}} )}
= N (\frac{α^{2}x + β^{2}μ}{α^{2} + β^{2}}, \frac{α^{2}β^{2}}{α^{2} + β^{2}} )

Note that all the above distibutions are Gaussian. b. The following Matlab code produced Figure 1:

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 g i v e n x = ( 1 . 0 / s q r t (2∗ pi ∗ var ) ) ∗ exp ( −((y−mean ) . ˆ 2 ) / ( 2 ∗ var ) )
var2 = aˆ 2 + bˆ 2
p y = ( 1 . 0 / s q r t (2∗ pi ∗ var2 ) ) ∗ exp ( −((y−m) . ˆ 2 ) / ( 2 ∗ var2 ) )
hold o f f
p l o t (y , p y given x , ’ b ’ )
hold on
p l o t (y , p y , ’ r ’ )
legend ( ’ p( y | x ) ’ , ’p( y ) ’ )
sy = s i z e ( y )
a x i s ( [ y ( 1 ) , y ( sy ( 2 ) ) , 0 , 0 . 2 ] )
x l a b e l ( ’ y ’ )
t e x t ( −70 ,0.14 , ’\mu=0 ’)
t e x t ( −70 ,0.12 , ’\ alpha =15 ’)
t e x t ( −70 ,0.1 , ’\ beta =3 ’)

Problem 4: Covariance Matrix

Let Λ_{X} = [ \begin{matrix}37 & −15 \\ −15 & 37 \end{matrix} ].

a. Verify that Λ_{X} is a valid covariance matrix. b. Find the eigenvalues and eigenvectors of Λ_{X} by hand. Show all your work. c. Write a program to find and verify the eigenvalues and eigenvectors of Λ_{X}.

Figure 1: 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.

d. We provide 200 data points sampled from the distribution N (0, Λ_{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?

Solution:

a. The matrix Λ_{X} is a valid covariance matrix if it is symmetric and positive semi-definite. Clearly, it is symmetric, since Λ_{X}^{T} = Λ_{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.

b. We can find the eigenvectors and eigenvalues of Λ_{X} by starting with the definition of an eigenvector. Namely, an vector e is an eigenvector of Λ_{X} if it satisfies

Λ_{X}e = λe

for some constant scalar λ, which is called the eigenvalue corresponding to e. This can be rewritten as

(Λ_{X} − λI)e = 0.

This is equivalent to

det(Λ_{X} − λI) = 0.

Thus, we require that

(37 − λ)^{2} − 15^{2} = 0

By inspection, this is true when λ = 52 and λ = 22, both of which are non-negative, thus confirming that Λ_{X} is indeed a positive semi-definite matrix.

To find the eigenvectors, we plug the eigenvalues back into the equation above to get

(Λ_{X}−52I)e = [ \begin{matrix}37 − 52 & −15 \\ −15 & 37 − 52 \end{matrix} ] [ \begin{matrix}a \\ b \end{matrix} ] = [ \begin{matrix}−15 & −15 \\ −15 & −15 \end{matrix} ] [ \begin{matrix}a \\ b \end{matrix} ] = [ \begin{matrix}0 \\ 0 \end{matrix} ] ,

which gives a = −b. Normalized, this results in the eigenvector

e_{1} = [ −\frac{1}{\sqrt{2}}^{\frac{1}{\sqrt{2}}} ] .

Similarly, λ = 39 gives

(Λ_{X}−22I)e = [ \begin{matrix}37 − 22 & −15 \\ −15 & 37 − 22 \end{matrix} ] [ \begin{matrix}a \\ b \end{matrix} ] = [ \begin{matrix}15 & −15 \\ −15 & 15 \end{matrix} ] [ \begin{matrix}a \\ b \end{matrix} ] = [ \begin{matrix}0 \\ 0 \end{matrix} ] ,

which gives a = b. Normalized, this results in the eigenvector

e_{1} = [ \frac{1}{\sqrt{2}} ; \frac{1}{\sqrt{2}} ] .

c. The following Matlab program prints out the eigenvectors and eigenvalues of Λ_{X}:

A = [ 3 7 −15; −15 3 7 ]
[V,D] = e i g (A)

d. The following Matlab program generated Figure 2:

t t = load ( ’ ps1 . txt ’ )
% o r i g i n a l c o r r e l a t i o n
t t c o r r = c o r r c o e f ( t t ( 1 , : ) , t t ( 2 , : ) )

f i g u r e

p l o t ( t t ( 1 , : ) , t t ( 2 , : ) , ’ . ’ )
x l a b e l ( ’ x ’ )
y l a b e l ( ’ y ’ )
% e i g e n v e c t o r
% A = [ 3 7 −15; −15 3 7 ]
% [V,D] = e i g (A)
V =[ 0.7071  0.7071
−0.7071  0 . 7 0 7 1 ]
zz = V’ ∗ t t ; % a x i s transformation
% c o r r e l a t i o n a f t e r transformation
z z c o r r = c o r r c o e f ( zz ( 1 , : ) , zz ( 2 , : ) )

f i g u r e

p l o t ( zz ( 1 , : ) , zz ( 2 , : ) , ’ . ’ )
x l a b e l ( ’ f i r s t eigenvector ’ )
y l a b e l ( ’ second eigenvector ’ )

The second plot in Figure 2 shows the data rotated to align with the eigenvectors of the data’s covariance matrix.

Problem 5: Distribution Linearity

Let X_{1} and X_{2} be i.i.d. according to

p(x_{i}) = \{ \begin{matrix}1, for 0 ≤ x_{i} ≤ 1 \\ 0, otherwise \end{matrix} for i = 1, 2

Let Y = X_{1} + X_{2}.

a. Find an expression for p(y). Plot p(y) for some reasonable range of y.

b. 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}.

c. Repeat the parts above, this time letting X_{1} and X_{2} be i.i.d. according to N (0, 1).

d. What was the point of this problem? Hint: check out the title.

Solution:

a. 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,

p_{y}(y) = (p_{x_{1}} ∗ p_{x_{2}})(y)
= ∫_{−∞}^{∞}p_{x_{1}}(x) p_{x_{2}}(y − x) dx
= ∫_{0}^{1}1 p_{x_{2}}(y − x) dx
= ∫_{0}^{1}p_{x_{2}}(y − x) dx
= ∫_{0}^{1} \{ \begin{matrix}1 & for 0 ≤ y − x ≤ 1 \\ 0 & otherwise \end{matrix} \} dx
= ∫_{0}^{1} \{ \begin{matrix}1 & for y − 1 ≤ x ≤ y \\ 0 & otherwise \end{matrix} \} dx

= ∫_{max(0,y−1)}^{min(1,y)}1 dx
= max{0, min(1, y) − max(0, y − 1)}

= \{ \begin{matrix}0 & for y ≤ 0 \\ y & for 0 ≤ y ≤ 1 \\ 2 − y & for 1 ≤ y ≤ 2 \\ 0 & for y ≥ 2 \end{matrix}

This p.d.f. is shown in Figure 3, which was produced using the following Python program:

from m a t p l o t l i b . numerix import ∗
from numarray import ∗
from pylab import plot , subplot , legend , axis , xlabel , ylabel , text , show
Error . setMode ( a l l=None , overflow =’warn ’ , underflow =’ ignore ’ , d i v i d e b y z e r o =’warn ’ , i n v a l i d =’warn ’ )
p l o t ( [ − 1 , 0 , 1 , 2 , 3 ] , [ 0 , 0 , 1 , 0 , 0 ] )
a x i s ([ −1 ,3 , −0.5 ,2])
x l a b e l ( ’ y ’ )
y l a b e l ( ’ p( y ) ’ )
show ( )

This p.d.f. is shown in Figure 3, which was produced using the following Matlab program:

hold o f f

subplot (111)
p l o t ( [ − 1 , 0 , 1 , 2 , 3 ] , [ 0 , 0 , 1 , 0 , 0 ] )
hold on
a x i s ([ −1 ,3 , −0.5 ,2])
x l a b e l ( ’ y ’ )
y l a b e l ( ’ p( y ) ’ )

b. Using Bayes’ Rule, we have

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

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

p_{x_{1}|y}(x_{1}|y) = \frac{p_{x_{2}}(y − x_{1}) p_{x_{1}}(x_{1})}{p_{y}(y)}

= \{ \begin{matrix}\frac{1}{y} & for 0 ≤ y ≤ 1 and 0 ≤ x_{1} ≤ y \\ \frac{1}{2−y} & for 1 ≤ y ≤ 2 and y − 1 ≤ x_{1} ≤ 1 \\ 0 & otherwise \end{matrix}

See Figure 4, which was produced by the following Python program:

from m a t p l o t l i b . numerix import ∗
from numarray import ∗
from pylab import plot , subplot , legend , axis , xlabel , ylabel , text , show
Error . setMode ( a l l=None , overflow =’warn ’ , underflow =’ ignore ’ , d i v i d e b y z e r o =’warn ’ , i n v a l i d =’warn ’ )
subplot (211)
f o r y in arange ( 0 . 1 , 1 . 1 , s t r i d e =0.1) :
x = array ( [ 0 , 0 , y , y , 1 , 1 ] )
Px given y = array ( [ 0 , 1 . 0 / y , 1 . 0 / y , 0 , 0 , 0 ] )
p l o t (x , Px given y )
x l a b e l ( r ’ $x 1$ ’ )
y l a b e l ( r ’ $p ( x 1 \ given \ 0<y<1)$ ’ )
a x i s ( [ − 0 . 1 , 1 . 1 , 0 , 1 2 ] )
subplot (212)
f o r y in arange ( 1 . 0 , 2 . 0 , s t r i d e =0.1) :
x = array ( [ 0 , 0 , y−1,y −1 ,1 ,1])

Px given y = array ([0 ,0 ,0 ,1.0/(2 − y ) ,1.0/(2 − y ) , 0 ] )
p l o t (x , Px given y )
x l a b e l ( ’ $x 1$ ’ )
y l a b e l ( ’ $p ( x 1 \ given \ 1<y<2)$ ’ )
a x i s ( [ − 0 . 1 , 1 . 1 , 0 , 1 2 ] )
show ( )

See Figure 4, which was produced by the following Matlab program:

hold o f f
subplot (211)
hold on
f o r y = 0 . 1 : 0 . 1 : 1
x = [ 0 , 0 , y , y , 1 ]
Px given y = [ 0 , 1 . 0 / y , 1 . 0 / y , 0 , 0 ]
p l o t (x , Px given y , ’ b ’ )

end

x l a b e l ( ’ x 1 ’ )
y l a b e l ( ’ p( x 1 | 0<y <1) ’)
a x i s ( [ − 0 . 1 , 1 . 1 , 0 , 1 2 ] )
subplot (212)
hold on
f o r y = 1 . 0 : 0 . 1 : 1 . 9
x = [ 0 , y−1,y −1 ,1 ,1]
Px given y = [0 ,0 ,1.0/(2 − y ) ,1.0/(2 − y ) , 0 ]
p l o t (x , Px given y , ’ r ’ )

end

x l a b e l ( ’ x 1 ’ )
y l a b e l ( ’ p( x 1 | 1<y <2) ’)
a x i s ( [ − 0 . 1 , 1 . 1 , 0 , 1 2 ] )

c. Repeating the above using normal distributions, we get

p_{y}(y) = (p_{x_{1}} ∗ p_{x_{2}})(y)
= ∫_{−∞}^{∞}p_{x_{1}}(x) p_{x_{2}}(y − x) dx
= ∫_{−∞}^{∞} (\frac{1}{\sqrt{2π}} e^{\frac{−x^{2}}{2}} ) (\frac{1}{\sqrt{2π}} e^{\frac{−(y−x)^{2}}{2}} ) dx
= ∫_{−∞}^{∞} \frac{1}{2π} e^{\frac{−(2x^{2}−2xy+y^{2})}{2}} dx

= ∫_{−∞}^{∞} \frac{1}{2π} e^{\frac{−(x^{2}−xy+ \frac{y^{2}}{4} − \frac{y^{2}}{4} + \frac{y^{2}}{2} )}{2 \frac{1}{2}}} dx

= ∫_{−∞}^{∞} \frac{1}{2π} e^{−((x− \frac{y}{2} )^{2}− \frac{y^{2}}{4} + \frac{y^{2}}{2} )} dx
= \frac{1}{\sqrt{4π}} e^{− \frac{y^{2}}{4}} ∫_{−∞}^{∞} \frac{1}{\sqrt{π}} e^{−(x− \frac{y}{2} )^{2}} dx
= \frac{1}{\sqrt{4π}} e^{− \frac{y^{2}}{4}}
= N (0, 2)

Similarly,

p_{x_{1}|y}(x_{1}|y) = \frac{p_{x_{2}}(y − x_{1}) p_{x_{1}}(x_{1})}{p_{y}(y)}

= \frac{(\frac{1}{\sqrt{2π}} e^{\frac{−(y−x_{1})^{2}}{2}} ) (\frac{1}{\sqrt{2π}} e^{\frac{−x^{2}_{1}}{2}} )}{(\frac{1}{\sqrt{4π}} e^{\frac{−y^{2}}{4}} )}

= \frac{1}{\sqrt{π}} e^{−^{\frac{y^{2}−4x_{1}y+4x^{2}_{1}}{4}}}
= \frac{1}{\sqrt{π}} e^{−(x_{1}− \frac{y}{2} )^{2}}
= N (\frac{y}{2} , \frac{1}{2})

See Figure 5, which was produced by the following Python program:

from m a t p l o t l i b . numerix import ∗
from numarray import ∗
from pylab import plot , subplot , legend , axis , xlabel , ylabel , text , show
Error . setMode ( a l l=None , overflow =’warn ’ , underflow =’ ignore ’ , d i v i d e b y z e r o =’warn ’ , i n v a l i d =’warn ’ )
import LinearAlgebra as l a
subplot (211)
y = arange ( −5 ,5 ,0.01)
p = ( 1 . 0 / s q r t (4∗ pi ) ) ∗ ( e ∗∗(−(y ∗∗2)/4))
p l o t (y , p)
x l a b e l ( r ’ $y$ ’ )
y l a b e l ( r ’ $p ( y ) $ ’ )
a x i s ([ −5 ,5 , −0.2 ,1.0])
subplot (212)
y = 1 . 6
x = arange ( −5 ,5 ,0.01)

p = ( 1 . 0 / s q r t ( pi ) ) ∗ ( e ∗∗( −((x−y / 2 ) ∗ ∗ 2 ) ) )
p l o t (x , p)
x l a b e l ( r ’ $x 1$ ’ )
y l a b e l ( r ’ $p ( x 1 \ given \ y=1.6)$ ’ )
a x i s ([ −5 ,5 , −0.2 ,1.0])
show ( )

See Figure 5, which was produced by the following Matlab program:

subplot ( 2 1 1 ) ;
y = − 5 : 0 . 0 1 : 5 ;
p = ( 1 . 0 / s q r t (4∗ pi ) ) ∗ ( exp(−(y . ˆ 2 ) / 4 ) ) ;
p l o t (y , p ) ;
x l a b e l ( ’ y ’ ) ;
y l a b e l ( ’ p( y ) ’ ) ;
a x i s ( [ − 5 , 5 , − 0 . 2 , 1 . 0 ] ) ;
subplot ( 2 1 2 ) ;
y = 1 . 6 ;
x = − 5 : 0 . 0 1 : 5 ;
p = ( 1 . 0 / s q r t ( pi ) ) ∗ ( exp ( −((x−y / 2 ) . ˆ 2 ) ) ) ;
p l o t (x , p ) ;
x l a b e l ( ’ x 1 ’ ) ;
y l a b e l ( ’ p( x 1 | y = 1 . 6 ) ’ ) ;
a x i s ( [ − 5 , 5 , − 0 . 2 , 1 . 0 ] ) ;

d. 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):

N (μ_{1}, σ_{1}^{2}) ∗ N (μ_{2}, σ_{2}^{2}) = N (μ_{1} + μ_{2}, σ_{1}^{2} + σ_{2}^{2})

Problem 6: Probabilistic Modeling

Let x ∈ {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 θ_{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 θ_{2}. a. Represent the joint probability distribution P (x, y|θ) for all x, y (a 2x2

matrix) as a function of the parameters θ = (θ_{1}, θ_{2}).

b. The Maximum Likelihood estimation criterion for the parameter θ is defined as:

θ̂_{ML} = arg max_{θ} L(t_{1}, ..., t_{n}; θ) = arg max_{θ} ∏_{i=1}^{n} p(t_{i}|θ)

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}; θ) =

∏_{i=1}^{n} p(t_{i}|θ). Likelihood is viewed as a function of the parameters, which de-

pends on the data. Since the above expression can be technically challenging, we maximize the log-likelihood log L(t_{1}, ..., t_{n}; θ) instead of likelihood. Note that any monotonically increasing function (i.e., log function) of the likelihood has the same maxima. Thus,

θ̂_{ML} = arg max_{θ} log L(t_{1}, ..., t_{n}; θ) = arg max_{θ} ∑_{i=1}^{n} log p(t_{i}|θ)

Suppose we get the following joint observations t = (x, y).

\begin{matrix}x & y \\ 1 & 0 \\ 1 & 1 \\ 0 & 0 \\ 1 & 1 \\ 1 & 0 \\ 0 & 1 \\ 0 & 0 \end{matrix}

What are the maximum-likelihood (ML) values of θ_{1} and θ_{2}? (Hint. Since P (x, y|θ) = P (y|x, θ_{2})P (x|θ_{1}), the estimation of the two parameters can be done separately in the log-likelihood criterion.)

Solution:

a. The probability mass function (pmf) of x ∈ {0, 1} is

P (x) = \{ \begin{matrix}θ_{1},  x = 0 \\ 1 − θ_{1},  x = 1 \end{matrix} \}

The conditional pmf of y ∈ {0, 1} given that x = 0 is

P (y|x = 0) = \{ \begin{matrix}θ_{2}, y = 0 \\ 1 − θ_{2}, y = 1 \end{matrix} \}

The conditional pmf of y given that x = 1 is

P (y|x = 1) = \{ \begin{matrix}1 − θ_{2}, y = 0 \\ θ_{2}, y = 1 \end{matrix} \}

Use P (x, y) = P (y|x)P (x) to tabulate the joint pmf of (x, y).

P (x, y) = ( \begin{matrix}P (0, 0) & P (0, 1) \\ P (1, 0) & P (1, 1) \end{matrix} ) = ( \begin{matrix}θ_{2}θ_{1} & (1 − θ_{2})θ_{1} \\ (1 − θ_{2})(1 − θ_{1}) & θ_{2}(1 − θ_{1}) \end{matrix} )

b. We select (θ_{1}, θ_{2}) to maximize the log-likelihood of the samples {(x_{i}, y_{i}), i = 1, ..., n} which may be expressed as

J(θ_{1}, θ_{2}) = ∑_{i} log P (x_{i}, y_{i})
= ∑_{i} (log P (y_{i}|x_{i}) + log P (x_{i}))

= (∑_{i} log P (y_{i}|x_{i})) + (∑_{i} log P (x_{i}))
= J_{2}(θ_{2}) + J_{1}(θ_{1})

Hence, we choose θ_{1} to maximize

J_{1}(θ_{1}) = ∑_{i} log P (x_{i})
= N(x = 1) log(1 − θ_{1}) + (n − N(x = 1)) log θ_{1}

where N(x = 1) = ∑_{i} x_{i}. Differentiating w.r.t. θ_{1} gives

\frac{∂J_{1}}{∂θ_{1}} = \frac{−N(x = 1)}{1 − θ_{1}} + \frac{n − N(x = 1)}{θ_{1}}

We set this derivative to zero and solve for θ_{1} to obtain

θ̂_{1} = 1 − \frac{N(x = 1)}{n}

Similarly, we choose θ_{2} to maximize

J_{2}(θ_{2}) = ∑_{i} log P (y_{i}|x_{i})
= N(x = y) log θ_{2} + (n − N(x = y)) log(1 − θ_{2})

where N(x = y) = ∑_{i} (x_{i}y_{i} + (1 − x_{i})(1 − y_{i})). Differentiating J_{2} w.r.t. θ_{2}, setting to zero and solving for θ_{2} gives

θ̂_{2} = \frac{N(x = y)}{n}

For the example data, θ̂_{1} = \frac{3}{7}, θ̂_{2} = \frac{4}{7}. Thus,

P̂(x, y) = ( \begin{matrix}θ̂_{2}θ̂_{1} & (1 − θ̂_{2})θ̂_{1} \\ (1 − θ̂_{2})(1 − θ̂_{1}) & θ̂_{2}(1 − θ̂_{1}) \end{matrix} )

The maximum likelihood of the data under this model is

∏_{i} P̂(x_{i}, y_{i}) = (\frac{4}{7})^{8} (\frac{3}{7})^{6} ≈ 7.044 × 10^{−5}

Problem 7: 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:

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.

• What choice should you make in order to maximize your chances of receiving the ring of power? Explain your answer.

• 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.

Solution:

• 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 | 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) + 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 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.

• Here is a sample simulation output for the Monty Hall problem:

a c t u a l : 1
guess1 : 2
r e v e a l : 3
swap : 0
guess2 : 2
a c t u a l : 3
guess1 : 3
r e v e a l : 1
swap : 0
guess2 : 3
a c t u a l : 2
guess1 : 3
r e v e a l : 1
swap : 0

guess2 : 3
swap  : 0
win  : 292
l o s e  : 708
win /( win+l o s e ) : 0.292
a c t u a l : 3
guess1 : 1
r e v e a l : 2
swap : 1
guess2 : 3
a c t u a l : 1
guess1 : 1
r e v e a l : 2
swap : 1
guess2 : 3
a c t u a l : 3
guess1 : 2
r e v e a l : 1
swap : 1
guess2 : 3
swap  : 1
win  : 686
l o s e  : 314
win /( win+l o s e ) : 0.686

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

from m a t p l o t l i b . numerix import ∗
from numarray import ∗
from pylab import plot , subplot , legend , axis , xlabel , ylabel , text , show , rand
Error . setMode ( a l l=None , overflow =’warn ’ , underflow =’ ignore ’ , d i v i d e b y z e r o =’warn ’ , i n v a l i d =’warn ’ )
from LinearAlgebra import ∗

f o r swap in range ( 2 ) :
win = 0
l o s e = 0
f o r i in range (1000) :
a c t u a l = i n t ( rand ()∗3)+1;
guess1 = i n t ( rand ()∗3)+1;
i f guess1 = a c t u a l :

r e v e a l = i n t ( rand ()∗2)+1;
i f r e v e a l = a c t u a l :
r e v e a l = r e v e a l + 1 ;
e l s e :
i f guess1 = 1 and a c t u a l = 2 :
r e v e a l = 3 ;
e l i f guess1 = 1 and a c t u a l = 3 :
r e v e a l = 2 ;
e l i f guess1 = 2 and a c t u a l = 1 :
r e v e a l = 3 ;
e l i f guess1 = 2 and a c t u a l = 3 :
r e v e a l = 1 ;
e l i f guess1 = 3 and a c t u a l = 1 :
r e v e a l = 2 ;
e l i f guess1 = 3 and a c t u a l = 2 :
r e v e a l = 1 ;
i f swap = 1 :
i f guess1 = 1 and r e v e a l = 2 :
guess2 = 3 ;
e l i f guess1 = 1 and r e v e a l = 3 :
guess2 = 2 ;
e l i f guess1 = 2 and r e v e a l = 1 :
guess2 = 3 ;
e l i f guess1 = 2 and r e v e a l = 3 :
guess2 = 1 ;
e l i f guess1 = 3 and r e v e a l = 1 :
guess2 = 2 ;
e l i f guess1 = 3 and r e v e a l = 2 :
guess2 = 1 ;
e l s e :
guess2 = guess1 ;
i f guess2 = a c t u a l :
win = win + 1 ;
e l s e :
l o s e = l o s e + 1 ;
# only p r i n t t r a c e f o r f i r s t 3 games
i f i < 3 :
p r i n t ’ a c t u a l : ’ , a c t u a l
p r i n t ’ guess1 : ’ , guess1
p r i n t ’ r e v e a l : ’ , r e v e a l
p r i n t ’ swap : ’ , swap
p r i n t ’ guess2 : ’ , guess2
# p r i n t r e s u l t s f o r each game play p o l i c y

p r i n t ’ swap  : ’ , swap
p r i n t ’ win  : ’ , win
p r i n t ’ l o s e  : ’ , l o s e
p r i n t ’ win /( win+l o s e ) : ’ , f l o a t ( win ) / f l o a t ( win + l o s e )

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

f o r swap = 0 : 1
win = 0 ;
l o s e = 0 ;
f o r i = 1:1000
a c t u a l = f l o o r ( rand ()∗3)+1;
guess1 = f l o o r ( rand ()∗3)+1;
i f guess1 = a c t u a l
r e v e a l = f l o o r ( rand ()∗2)+1;
i f r e v e a l = a c t u a l
r e v e a l = r e v e a l + 1 ;

end e l s e

i f guess1 = 1 && a c t u a l = 2
r e v e a l = 3 ;
e l s e i f guess1 = 1 && a c t u a l = 3
r e v e a l = 2 ;
e l s e i f guess1 = 2 && a c t u a l = 1
r e v e a l = 3 ;
e l s e i f guess1 = 2 && a c t u a l = 3
r e v e a l = 1 ;
e l s e i f guess1 = 3 && a c t u a l = 1
r e v e a l = 2 ;
e l s e i f guess1 = 3 && a c t u a l = 2
r e v e a l = 1 ;

end end

i f swap = 1
i f guess1 = 1 && r e v e a l = 2
guess2 = 3 ;
e l s e i f guess1 = 1 && r e v e a l = 3
guess2 = 2 ;
e l s e i f guess1 = 2 && r e v e a l = 1
guess2 = 3 ;
e l s e i f guess1 = 2 && r e v e a l = 3
guess2 = 1 ;
e l s e i f guess1 = 3 && r e v e a l = 1
guess2 = 2 ;
e l s e i f guess1 = 3 && r e v e a l = 2
guess2 = 1 ;

end e l s e

guess2 = guess1 ;

end

i f guess2 = a c t u a l
win = win + 1 ;

e l s e

l o s e = l o s e + 1 ;

end

%% only p r i n t t r a c e f o r f i r s t 3 games
i f i <= 3

a c t u a l guess1 r e v e a l swap guess2

end end

%% p r i n t r e s u l t s f o r each game play p o l i c y

swap

win / ( win + l o s e )

end Figure 2: The original data and the data transformed into the coordinate system defined by the eigenvectors of their covariance matrix.

Figure 3: The probability density function of the sum of two independent uniform random variables.

Figure 4: 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.

Figure 5: 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. N (0, 1) random variables.