Mathematica Supplement on Probability
U. Minn. 
Daniel Kersten

<<Graphics`Graphics`

Goals

Review the basics of probability distributions, cumulative distributions.
Introduce statistical inference in the context of discriminating two signals (bright or dim switch setting), given a measurement x (e.g. photon count).

Probability

Random variables, discrete probabilities, probability densities, cumulative distributions

Discrete: random variable X can take on a finite set of discrete values

X = {x(1),...,x(N)]

\(TraditionalForm
    \`∑\+\(i = 1\)\%N p\_i = \(∑\+\(i = 1\)\%N p(X = x(i)) = 1\)\)

Densities: X takes on continuous values, x, in some range.

Density :  p(x)\(TraditionalForm
    \`prob(x < X < dx + x) = 
      ∫\_X\%\(dX + X\)\(p(x)\) dx\)prob(x < X < dx + x) ≃ p(x) dx\(TraditionalForm
    \`∫\_\(-∞\)\%∞\( p(x)\) dx = 
      1\)

Cumulative  distribution:

\(TraditionalForm
    \`prob\ \((X < x)\) = 
      ∫\_\(-∞\)\%x\( p(X)\) dX\)

Densities of discrete random variables

The Dirac Delta function, δ[•], allows us to use the mathematics of continuous distributions for discrete ones, by defining the density as:

p[x]=(∑_(i = 1)^N) (Subscript[p, i])δ[x - x[i]], where δ[x - x[i]] ={({{∞}, {0}})({{ for  x =  x[i]}, {  for  x  ≠ x[i]}})

Think of the delta function, δ[•], as ϵ wide and 1/ϵ tall, and then let ϵ -> 0, so that:

\(TraditionalForm
    \`∫\_\(-∞\)\%∞\( δ(y)\) 
          dy = 1\)

The density, p[x], is a series of spikes. It is infinitely high only at those points for which x = x[i], and zero elsewhere. But "infinity" is scaled so that the local mass or area around each point x[i], is (Subscript[p, i]).

Joint probabilities

Prob (X AND Y) = p(X, Y)Joint density :  p(x, y)

Three basic rules of probability

Rule 1: Conditional probabilities from joints: The product rule

Probability about an event changes when new information is gained.
Prob(X given Y) = p(X|Y)

p(X | Y) = ((p(X, Y))/(p(Y)))

p(X, Y) = p(X | Y) p(Y)

The form of the product rule is the same for densities as for probabilities.

Rule 2: Lower dimensional probabilities from joints: The sum rule

p(X) = (∑_(i = 1)^N) p(X, Y(i))

\(TraditionalForm
    \`p(x) = ∫\_\(-∞\)\%∞\( p(x, y)\) 
          dx\)

Rule 3: Bayes' rule

From the product rule, and since p[X,Y] = p[Y,X], we have:

\(TraditionalForm
    \`p(Y | X) = 
      \(\(\(p(X | Y)\)\ \(p(Y)\)\)\/\(p(X)\) = 
        \(\(p(X | Y)\)\ \(p(Y)\)\)\/\(∑\+Y p(X, Y)\)\)\)

Bayes Terminology

p[S | I] = ((p[I | S] p[S])/p)[I]

Usually, we will be thinking of the Y term as a random variable over the hypothesis space, and X as data. So for visual inference, Y = S (the scene), and X = I (the image data), and I = f(S).
We'd like to have:
p(S|I) is the posterior probability of the scene given the image
-- i.e. what you get when you condition the joint by the image data. The posterior is often what we'd like to base our decisions on, because as we discuss below,  picking the hypothesis S which maximizes the posterior (i.e. maximum a posteriori or MAP estimation) minimizes the average probability of error.
p(S) is the prior probability of the scene.
p(I|S) is the likelihood of the scene. Note this is a probability of I, but not of S.

Independence

Knowledge of one event doesn't change the probability of another event. 
p(X)=p(X|Y)
p(X,Y)=p(X)p(Y)

Statistics

Expectation & variance

Definition of expectation or average:

\(Average[X]\  = \ 
      \(X\&_ = \(E[X] = Σ\ x[i]\ p[x[i]]\)\)\)

μ  =  E[X] = ∫x p(x) dx

Some rules:

E[X+Y]=E[X]+E[Y]
E[aX]=aE[X]
E[X+a]=a+E[X]

Definition of variance:

(σ^2)  =  Var[X] = E[[X-μ]^2]\(TraditionalForm
      \`\( = ∑\+\(j = 1\)\%N\((\((p(x(j)))\) \((x(j) - μ)\))\)\^2
            \)\)\(TraditionalForm
      \`\(\( = ∑\+\(j = 1\)\%N\( p\_j\)(x\_j - μ)\))\)\^2\)

Var[X] = ∫(((x - μ))^2) p(x) dx

Standard deviation:

σ = (Sqrt[Var])[X]

Some rules:

Var[X] = E[(X^2)] - E[X]^2Var[aX] = (a^2) Var[X]

Covariance & Correlation

Covariance:

Cov[X,Y] =E[[X - (Subscript[μ, X])] [Y - (Subscript[μ, Y])] ]

Correlation:

ρ[X, Y] = Cov[X, Y]/((Subscript[σ, X]) (Subscript[σ, Y]))

Independent random variables

E[X Y] = E[X] E[Y]Cov[X, Y] = ρ[X, Y] = 0Var[X + Y] = Var[X] + Var[Y]

Degree of  belief vs., relative frequency

What is the probability that the Vikings will win the Superbowl in 2001? Assigning a number between 0 and 1 is assigning a degree of belief. These probabilities are also called subjective probabilities. "Odds" determine subjective probabilities, where the "odds of x to y" means probability = x/(x+y).
What is the probability that a coin will come up heads? In this case, we can do an experiment. Flip the coin n times, and count the number of heads, say h[n], and the set the probability, p = h[n]/n -- the relative frequency . Of course, if we did it again, we may not get the same estimate of p. One solution often given is:

p = (lim_(n -> ∞))((h(n))/n)A problem with this, 
    is that there is no guarantee that a well - 
      defined limit exists.

In visual inference, we can at least in principle, measure statistics, and model probabilities of both image and scene variables. So the relative frequency interpretation seems reasonable. In practice, the dimensions of vision problems are so high, that it is impractical to do this. Once we use the statistical framework to model perception, say of a particular shape from shading information, then probabilities are more like  "subjective unconscious beliefs".

Principle of insufficient reason

Principle of symmetry

Suppose we have N events, x[1],x[2],x[3],...,x[N] that are all physically identical except for the label. The we assume that

\(TraditionalForm
    \`prob(x(1)) = \(prob(x(2)) = \(prob(x(3)) = \(prob(x(N)) = 1\/N\)\)\)\)

In other words,if we have no additional information about the events,we should assume that they are uniformly distributed.

Maximum entropy

Information theory provides a powerful extension to the principle of symmetry. Information of event X is:

Information[X] = -(Subscript[log, 2])(p(X))

Using the definition of expectation above, we can specify the expectation of information, which is called entropy. Entropy of a random variable X with probability distribution p[X] is:

\(TraditionalForm
    \`H(X) = \(Average(Information[X]) = 
        \(-\(∑\+X\( p(X)\)\ \(\(log\_2\)(p(X))\)\)\)\)\)

It can be shown that out of all possible probability distributions, p(X) is biggest for the uniform distribution, p(X)=1/N. It turns out that a more powerful formulation of the principle of symmetry is maximum entropy. For example, out of all possible probability distributions of a random variable with infinity range, but with a specific mean and standard deviation, the Gaussian is unique in having the largest entropy. If the range goes from zero to infinity, and we know the mean, the maximum entropy distribution is an exponential (Cover and Thomas).
An interesting application of the maximum entropy principle is to learning image textures joint probabilities: p(I[1],...,I[N]), where N is very big, but where one has only a relatively small number of measured statistics relative to the number of possible images (which is really huge). The measurements underdetermine the dimesionality of the probability space--i.e. there are many different probability distributions which give the same statistics. So the principle of symmetry, or insufficient reason, says to choose the one with the maximum entropy.

Marginalization and conditioning: A small dimensional example using Mathematica

A discrete joint probability

All of our knowledge regarding the signal discrimination 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, ...}

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

p = {{(1/12), (1/12), (1/6)}, {(1/3), (1/6), (1/6)}}

TableForm[p, 
      TableHeadings
         -> {{"H=S1",  "H=S2"}, {"x=1", "x=2", "x=3"}}]

The total probability should sum up to one. Let's test to make sure. We first turn the list of lists into a singel list of scalars using Flatten[]. And then we can sum either with Apply[Plus,Flatten[p]].

Plus @@ Flatten[p]

We can pull out the first row of p like this:

p[[1]]

{(1/12), (1/12), (1/6)}

Is this the probability of x? No. For a start, the numbers don't sum to one. But we can get it through the two processes of marginalization and conditioning.

Marginalizing

What are the probabilities of the data, p(x)? To find out, we use the sum rule to sum over the columns:

px  =  Apply[Plus, p]

"Summing over "is also called marginalization or "integrating out".  Note that marginalization turns a probability function with higher degrees of freedom into one of lower degrees of freedom.

What are the prior probabilities? p(H)? To find out, we sum over the rows:

pH  =  Apply[Plus,  Transpose[p]]

Conditioning

Now that we have the marginals, we can get use the product rule to obtain the conditional probability through conditioning of the joint:

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

In the Exercises, you can see how to using Mathematica to do the division for conditioning. The syntax is simple:

pxH  =  p/pH

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. The posterior probability is obtained by conditioning on x:

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

Syntax here is a bit more complicated, because the number of columns of px don't match the number of rows of p. We use Transpose[] to exchange the columns and rows of p before dividing, and then use Transpose again to get back the 2x3 form:

pHx  =  Transpose[Transpose[p]/px]

Plotting the joint

The following BarChart[] graphics function requires in add-in package (<< Graphics`Graphics`), which is specified at the top of the notebook. You could also use ListDensityPlot[].

BarChart[p[[1]], p[[2]]];

Basic Inference: SDT Poisson signal discrimination example

Usually for the problems of visual inference, it is hard to directly specify the posterior probability, and it is easier to model the prior probabilities and the likelihoods.  But optimal performance is determined by the posterior probability, and that is where Bayes' rule comes to the rescue. But first, let's consider the roles the priors and likelihoods play separately.

Prior only

Consider the light discrimination problem of lecture 3. Suppose the light transmission switch is set to S1 1/3 of the time, and to S2 2/3 of the time. But the photon counter isn't working at all--so there are no data. What is the observers's optimal strategy? First we have to specify the receiver's goal, and then it makes sense to answer this question. If the receiver wants to minimize the average probability of error, then it should always pick S2.

In maximum a posteriori estimation, the rule is the minimum error rule:
	choose S2 if p(S2) > p(S1)
	choose S1 if p(S1) > p(S2)
So if p(S2) > p(S1), the rule says to always choose S2. Interesting enough, people usually don't do this. With feedback, they tend to match probabilities.
Following the above rule, S1 is presented p(S1) fraction of the time--so this is exactly the error rate, i.e. the fraction of the time that the observer gets the wrong answer using the minimum error rule.
Exercise: What is the probabiliyt of error if the probabilities are matched?

Likelihood

We use Mathematica to define p(x|H) = poisson[x,a], where the hypothesis or switch setting H, determines the average number of samples, x:

poisson[x_,a_] := Exp[-a] a^x / Factorial[x]

Then we make two lists corresponding to means of 4 and 9:

\(p1\  = \ Table[N[poisson[x, 4]], \ {x, 0, 20}]; 

    p2 = \ Table[N[poisson[x, 9]], \ {x, 0, 20}]; \)

And use a function in the Graphics package to plot the distributions.

BarChart[p1, p2, PlotRange -> {0, .25}];

Let's put the two likelihoods together into one 2x21 array:

poispxH  =  {p1, p2};

As we noted in the lecture 3 notes, choosing the hypothesis with the largest likelihood will minimize the probabilty of error when the prior probabilities are uniform--i.e. maximum likeihood decisions are the same as maximum a posterior decisions in this case. But what if the priors are not uniform (i.e. the same constant probability for all outcomes)?

Combining the prior and likelihood: Maximum a posteriori (MAP) rule

Suppose (as we assumed for the light discrimination example) that we know the likelihoods: p(x|S1) and p(x|S2). Bayes' rule tells us how to combine them to obtain the a posteriori probabilbity of the hypotheses conditional on the data:

\(TraditionalForm
    \`p(H | x) = 
      \(\(\(p(x | H)\) \(p(H)\)\)\/\(p(x)\) = 
        \(\(p(x | H)\) \(p(H)\)\)\/\(∑\+H\( p(x | H)\) \(p(H)\)\)\)\)

posterior  =  
      Transpose[Transpose[(poispxH pH)]/Plus@@(poispxH pH)];

Below is a plot of the posterior probabilities, p(S2|x) and  p(S1|x),  as function of the data x:

BarChart[posterior[[1]], posterior[[2]], 
      PlotRange -> {0, 2}];

In maximum a posteriori estimation, the rule is:
	choose S2 if p(S2|x) > p(S1|x)
	choose S1 if p(S1|x) > p(S2|x)

More generally, a function which chooses the value of argument (H = S1 or S2) that maximizes a function is sometimes written as:
 (Underscript[argmax, H]) p(H|x)

MAP minimizes the average probability of error

If we just  studied the the results of decisions using the minimum error rule for a fixed measurement x, we'd find the probability of error to be, say p(error|x). We can get the total probability of error by integrating out x from the joint probability. Using the sum and product rules:

\(TraditionalForm
    \`p(error) = 
      \(∫\_\(-∞\)\%∞\( p(error, x)\) 
            dx = 
        ∫\_\(-∞\)\%∞\( p(error | x)\) \(p(x)\) 
            dx\)\)

So if we always choose the hypothesis that maximizes the posterior p(H|x), p(error|x) would be as small as possible, and thus the integral would be as small as possible (p(x)) isn't affected by the decision strategy we use).

Exercises

A note on list arithmetic

We haven't done standard matrix/vector operations above to do conditioning. We've take advantage of how  Mathematica divides a 2x3 array by a 2-element vector:

M=Array[m,{2,3}]
X = Array[x,{2}]

M/X

Putting the probabilities back together again to get the joint

Transpose[Transpose[pHx] px]

pxH pH

Getting the posterior from the priors and likelihoods:

One reason Bayes' theorem is so useful is that it is often easier to formulate the likelihoods (e.g. from a causal model of how the data could have occurred), and the priors (often from heuristics, or in computational vision empirically testable models of the external visual world). So let's use Mathematica to derive p(H|x) from p(x|H) and p(H) , (i.e. pHx from pxH and pH ).

px2  =  Plus@@(pxH pH)

Transpose[Transpose[(pxH pH)]/Plus@@(pxH pH)]

Show that this joint probability has a uniform prior (i.e. both priors equal).

p = {{1/8, 1/8, 1/4}, {1/4, 1/8, 1/8}}

References

Applebaum, D. (1996). Probability and Information . Cambridge, UK: Cambridge University Press.
Duda, R. O., & Hart, P. E. (1973). Pattern classification and scene  analysis . New York.: John Wiley & Sons.

Notes on getting started with Mathematica

Numerical Calculations. You can do arithmetic. For example, type 5+7 as shown in the cell below, and then hit the "enter" key. Try other operations, 5^3, 4*3 (note that 4 3,
where a space separates the digits is also interpreted as multiplication). Note that if you
try division, e.g. 2/3, you get the exact answer back. To get a decimal approximation,
type N[2/3].

5+7

N[2/3]

You can go back and select an expression by clicking on the brackets on the far right. These brackets are features of the Macintosh interface and serve to organize text and calculations into a Notebook with outlining features. You can group or ungroup cells for text, graphs, and expressions in various ways to present your calculations. Explore these options under Cell in the menu.  You can see the possible cell types under the Style menu.

Built-in functions. Mathematica has a very large library of built-in functions. They all begin with an uppercase letter. You can get information about a function, e.g. for the exponential of a function, or for plotting graphs:

?Exp

?Plot

If you type two question marks before a function, ??Plot, you'll get more information. Try it. What does the Random function do?