Computational Vision
U. Minn. Psy 5036W

Problem set 2:

Probability review

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

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

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.

Try to use functional programming methods and avoid loops.

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 purpose of these exercises is to better understand marginalization and conditioning.

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

For example, H could represent two possible signals (H= s1or s2), that are associated with imperfect sensory data  (x = 1,2, or 3).

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

x ϵ{1,2,3}

The probability of any pair H and x values occurring together is represented by the joint, which could be represented as a 2x3 matrix of probabilities, as shown below. 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"}}]

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

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

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

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.

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.

p[ H|x]=((p[H, x])/(p[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?

Ideal observer analysis

In this section, you will simulate  the signal-known-exactly ideal (SKE).  This observer has a built-in template that matches the signal that it is looking for. The template would be like an image that you might store in memory, which you then use to compare with an incoming pattern. Alternatively, you can think of the template as the image pattern or feature that a neuron is most responsive to. 

The signal is embedded in “white gaussian noise”.   “white” means the pixels are not correlated with each other--intuitively this means that you can’t reliably predict what any one  pixel’s value is from any of the others.

There are some neurons in the the primary cortex of the visual system called “simple cells”. If there are no internal sources of noise, these cells can be modeled as ideal detectors for the patterns that match their receptive fields. Given the right orientation and size, an image stimulus consisting of an sinusoidal grating “patch” (a “gabor patch”)  can be a good match for a given simple cell.  In actual practice, simple cells, and neurons in general are not noise-free, and are not linear. In particular, the response in terms of spike rate can’t be negative.

In a classic Science paper published in 1981, Burgess and colleagues showed that human observers were almost as good as the ideal pattern discriminator for gabor patches of the right spatial frequency and size (See Burgess, A. E., Wagner, R. F., Jennings, R. J., & Barlow, H. B. (1981). Efficiency of human visual signal discrimination. Science, 214 (4516), 93-94.) They interpreted their results in terms of human decisions using visual neurons that were well-tuned to these patterns.

Before you get started

In the exercises below, you have a signal which is represented as a size x size matrix. But we will do our calculations treating the image signal as a vector. How do you get back and forth between images as matrices vs. vectors? You saw the answer in Problem Set 1: use Flatten[] and Partition[].

The target

The target is given to you  below. It is a low contrast sine-wave grating "patch" oriented at 45 deg and vignetted by a gaussian window. This patch is an example of what is called a 2D gabor function.
The x-y dimensions are size x size. We will define both the target signal and the noise to have zero means so that when we add the noise as contrast noise, the mean doesn't change.

Grating[x_,y_,fx_,fy_] := Cos[2.0 Pi (fx x + fy y)];
GratingPatch[x_,y_,fx_,fy_] := Exp[-(fx x)^2 - (fx y)^2]*Grating[x,y,fx,fy]

We'll plot a cross-section to see what it should look like:

Plot[GratingPatch[x, 0, 1, 1], {x, -2, 2}]

OK, now let's make the 2D target image.

size = 32; amplitude = 0.0625;
fx= 1.0; fy = 1.0;
x1=-2; x2=2;
signal = Table[amplitude*GratingPatch[x,y,fx,fy],
{x,x1,x2,(x2-x1)/(size-1)},{y,x1,x2,(x2-x1)/(size-1)}];

And it looks like this at full contrast:

ArrayPlot[signal,Mesh->False,Frame->False,ColorFunction->"GrayTones"]

...and like this at 6.25% contrast:

ArrayPlot[signal,Mesh->False,Frame->False,ColorFunction->"GrayTones", PlotRange->{-1,1}]

(Note: On an 8-bit display, graylevels go in integer steps from 0 to 255.  We use the option PlotRange->{-1,1} to map our signal values {-1,1} to pixel intensity values L, {0,255}. Contrast of a grating is defined as δL/(Subscript[L, 0]), where  (Subscript[L, 0]) is the mean intensity, and δL is the amplitude of the grating.  Zero maps to graylevel (Subscript[L, 0])= 128 (rounding up). A signal amplitude of 0.0625 maps to 0.0626*128 = 8. So the contrast, δL/(Subscript[L, 0])= 8/128 ~ 0.0626*128/128 = 0.0625. 

If you were actually doing psychophysics, it would be important to check to make sure that the computer monitor can actually supply the graylevel differences you want. E.g.if the amplitude gets below 1/128, an 8-bit display card won't be able to deliver--and human contrast thresholds can be lower than this. Another note, is that the monitor's intensity isn't necessarily linearly related to the requested graylevel. This problem can be fixed by doing "gamma correction", which involves setting the lookup table in your computer's graphics card to compensate for the particular display you may be using.)

Define a function noise that will generate a size x size array of Gaussian noise with standard deviation of sigma= 0.25, and a mean of 0.  (Use the built-in function, NormalDistribution[]). Make sure that you define it in such a way that each time you call noise, it gives you a new random image.

Add the signal to the noise to generate an observation vector called observation, and use ArrayPlot[] (with the options used above) to display a sample observation.

Create an input cell corresponding to the formula:

dprime[s_,sigma_] := Sqrt[s.s]/sigma

to calculate the ideal's d' using Mathematica.

Psychophysicists sometimes design their experiments to look for the signal strength (e.g. contrast) that gives 76% correct (rather than, say 75 which is half way between 50 and 100%). The reason is that 76% is very close to 76.025%, for which dprime2 is exactly 1. Check it out by calculating dprime2[.76025]

Generate ntrials = 100 observation pairs of signal+noise, and  noise. For each pair compute an ideal decision variable (i.e. the output of a cross-correlator, which is the dot product of the observation with the signal). Pick which observation has the higher value of the decision variable. Decide whether this decision is the right answer or not. Total up the proportion of right answers and compute d' using dprime2[]. Repeat this four times (total of 400 trials) to get four estimates of d'. How does the average of these four compare to your result from the closed form prediction above?

fx= 1.0; fy = 1.4;
x1=-2; x2=2;
wrongtemplate = Table[amplitude*GratingPatch[x,y,fx,fy],
{x,x1,x2,(x2-x1)/(size-1)},{y,x1,x2,(x2-x1)/(size-1)}];

ArrayPlot[wrongtemplate,Mesh->False,Frame->False,ColorFunction->"GrayTones"]

You don’t need to write anything for this section. Make sure you understand and can define the following terms.

hit rate, false alarm rate, ROC, statistical efficiency, kurtosis, maximum likelihood estimation, signal-known-exactly ideal, cross-correlation, joint probability, conditional probability, conditioning, marginalization