Introduction to Neural Networks

Problem Set 4

Given the three discrete probability distributions for the random variables X ∈ {1,2,3,4,5},

pdist1 = N[{11,3,2,1,10}/Total[{1,10,2,3,1}]];
pdist2 = N[{1,5,2,3,3}/Total[{1,10,2,3,1}]];
pdist3 = N[Reverse[{1,5,2,3,3}]/Total[{1,10,2,3,1}]];

GraphicsRow[{BarChart[pdist1,AxesLabel->{X,p},ChartLabels->{Range[5]}],BarChart[pdist2,AxesLabel->{X,p},ChartLabels->{Range[5]}],BarChart[pdist3,AxesLabel->{X,p},ChartLabels->{Range[5]}]}]

a) Find the most probable X values (modes). You can use: argmax[x_]:=Position[x,Max[x]]

b) Calculate the expectations (mean) values of X for each distribution. Note how the modes differ from the means.

c) Calculate the variances of X. You can use the syntax: ({1,2,3,4,5}-mean)^2

d) Like variance, entropy is a useful measure of uncertainty, or how “random” a distribution is. In some limited cases, entropy is monotonically related to variance. However, because entropy doesn’t depend on the values of X, only on their probabilities, one has to be careful about assuming a relationship.

If all five of the X’s had probability = 1/5, the entropy would be maximum. If  one of the X’s had a probability of 1, and the others 0, the entropy would be minimum.

Calculate the entropies of each of the pairs pdist1,pdist2, and pdist3 (Use Log base e, so the entropy is in the units of nats, rather than bits). Note the differences between the variances and the entropies.

e) Recall that KL divergence was used Boltzmann machine learning to measure how far apart the measured environmental distribution was from the simulated (“dreamed”) distribution. 

Calculate the KL divergences (relative entropies) between all pairs (Use Log base e, i.e. Log[x] has base e as the default).

You should be able to see why the KL measure is called a “divergence” rather than a distance --- true distances are non-negative, and symmetric.

Try to use functional programming methods and avoid loops.

(Subscript[𝒟, 1]) =PDF[NormalDistribution[(Subscript[μ, 1]),σ],x];
(Subscript[𝒟, 2]) =PDF[NormalDistribution[(Subscript[μ, 2]),σ],x];

Hint 1: You’ll need to define and integrate an expression, “myexp”. Rather than trying to directly integrate  Integrate[“myexp”,{x,-Infinity,Infinity}}], help Mathematica out by first using PowerExpand[“myexp”] to simplify the expression for “myexp”.

Hint 2: To avoid having Mathematica return conditions having to do with whether σ is a real number or possibly negative, use: Integrate[“myexp”,{x,-Infinity,Infinity},Assumptions->{σ∈Reals,σ>0}]

Write a function, nsamples[n_],  to draw n samples from a uniform distribution between -1 and 1, and which then sums the n samples to make a new sample.  Write a second function that calls nsamples m=1000 times. Use Manipulate[ ] and HistogramDistribution[ ], and DiscretePlot[ ] to plot up estimates of the probability distribution function of nsamples[n], with n as a slider going from 1 to  20.

Calculate what the standard deviation, σ,  should be as a function of n (using the rule for the sum of variances of i.i.d samples), use this σ, in PDF[NormalDistribution[0,σ],x],  and plot this normal curve on top of the DiscretePlot[].

Note the variance of the UniformDistribution:

Variance[UniformDistribution[{-1,1}]]

And that the variance of n samples (below either 10 or 100) is proportional to n:

Variance[Table[Total[RandomVariate[UniformDistribution[{-1,1}],10]],{10000}]]
Variance[Table[Total[RandomVariate[UniformDistribution[{-1,1}],100]],{10000}]]

Your output should look something like this:

The Appendix in Lecture 14 goes through the example below using list manipulations, which hide the relationship of the computations to the formal definitions of marginalization and conditioning. This exercise uses Mathematica functions to reformulate the problem.

A discrete joint probability

All of our knowledge regarding an inference problem can be described in terms of the joint probability of the hypotheses, H and the possible data measurements, X. The probability function assigns a number to all possible combinations:

p[H, x]

That is, we are assuming that both the hypotheses and the data are discrete random variables.

H={({{S1}, {S2}})

x ϵ{1,2,3}

For example, H could represent two possible signals (H= s1or s2), that are associated with imperfect sensory data  (x = 1,2, or 3). The probability of any pair H and x values occuring together is represented by the joint. Often, such signal discrimination problems are set up assuming priors and likelihoods. Here we’ll start off assuming we know the joint distribution.

Let's assume that X can only take on one of three values, 1, 2, or 3. And suppose the joint probability is:

p[H_,x_]:=Which[H==1 && x ==1,1/12,H==1 && x ==2,1/12,H==1 && x ==3,1/6,H==2 && x ==1,1/3,H==2 && x ==2,1/6, H==2 && x ==3,1/6];

TableForm[Table[p[H,x],{H,1,2},{x,1,3}],TableHeadings->{{"H=s1", "H=s2"},{"X=1","X=2","X=3"}}]

The total probability should sum up to one. Use Sum[ ] to verify

Marginalizing to obtain p(x) and the prior, p(H)

What are the probabilities of the data, p(x)? 

Calculate p(x) by defining a function px[x_] which sums p[H,x] over values of H. 
(If we were using list manipulations we could use the sum rule to sum over the columns of the matrix, as in the Appendix of Lecture 14.) Again, use Sum[ ] to do the  “summing over”, also called marginalization or “integrating out”.

What are the prior probabilities, p(H)? Define a function: pH[H_], and output the values.

Conditioning to obtain the likelihood p[x|H] and posterior p[ H|x]

Use the product rule to obtain the conditional probability through conditioning of the joint:

p[ x|H]=((p[H, x])/(p[H]))

Use the Mathematica function definitions above to do the conditioning. 
Let pxH[H_,x_] represent the probability of x conditioned on H: p[x|H]. 

Show the results in TableForm using: 

TableForm[Table[pxH[H,x],{H,1,2},{x,1,3}], TableHeadings->{{“H=s1”, “H=s2”},{“X=1”,”X=2”,”X=3”}}]

With this table, you can do maximum likelihood estimation. 
What is the most likely value of H given X =1?

But this ignores the difference in the priors. 

Note that the probability of x conditional on H sums up to 1 over x, i.e. each row adds up to 1. But, the columns do not. 
	p[x|H] is a probability function of x, but a likelihood function of H.

Now calculate the posterior probability by conditioning the joint on x. Let pHx[H_,x_] represent p[H|x].  Show the results in TableForm as above.

With this table, you can do MAP estimation.
What is the most probable value of H given X = 1?

p[ H|x]=((p[H, x])/(p[x]))

Multivariate gaussian distributions are often inadequate to model real-life problems, that for example might involve more than one mode. One solution is to approximate more general distributions by a sum or mixture of gaussians.

Clear[mix];

r1=0.4*{{1,.6},{.6,1}};
r2=0.4*{{1,-.6},{-.6,1}};
m1 = {1,.5}; m2 = {-1,-.5};
ndist1 = MultinormalDistribution[m1, r1];
ndist2 = MultinormalDistribution[m2, r2];

mix[x_] := 0.4 PDF[ndist1, x] + .6 PDF[ndist2, x];

gg1=ContourPlot[mix[{x1,x2}],{x1,-2,2},{x2,-2,2},PlotRange->Full,ImageSize->Small]

Marginals for mixture

Write Mathematica code to calculate and then plot the two marginal distributions for the above mixture distribution. Your plot should look like:

By inspecting your plot, estimate the modes for each of the marginal distributions. How do they compare with the mode of the joint distribution?

One method used in exploratory statistics of high dimensional distributions, where one can’t visualize the structure, is to look for “interesting” marginal projections on to lower dimensions, where one can visualize the results. But even without visualization, “interesting” can be defined as departures from gaussianity. In our cased, the x dimension is more interesting than the y direction because the red curve is not at all gaussian, but the green curve looks somewhat gaussian.

A recurring theme in probabilistic computations is to combine information based on knowledge of how the variables are distributed probabilistically. We’ve seen this in the context of combining data with priors, for example to update the prior during learning, and also when passing messages in belief propagation. 

Work in theoretical neuroscience by Ma et al. (2006) showed conditions under which neural populations could achieve optimal estimates by summing spikes ( Bayesian inference with probabilistic population codes. Nature Neuroscience, 9, 1423-8”). Their model was motivated in part by psychophysical observations showing that humans often combine sensory cues optimally.

Suppose there are two cues ((Subscript[x, 1]), (Subscript[x, 2])) to the depth, μ,  of an object. The generative model can be described by: (Subscript[x, i]) = μ + (Subscript[n, i]), where (Subscript[n, i]) ~ N[0, (Subscript[σ, i])], for i = 1, 2. In other words, assume the Gaussian case, with conditional independence (as represented in graph below). 

Assume you’ve separately measured the means,  (Subscript[μ, i]) and standard deviations, (Subscript[σ, i]) of your depth estimates for each cue. Show that the maximum a posteriori estimate of the depth using the combined cues is:

μ = ((Subscript[r, 1])/((Subscript[r, 1]) + (Subscript[r, 2])))(Subscript[μ, 1]) +  ((Subscript[r, 2])/((Subscript[r, 1]) + (Subscript[r, 2])))(Subscript[μ, 2])

where (Subscript[r, i]) = 1/(Subscript[σ, i]^2), and (Subscript[μ, 1]) and (Subscript[μ, 2]) are the estimates of the means obtained separately. 

HINT: See Lecture 16 on computing the MAP estimate of the mean given several measurements. Assume the mean of the prior is zero. Obtain an expression for μ  as a function of two measurements, (Subscript[x, 1]) and (Subscript[x, 2]), then take the expectation of both sides of the expression.