Introduction to Neural Networks
U. Minn. Psy 5038

Gaussian generative models, learning, and inference

Initialize standard library files:

Off[General::spell1];

Last time

Review of probability and statistics

Basic rules of probability

Basic rules of inference: Given a joint distribution, condition on what you know, and integrate out what you don’t care about.

...but easier said than done.

Rationale: The tools of probabilistic modeling can provide a deeper understanding of neural networks. They also provide quantitative descriptions closer to functional behaviors.

Today

Examples of computations on continuous probabilities
Examples of computations on discrete probabilities
Low dimensional sampling
Introduction to Bayes learning

Mathematica functions for multivariate distributions & exploring marginals

Multivariate gaussian probability density

An n-variate multivariate gaussian (multinormal) distribution with mean vector μ and covariance matrix Σ is denoted (Subscript[N, n])(μ,Σ). The density is:

p(x)=(1/((((2π))^(n/2))((Det[Σ])^(1/2))))Exp[-(1/2)(((x-μ))^T)(Σ^(-1))(x-μ)]

With Mathematica’s built-in function you can, for example, define the bivaraite probability density function with mean vector {(Subscript[μ, 1]),(Subscript[μ, 2])}, and covariance matrix {{((Subscript[σ, 11])^2),ρ*(Subscript[σ, 11])*(Subscript[σ, 22])},{ρ*(Subscript[σ, 11])*(Subscript[σ, 22]),((Subscript[σ, 22])^2)}}, where ρ parameterizes correlation.

Σ = {{((Subscript[σ, 11])^2),ρ*(Subscript[σ, 11])*(Subscript[σ, 22])},{ρ*(Subscript[σ, 11])*(Subscript[σ, 22]),((Subscript[σ, 22])^2)}};
PDF[MultinormalDistribution[{(Subscript[μ, 1]),(Subscript[μ, 2])},Σ ],{x,y}]

Σ

Covariance[MultinormalDistribution[{(Subscript[μ, 1]),(Subscript[μ, 2])},Σ]]

Examples of PDF, CDF

A specific example, MultinormalDistribution[μ,Σ]

m1 = {1,1/2};
r=(1/2)*{{1,2/3},{2/3,4}};
ndist = MultinormalDistribution[m1, r];
pdf = PDF[ndist, {x1, x2}]

GraphicsRow[{g1=ContourPlot[PDF[ndist,{x1,x2}],{x1,-3,3},{x2,-3,3}],g13D = Plot3D[PDF[ndist,{x1,x2}],{x1,-3,3},{x2,-3,3}]}]

Calculating probabilities using the CDF

What is the probability of (Subscript[x, 1]) and  (Subscript[x, 2]) taking on values in the region (Subscript[x, 1])<.5⋂(Subscript[x, 2])<2? 
Recall that the cumulative distribution function is given by:

CDF((Subscript[x, 1]),(Subscript[x, 2])) = (∫_(-∞)^(Subscript[x, 2]))(∫_(-∞)^(Subscript[x, 1]))p((Subscript[x, 1]),(Subscript[x, 2]))d(Subscript[x, 1])d(Subscript[x, 2])

So the answer is the volume under the PDF shown below.

grp=RegionPlot[x1<.5 && x2<2,{x1,-4,4},{x2,-4,4},PlotStyle->Directive[Opacity[.25],EdgeForm[],FaceForm[Gray]]];

Show[{ g1,grp},ImageSize->Small]

We could numerically integrate to find the volume. Alternatively, the answer can be found from the built-in cumulative distribution function, or CDF. The value corresponds to the height of the contour at {(Subscript[x, 1]), (Subscript[x, 2])} = {.5, 2.0}. Move your mouse cursor over the contour plot below near the {.5, 2.0} point.

gcdf = ContourPlot[CDF[ndist, {x1, x2}], {x1, -4, 4}, {x2, -4, 4},ImageSize -> Small];
Show[{ gcdf, grp}, ImageSize -> Small]

A more precise answer is:

CDF[ndist, {.5, 2.0}]

Finding the mode

The mode is the value of the random variable with the highest probability or probability density. For discrete distributions, think of it as the most frequent value.

(Sometimes the word “mode” is used to refer to a local maximum in a density function. Then the distribution is called multimodal. If it has two modes, it is called bimodal.) There may not be a unique mode--the uniform distribution is an extreme case of this. 

For the Gaussian case, the mode vector corresponds to the mean vector. But we can pretend we don't know that, and use the FindMaximum[] function to find the maximum and the coordinates where the max occurs:

FindMaximum[PDF[ndist,{x1,x2}],{{x1,0},{x2,0}}]

Marginals

Suppose we want the distribution for just (Subscript[x, 1]), PDF[(Subscript[x, 1])] = marginal[(Subscript[x, 1])]. This is called calculating the marginal distribution of (Subscript[x, 1]). How do we find it?  It is obtained by integrating out the other variable, (Subscript[x, 2]). I.e.Integrate PDF[(Subscript[x, 1]), (Subscript[x, 2])] with respect to (Subscript[x, 2]). Similarly, we can calculate PDF[ (Subscript[x, 2])].

Clear[x1,x2];
marginal[x1_]:=(∫_(-∞)^∞)PDF[ndist,{x1,x2}]\[DifferentialD]x2;
marginal2[x2_]:=(∫_(-∞)^∞)PDF[ndist,{x1,x2}]\[DifferentialD]x1;

What is the mode of PDF[(Subscript[x, 2])]?

FindMaximum[marginal[x1],{{x1,0}}]
FindMaximum[marginal2[x2],{{x2,0}}]

Now let’s plot up the marginals on top of the contour plot of the joint distribution.

mt=Table[{x1,marginal[x1]},{x1,-3,3,.2}];
g2=ListPlot[mt,Joined->True,PlotStyle->{Red,Thick},Axes->False];

mt2=Table[{x2,marginal2[x2]},{x2,-3,3,.4}];
g3=ListPlot[mt2,Joined->True,PlotStyle->{Green,Thick},Axes->False];

theta=Pi/2;
Show[g1,Epilog->{Inset[g2,{0,-3},{0,0}],Inset[g3,{-3,0},{0,0},Automatic,{{Cos[theta],Sin[theta]},{Sin[theta],-Cos[theta]}}]}]

Drawing samples

As we've used in earlier lectures, drawing samples is done by:

RandomVariate[ndist]

(RandomReal[ndist] also works, but RandomVariate[] is preferred.)

Here’s an example of a “white gaussian noise” image. It is called “white” because there are no correlations between the pixel intensities. Each one is drawn independently of any of the others.

width = 32;
data=RandomVariate[NormalDistribution[0,1],width*width];
Image[Partition[data,width]]//ImageAdjust

Mixtures of gaussians  with MultinormalDistribution[]

Multivariate gaussian distributions are often inadequate to model real-life problems, that for example might involve more than one mode or might have non-gaussian properties. 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_,mp_] := mp*PDF[ndist1, x] + (1-mp)*PDF[ndist2, x];

Manipulate[gg1=Plot3D[mix[{x1,x2},mp],{x1,-2,2},{x2,-2,2},PlotRange->Full,ImageSize->Small],{{mp,.2},0,1}]

Mathematica has a built-in function for defining mixture distributions:

𝒟[w_]=MixtureDistribution[{w,1-w},{MultivariatePoissonDistribution[7,{9,10}],MultivariatePoissonDistribution[2,{5,4}]}];
Manipulate[DiscretePlot3D[PDF[𝒟[w],{x,y}],{x,0,30},{y,0,30},ExtentSize->Full],{w,0,1}]

See Zoran and Weiss (2011) for an interesting application of mixture models to discovering visual features in natural images.

Sampling in low dimensions

Why sample? Rationale for Monte Carlo. 

	Inferring means, marginalization require integrations or sums. 
		But Integrations, especially in high dimensional spaces can be computationally hard. 
	Data synthesis can provide a “reality” check.
	
Later we’ll look at the problem of sampling in high dimensions. 

Let’s look at some methods of sampling. We first go over the simple one dimensional (univariate) case for densities.

Density mapping theorem

Suppose we have a change of variables that maps a discrete set of x's uniquely to y's:  X->Y.

Discrete random variables

No change to probability function. The mapping just corresponds to a change of labels, so the probabilities p(X)=p(Y).

Continuous random variables

In this case, the form of probability density function changes because we require the probability "mass" to be unchanged: 
p(x)dx = p(y)dy

Suppose, y=f(x)

(Subscript[p, Y])(y)δy=(Subscript[p, X])(x)δx

Transformation of variables is useful in making random number generators for probability densities other than the uniform distribution, such as a Gaussian.

Below we’ll need to use the cumulative distribution function: CDF(x) = prob (X<x)=(∫_(-∞)^x)p(X)dX

Univariate sampling

Making a univariate (scalar) gaussian random number generator:

We assume we have a random number generator that provides uniformly distributed numbers between 0 and 1. How can we get numbers that are Gaussian distributed?

Well, the easiest way is to use a built-in function:

RandomVariate[NormalDistribution[0,1]]

but we’d like to better understand some principles behind generating random numbers for a specified distribution.

Method 1: Just for Gaussian. Use Central Limit Theorem

If all we want to do is make a Gaussian random number generator from a uniformly distributed generator, we can use the Central Limit Theorem. The Central Limit Theorem says that the sum of a sufficiently large number of independent random variables drawn from the same underlying distribution (with finite mean and variance), will be approximately normally distributed. Note: These draws are said to be “i.i.d.”, meaning “independent, identically distributed”.

The approximation gets better as the number of samples increases.

Try the cell below with nusamples = 1, 2, .., 10,...

Manipulate[zl=Table[(∑_(i=1)^nusamples)RandomReal[]-(nusamples/2),{1000}];
Histogram[zl,ImageSize->Small],{nusamples,1,30,1}]

Method 2:  Use Density Mapping theorem. More general.

We'll use the density mapping theorem to turn uniformly distributed random numbers RandomReal[] into gaussian distributed random numbers with mean =0 and standard deviation =1.

(Subscript[p, Y])(y)δy=(Subscript[p, X])(x)δx

(Subscript[p, Y])(y)((δy)/(δx))=(Subscript[p, X])(x)

Suppose (Subscript[p, Y])(y)=1 (over the unit interval, but zero elsewhere). Then

y(x)=(∫_(-∞)^x)(Subscript[p, x])(x')dx'=P(x)

Thus if we sample from the uniform distribution to get y, x should be distributed according to (Subscript[p, X])(x). 
To do this, we need a mapping from y → x. This is given by the inverse cumulative distribution, i.e. (P^(-1))(y).
Let's implement this. The quick way is to use Mathematica's built-in function, InverseErf[ ],  to get the inverse cumulative.

Method 2a: Applied to Gaussian

InverseErf[ ] is the inverse of the classic “error function”:

erf(z)=(2/(Sqrt[π]))(∫_0^z)(e^(-(t^2)))dt

We can use this to define a function for the inverse cumulative of a gaussian:

Clear[z];
z[p_]:=(Sqrt[2]) InverseErf[1-2 p];
Plot[z[y],{y,0,1}]

binsize=0.1;
zl=Table[z[RandomReal[]],{5000}];
freq=BinCounts[zl,{-3,3,binsize}];
ListPlot[freq,Filling->Axis]

Method 2b: Cumulative from scratch: Works for almost any low dimensional distribution.

Suppose we have a discrete approximation of any cumulative distribution. How can we generate samples? For illustration purposes, we'll illustrate the method with a discretization of the Gaussian. 

Our first goal is to produce a discrete approximation to the cumulative gaussian. To review where things come from, we'll start with the definition of a Gaussian, and make sure it is normalized.

Integrate[Exp[-(x-x0)^2/(2*σ^2) ],{x,-Infinity,Infinity}]

Let x0=0 and σ=1:

Plot[((E^(-((x1^2)/2)))/(Sqrt[2 π])),{x1,-4,4}]

Note that Plot[PDF[NormalDistribution[0,1],x1],{x1,-4,4}]; gives the same thing using the built-in normal distribution function.

Cumulative gaussian

Clear[cumulgauss,x,x1];
cumulgauss[x_]:=NIntegrate[Exp[-(x1^2 )/2]/(Sqrt[2*Pi]),{x1,-Infinity,x}]

cumulgauss[Infinity]

We can plot up cumulgauss:

Plot[cumulgauss[x],{x,-4,4}]

Now make a discrete version of the cumulative distribution:

lcumulgauss=Table[{x,cumulgauss[x]},{x,-4.`,4.`,0.25`}];
ListPlot[lcumulgauss,Filling->Axis]

Remember the main reason we did the above is to now show how we can go from an arbitrary histogram to drawing samples. I.e. it doesn’t have to be a gaussian at this stage, the cumulative could have come from some any one dimenstional histogram from data gathered elsewhere.

Make inverse cumulative gaussian table

Here is a useful trick whenever you want an inverse function, given a discrete representation.

invlcumulgauss=RotateLeft[lcumulgauss,{0,1}];

To see what this does, evaluate:

{{x1,y1},{x2,y2},{x3,y3}}
RotateLeft[{{x1,y1},{x2,y2},{x3,y3}},{0,1}]

ListPlot[invlcumulgauss,Filling->Axis]

Make a smooth interpolated function from the discrete inverse cumulative

Another useful trick.

Interpolation works by fitting polynomial curves to the data. Try the test below with various interpolation orders (the default is 3)

test=Interpolation[{{1,2.},{2,4},{3,9},{4,16}},InterpolationOrder->2];
Plot[test[x],{x,1,4}]

interinvlcumulgauss=Interpolation[invlcumulgauss];

Plot[interinvlcumulgauss[x],{x,0.01,0.99}]

Draw samples with a standard deviation of Sqrt[10]

Round[10 interinvlcumulgauss[RandomReal[]]]

Draw a bunch of samples, and plot up histogram

z=Table[Round[10 interinvlcumulgauss[RandomReal[]]],{10000}];
domain=Range[-20,20];
Freq=(Count[z,#1]&)/@domain;

ListPlot[Freq,Filling->Axis]

digression...a quick & dirty way to smooth is to do a moving average

ListPlot[MovingAverage[Freq,6],Filling->Axis]

Plot up cumulative histogram

CumFreq=FoldList[Plus,0,Freq];
ListPlot[CumFreq,Filling->Axis]

Same thing, with Accumulate[] :

CumFreq=Accumulate[Freq];
ListPlot[CumFreq,Filling->Axis]

Graphical Models

Causal structure and conditional independence

In general, we’d like to be able to specify natural  patterns such as images by a high-dimensional joint probability. This is usually too difficult in practice, and requires simplification and approximations. If we have some idea of the conditional relationships between various causes and data, this knowledge can be represented in a directed graph.  (Markov Random Fields are represented by undirected graphs.) The idea is to represent the probabilistic structure of the joint distribution P(S,L,I) by  a graphical model that expresses how variables influence each other.  Random variables are represented by nodes, and the links or arrows between the nodes represent conditional relationships.

The above figure illustrates several types of graphical structures. A. undirected graph as in a Markov Random Field. B-E illustrate directed graphs: B. Simple influence relationships. C. Markov process. D. Simple hierarchical model. E. Common vs. independent causes.

Conceptual illustration from vision

Image[]

Figure from: Kersten, D., & Yuille, A. (2003). 150–158.

Consider  three basic building blocks: converging, diverging, and intermediate nodes in the figure below.  Multiple (e.g. scene) variables causing a given image measurement (left), a single variable producing multiple image measurements (middle), or a cause indirectly influencing an image measurement through an intermediate variable (right). These are directed graphs because the arrow shows the direction of influence.

For example, these nodes could correspond to:multiple (scene) causes {shape S1, illumination S2 giving rise to the same image measurement, I ; one cause, S influencing more than one image measurement, {color, I1, brightness, I2}; a scene (or other) cause S, {object identity, S} influencing an image measurement (image contour) through an intermediate variable L (3D shape) .

The arrows tell us how to factor the joint probability into conditionals. So for the three examples above, we have:
	p(S1, S2, I) = p(I | S1, S2) p(S1) p(S2)
	p(S, I1, I2) = p(I1 | S) p(I2 | S) p(S)
	p(S, L,  dI) = p(I | L) p(L|S) p(S)

Primary, secondary variables.

Influences between variables are represented by conditioning, and a graphical model expresses the conditional independencies between variables.  Two random variables may only become independent, however, once the value of some third variable is known.  This is called conditional independence. Recall  that two random variables are independent if and only if their joint probability is equal to the product of their individual probabilities.  Thus, if p(A,B) = p(A)p(B), then A and B are independent.  If p(A,B|C) = p(A|C)p(B|C), then A and B are conditionally independent.  

Here’s an example to help build intuition. When corn prices drop in the summer, hay fever (allergies) incidence goes up.  However, if the joint on corn price and hay fever is conditioned on ``ideal weather for corn and ragweed'', the correlation between corn prices and hay fever drops.  This is because corn price and hay fever symptoms are conditionally independent.

You might have heard the saying “correlation is not causation”.  There is a correlation between eating ice cream and drowning. Why? What event could you condition on to make the dependence go away?

What is noise? Primary and secondary variables.

Noise is whatever variables you don't care to estimate, but contributes to the data. These variables are also called nuisance variables or confounding variables. Yuille and Kersten describe variables as primary and secondary in the context of visual inference.

The task determines what is important and what is not. If the task doesn’t require you to estimate a random variable (secondary), then it gets integrated out of the joint distribution to produce a marginal distribution over the variable that you are interested in (primary). 

One counter-intuitive consequence of this is that the most probable (MAP) estimates are not necessarily consistent across tasks.

This is illustrated in the next section.

Optimal Inference and task dependence: Fruit example

(due to James Coughlan; see Yuille, Coughlan, Kersten & Schrater).

The graph specifies how to decompose the joint probability:
 p[F, C, Is, Ic ] = p[ Ic | C ] p[C | F ] p[Is | F ] p[F ]

Generative model: The prior model on hypotheses, F & C

More apples (F=1) than tomatoes (F=2), and:

ppF[F_]:=If[F==1,9/16,7/16];
TableForm[Table[ppF[F],{F,1,2}],TableHeadings->{{"F=a", "F=t"}}]

The conditional  probability cpCF[C|F]:

cpCF[F_,C_]:=Which[F==1 && C ==1,5/9,F==1 && C ==2,4/9,F==2 && C ==1,6/7,F==2 && C ==2,1/7];
TableForm[Table[cpCF[F,C],{C,1,2},{F,1,2}],TableHeadings->{{"C=r","C=g"},{"F=a", "F=t"}}]

The above conditional is a probability distribution on C. So would you expect the sum over a row to be 1? Over a column?

So by the product rule the joint is:

jpFC[F_,C_]:=cpCF[F,C] ppF[F];
TableForm[Table[jpFC[F,C],{F,1,2},{C,1,2}],TableHeadings->{{"F=a", "F=t"},{"C=r","C=g"}}]

Note that now all the entries sum to 1.

We can marginalize to get the prior probability on color alone:

ppC[C_]:=(∑_(F=1)^2)jpFC[F,C]

Which color is a priori more probable?

Is fruit identity independent of material color--i.e. is F independent of C?

Check whether the joint probability on Fruit and Color can be factored into the product of the prior probabilities on Fruit and Color.

Answer

No.

TableForm[Table[jpFC[F,C],{F,1,2},{C,1,2}],TableHeadings->{{"F=a", "F=t"},{"C=r","C=g"}}]
TableForm[Table[ppF[F] ppC[C],{F,1,2},{C,1,2}],TableHeadings->{{"F=a", "F=t"},{"C=r","C=g"}}]

Generative model: The likelihood model--probabilities of measurements, i.e. some features given hypotheses

Suppose that we have gathered some "image statistics" which provides us knowledge of how the image measurements for shape (Is), and for color (Ic) depend on the type of fruit F, and material color, C. For simplicity, our measurements are discrete and binary , say Is = {am, tm}, and Ic = {rm, gm}. (In a more realistic case, they would have continuous values).

P(I_S= am, tm | F=a) = {11/16, 5/16}
P(I_S= am, tm | F=t) = {5/8, 3/8}
P(I_C= rm, gm | C=r) = {9/16, 7/16}
P(I_C= rm, gm | C=g) = {1/2, 1/2}
We use the notation am, tm, rm, gm because the measurements are already suggestive of the likely cause. So there is a correlation between apple and apple-like shapes, am; and between red material, and "red" measurements, rm. 
(This may sound too artificial, but a red apple can result in a greenish measurement if illuminated with a greenish light.)
In general, there may not be an obvious correlation like this.

We define a function for the  probability of Ic given C,  cpIcC[Ic | C]:

cpIcC[Ic_,C_]:=Which[Ic==1 && C ==1,9/16,Ic==1 && C ==2,7/16,Ic==2 && C ==1,1/2,Ic==2 && C ==2,1/2];
TableForm[Table[cpIcC[Ic,C],{C,1,2},{Ic,1,2}],TableHeadings->{{"Ic=rm","Ic=gm"},{"C=r", "C=g"}}]

The  probability of Is conditional on F is cpIsF[Is | F]:

cpIsF[Is_,F_]:=Which[Is==1 && F ==1,11/16,Is==1 && F ==2,5/8,Is==2 && F ==1,5/16,Is==2 && F ==2,3/8];
TableForm[Table[cpIsF[Is,F],{Is,1,2},{F,1,2}],TableHeadings->{{"Is=am","Is=tm"},{"F=a", "F=t"}}]

The total joint probability

We now have enough information to put probabilities on the 2x2x2 "universe" of possibilities, i.e. all possible combinations of fruit, color, and image measurements. Looking at the graphical model makes it easy to use the product rule to construct the total joint, which is:
 p[F, C, Is, Ic ] = p[ Ic | C ] p[C | F ] p[Is | F ] p[F ]:

jpFCIsIc[F_, C_, Is_, Ic_ ] := cpIcC[ Ic, C ] cpCF[F,C ] cpIsF[Is,F ] ppF[F ]

Usually, we don't need the probabilities of the image measurements (because once the measurements are made, they are fixed and we want to compare the probabilities of the hypotheses. But in our simple case here, once we have the joint, we can calculate the probabilities of the image measurements through marginalization p(Is,Ic)=(∑_C)(∑_F)p(F,C,Is,Ic), too:

jpIsIc[Is_,Ic_]:=(∑_(C=1)^2)(∑_(F=1)^2)jpFCIsIc[F, C, Is, Ic ]

Three MAP tasks

We are going to show that the best guess depends on the task. 
In other words, given measurements Is, Ic, the most probable choices of fruit and/or color depend on what which combination of hypotheses we care about.

Define argmax[] function:

argmax[x_]:=Position[x,Max[x]];

This returns the position of the biggest value in a list.

TASK 1: Pick most probable fruit AND color--Answer "red tomato"

We are given some data--i.e. values of Is and Ic and want to draw some conclusions about what kind of fruit we are looking at, and its material color. First, suppose the task is to make the best bet as to the fruit AND material color from the measurements given. To make it concrete, suppose that we see an "apple-ish shape" with a reddish color, i.e., we measure Is=am=1, and Ic = rm=1. The measurements suggest "red apple", but to find the most probable, we need to take into account the priors too in order to make the best guess. 
p(F,C | Is, Ic) is given by:

FCcIsIc[F_,C_,Is_,Ic_]:=jpFCIsIc[F, C, Is, Ic ]/jpIsIc[Is,Ic]

TableForm[FCcIsIcTable=Table[FCcIsIc[F, C, 1, 1 ],{F,1,2},{C,1,2}],TableHeadings->{{"F=a", "F=t"},{"C=r","C=g"}}]
Max[FCcIsIcTable]
argmax[FCcIsIcTable]

So looking at the table we can see that our best bet is “red tomato”.

TASK 2: Pick most probable color--Answer "red"

Same measurements as before. But now suppose we only care about the true material color, and not the identity of the object. Then we want to integrate out or marginalize with respect to the shape or fruit-type variable, F. In this case, we want to maximize the posterior:

p(C | Is=1, Ic=1)=(∑_(F=1)^2)p(F,C|Is=1,Ic=1)

pC[C_,Is_,Ic_]:=(∑_(F=1)^2)FCcIsIc[F, C, Is, Ic ]

pCTable=Table[pC[C,1,1],{C,1,2}];
TableForm[pCTable,TableHeadings->{{"C=r","C=g"}}]
Max[pCTable]
argmax[pCTable]

Answer is that the most probable material color is C = r, "red".

TASK 3: Pick most probable fruit--Answer "apple"

Same measurements as before, Is=am=1, and Ic = rm=1. But now, we don't care about the material color, just the identity of the fruit. Then we want to integrate out or marginalize with respect to the material variable, C. In this case, we want to maximize the posterior:

p(F | Is, Ic)

pF[F_,Is_,Ic_]:=(∑_(C=1)^2)FCcIsIc[F, C, Is, Ic ]

pFTable=Table[pF[F, 1, 1 ],{F,1,2}];
TableForm[pFTable,TableHeadings->{{"F=a","F=t"}}]
Max[pFTable]
argmax[pFTable]

The answer is "apple". So to sum up, for the same data measurements, the most probable fruit AND color is "red tomato", but the most probable fruit is "apple"!

Important "take-home message": Optimal inference depends on the precise definition of the task

Further, the above example shows that that most probable answers from the marginals does not have to be consistent with the most probable answers from the joint distribution.

Try expressing the task-dependent consequences using the frequency interpretation of probability.

How to sample from this model?

Appendix

Exercises

Question: Is p(F, C, Is=1, Ic=1) a joint probability on F and C? (Answer: No.)

Redefine jpFCcIsIcTable and jpFCIsIcTable, without the TableForm[] wrapper and calculate:

Total[jpFCcIsIcTable,{2}]
Total[jpFCIsIcTable,{2}]

Note that we don’t always have to calculate the conditional. For example in Task 1 above.

We can use the fact that p(F,C | Is, Ic) = ((p(F,C,Is,Ic))/(p(Is,Ic)))∝ p(F, C , Is=1, Ic=1). I.e. the conditional is proportional to the function specified by the joint evaluated at Is=1, and Ic=1.

jpFCcIsIcTable=Table[jpFCIsIc[F, C, 1, 1 ],{F,1,2},{C,1,2}];
TableForm[jpFCIsIcTable,TableHeadings->{{"F=a", "F=t"},{"C=r","C=g"}}]
Max[jpFCIsIcTable]
argmax[jpFCIsIcTable]

Using Mathematica lists to manipulate discrete priors, likelihoods, and posteriors

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 or generativemodel 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}}

Marginalization and conditioning: A small dimensional example using list manipulation in 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 use 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]]]

Marginalization and conditioning: An example using Mathematica functions

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[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. 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]].

Sum[p[H,x],{H,1,2},{x,1,3}]

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

Table[p[1,x],{x,1,3}]

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[x_]:=Sum[p[H,x],{H,1,2}];

Table[px[x],{x,1,3}]

"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[H_]:=Sum[p[H,x],{x,1,3}];

Table[pH[H],{H,1,2}]

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

We use function definition in Mathematica to do the division for conditioning. The syntax is simple:

pxH[H_,x_]:=p[H,x]/pH[H];

Table[pxH[H,x],{H,1,2},{x,1,3}]

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

pHx[H_,x_]:=p[H,x]/px[x];

Table[pHx [H,x],{H,1,2},{x,1,3}]

Plotting the joint

We use ArrayPlot[].

ArrayPlot[Table[p[H,x],{H,1,2},{x,1,3}]]

References

Applebaum, D. (1996). Probability and Information . Cambridge, UK: Cambridge University Press.
Cover, T. M., & Joy, A. T. (1991). Elements of Information Theory. New York: John Wiley & Sons, Inc.
Duda, R. O., & Hart, P. E. (1973). Pattern classification and scene  analysis . New York.: John Wiley & Sons.
Golden, R. (1988). A unified framework for connectionist systems. Biological Cybernetics, 59, 109-120.
Intrator, N.  Combining Exploratory Projection Pursuit and Projection Pursuit Regression.  Neural Computation (5):443-455, 1993. http://www.physics.brown.edu/users/faculty/intrator/papers/epp-ppr.ps.gz
Kersten, D. and P.W. Schrater (2000), Pattern Inference Theory: A Probabilistic Approach to Vision, in Perception and the Physical World, R. Mausfeld and D. Heyer, Editors. , John Wiley & Sons, Ltd.: Chichester. (pdf)
Kersten, D., Mamassian P & Yuille A (2004) Object perception as Bayesian inference. Annual Review of Psychology. (pdf, http://arjournals.annualreviews.org/doi/pdf/10.1146/annurev.psych.55.090902.142005)
Kersten, D., & Yuille, A. (2003). Bayesian models of object perception. Current Opinion in Neurobiology, 13(2) http://gandalf.psych.umn.edu/~kersten/kersten-lab/papers/KerstenYuilleApr2003.pdf
Ma, W. J., Beck, J. M., Latham, P. E., & Pouget, A. (2006). Bayesian inference with probabilistic population codes. Nat Neurosci, 9(11), 1432-1438.
Murphy, K. P. (2012). Machine Learning: a Probabilistic Perspective. MIT Press.
Ripley, B. D. (1996). Pattern Recognition and Neural Networks. Cambridge, UK: Cambridge University Press.
Yuille, A., Coughlan J., Kersten D.(1998) (pdf)
Zhu, S. C., Wu, Y., & Mumford, D. (1998). Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling. International Journal of Computer Vision, 27(2), 107–126.
Zoran, D., & Weiss, Y. (2011). From learning models of natural image patches to whole image restoration, 479–486.
http://www.cs.huji.ac.il/~daniez/EPLLICCVCameraReady.pdf