Introduction to Neural Networks
Daniel Kersten
Belief Propagation

Initialize

Off[General::spell1];
Needs["ErrorBarPlots`"]

Last time

Bayesian inference and conjugate priors

Graphical models

(Bias/variance trade-off and over-fitting)

Today

Bayesian learning of the mean of a gaussian

Belief propagation on a graph

Final projects

Bayesian learning of univariate Gaussian mean: MAP

We’ll spend more time on unsupervised learning later, but it is useful to see how to apply what we’ve learned about to one of the simplest applications of Bayes to unsupervised learning.

From a statistical point of view, one form of unsupervised learning is "density estimation" which can be done (non-parametrically) from histogram measurements. In high dimensions this is hard, but is easier if we have a low-dimensional parametric model for the density--i.e. the density is modeled in terms of a few parameters. So for example, the 1D Gaussian could be approximated by a  huge list of numbers ("statistics")--a frequency of occurrence for each bin with a range Δx, each number representing an estimate of the probability of the value of the random variable falling in that bin. But because it is Gaussian, we can be more efficient by representing the density in terms of just two numbers (also "statistics", but just the mean and variance), and a formula.  In this context, learning becomes parameter estimation.

A Bayesian learning example: Suppose we believe the data comes from a Gaussian generative process, but we don't know the mean.

This example follows the same logic as the coin example in the previous lecture, except now we work with continuous random variables, and gaussian rather than the binomial distribution.

For simplicity, suppose we have a set of samples that come from a Gaussian distribution with known variance (σ^2), but unknown mean μ.

(Subscript[x, i])=noise, where noise ~N[μ,σ], or equivalently
			(Subscript[x, i])=μ+noise, where noise ~N[0,σ]

ndist0=NormalDistribution[μ,σ];

This is the forward model. We now want to do “inverse probability”, i.e. estimate the mean and its distribution based on prior assumptions and data.

So lets treat the unknown mean,  μ as a random variable, and assume a Gaussian prior on it:

p(μ), with μ~N[(Subscript[μ, 0]),(Subscript[σ, 0])]

ndistμ=NormalDistribution[(Subscript[μ, 0]),(Subscript[σ, 0])];
PDF[ndistμ,μ]

We can think of  μ0 as an initial guess of the mean's mean, with standard deviation (σ0). But we are willing to change our estimate of the mean given new data. If we are really uncertain at the beginning,, we can start of assuming a large standard deviation (Subscript[σ, 0]), and as we gather data, the uncertainty about the value of the mean will decrease.

Suppose the generative model N[μ,σ] produces three i.i.d. (independent, identically distributed) samples (Subscript[x, 1]), (Subscript[x, 2]),  (Subscript[x, 3]).What is the MAP estimate of μ? Which value of  μ makes the posterior biggest? We use Bayes rule:

p(μ|(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]))=((p((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])|μ)p(μ))/(p((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]))))

p((Subscript[x, 1])|μ)is given by:

PDF[ndist0,(Subscript[x, 1])]

Because the samples are drawn independently, the p((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])|μ) is the product of three terms, so the numerator is p((Subscript[x, 1])|μ)p((Subscript[x, 2])|μ)p((Subscript[x, 3])|μ)times the prior p(μ):

PDF[ndist0,(Subscript[x, 1])]*PDF[ndist0,(Subscript[x, 2])]*PDF[ndist0,(Subscript[x, 3])]*PDF[ndistμ,μ]

Calculating the MAP estimate of mean

To find the value of the mean that is largest given our three samples, and our prior assumption, we need to find μ where p((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])|μ)p(μ)is biggest. This is the same as finding the value where the log of  p((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])|μ)p(μ) is biggest--a simpler expression. And with calculus, we can do that by finding where the slope of the derivative of the log with respect to μ  is zero.

g=PDF[ndist0,(Subscript[x, 1])]*PDF[ndist0,(Subscript[x, 2])]*PDF[ndist0,(Subscript[x, 3])]*PDF[ndistμ,μ];
t=Log[g]

t1=PowerExpand[t]

t2=D[t1,μ]

...Mathematica is having problems with subscripts, so get rid of them:

Solve[-((-μ+(Subscript[x, 1]))/(σ^2))+((-μ+(Subscript[x, 2]))/(σ^2))+((-μ+(Subscript[x, 3]))/(σ^2))-(((μ-μ0))/((σ0)^2))==0,μ]

Divide the above expression by ((σ0)^2)(σ^2).

{{μ->((((μ0)/((σ0)^2))+(1/(σ^2))((Subscript[x, 1])+ (Subscript[x, 2]) + (Subscript[x, 3])))/((3/(σ^2))+(1/((σ0)^2))))}}

We've done the calculation with just three data points, but one can see the pattern. So in general, μ can be estimated from n samples in batch mode by the following rule:

{{μ->((((μ0)/((σ0)^2))+(1/(σ^2))(∑_(i=1)^n)(Subscript[x, i]))/((n/(σ^2))+(1/((σ0)^2))))}}

With a bit more effort, one can also derive the solution to estimate the standard deviation too, when it is unknown.

So to sum up, we have a method to estimate the probability distribution of the mean given the data.  (For the multi-variate case, cf. Duda and Hart.)

Interpolation and perceptual surface completion revisited

Recall this previous example of your visual system interpolating a smooth surface from sparse data (cf. Nakayama and Shimojo, 1992). The data is sparse because the information about depth comes from the disparities in the left and right vertical edges of the horizontal rectangle. The random dot stereogram illustrated here is “dense” rather than sparse. The data is dense because there were lots of potential features to match throughout the background and the square that floated out in depth. Here’s the example of a sparse random dot stereogram seen earlier:

GraphicsRow[{Image[],Image[]}]

One sees the white points on dark, central surface floating in front of a background surface, also with a white point texture. How can one model processes of interpolation? Earlier we approached the problem by constructing a cost function to be minimized.

Here’s another demonstration.  When you view the left and right images below in stereo, you may see a horizontal rectangle floating out in front of the vertical rectangle in the back:

backwidth=50;backheight=50;
{x0,y0,vwidth,vheight,vxoff,vyoff,hwidth,hheight,hxoff,hyoff}={(backwidth/2),(backheight/2),(backwidth/4),(backheight/8),(backwidth/2),(backheight/2),(backwidth/8),(backheight/4),(backwidth/2),(backheight/2)};

gleft=Show[Graphics[{Hue[0.6`],Rectangle[{0,0},{backwidth,backheight}],Hue[0.27`],Rectangle[{x0-vwidth,y0-vheight},{x0+vwidth,y0+vheight}],Hue[0.27`],Rectangle[{x0-hwidth,y0-hheight},{x0+hwidth,y0+hheight}]}],AspectRatio->Automatic];
gright=Show[Graphics[{Hue[0.6`],Rectangle[{0,0},{backwidth,backheight}],Hue[0.27`],Rectangle[{x0-vwidth-1,y0-vheight},{x0+vwidth-1,y0+vheight}],Hue[0.27`],Rectangle[{x0-hwidth,y0-hheight},{x0+hwidth,y0+hheight}]}],AspectRatio->Automatic];
gappearance=Show[Graphics[{Hue[0.6`],Rectangle[{0,0},{backwidth,backheight}],Hue[0.42`],Rectangle[{x0-hwidth,y0-hheight},{x0+hwidth,y0+hheight}],Hue[0.27`],Rectangle[{x0-vwidth,y0-vheight},{x0+vwidth,y0+vheight}]}],AspectRatio->Automatic,ImageSize->Small];
Show[GraphicsRow[{gright,gleft,"",gappearance}]]

If you can cross your eyes, so that the left image is in the right eye, and the right image in the left, you will see a green horizontal bar floating out in front of a green vertical bar. So-called "free-fusing" isn't easy, but when you've got it, you should see a total of three green crosses. The one in the middle is the one in which the two images are fused by your brain--and this is the one we are talking about. The interesting point here is that even though there is no local information in the image to support the percept of a horizontal occluding bar, observers still see an illusory completion.

It looks a bit like the figure on the far right, except that the color of the horizontal bar is changed here slightly just for illustration.

One’s first guess might be that observers should not see the horizontal bar as a plane, but as shown in the figure, rather they would interpolate a surface that on the left side is close to the viewer, but then descends back towards the depth of the vertical bar, stays at that constant far depth  left to right, and then comes back towards the viewer on the right (see Nakayama, K., & Shimojo, S., 1992.) This suggests a Bayesian prior that captures long-range, or global constraints, in addition to the local ones studied in this lecture.

Below we first review how to solve the interpolation problem using standard gradient descent on a cost function, and then re-formulate the problem on a graphical model to introduce  Bayesian Belief Propagation.

Review of interpolation using smoothness: Gradient descent

For simplicity, we'll assume 1-D as in the lecture on sculpting the energy function. In anticipation of formulating the problem in terms of a graph that represents conditional probability dependence, we represent observable depth cues by (SuperStar[y]), and the true ("hidden") depth estimates by y.

Figure from Weiss (1999).

We want to estimate the underlying function using only a local smoothness constraint. This is in contrast to linear regression, which assumes fitting the data with a straight line--a global constraint.

First-order smoothness

Earlier we saw that under specific assumptions, biologically plausible neural updating can be seen to decrease the value of an energy or cost function. One can also start off with an assumed cost function, determined by a set of constraints, and use gradient descent to derive an update rule that minimizes the cost.

For an interpolation problem, we can write the energy or cost function by:

where (Subscript[w, k])(= xs[[k]] below) is an "indicator function", and (Subscript[y, k]^( * ))= d, are the data values. The indicator function is 1 if where is data available, and zero otherwise. The second sum represents the sum of the squared differences between neighboring y-values. Minimizing the first sum encourages the estimates of y to be close to the measured values. Minimizing the second sum encourages nearby y-values to be the same. Thus minimizing J(Y) encourages fidelity to the data where present, and similarity to nearby values where there is no data. 

Gradient descent gives the following local update rule:

As before, λ is a free parameter that controls the degree of smoothness, i.e. smoothness at the expense of fidelity to the data.

At each step, the rule encourages the estimates to get closer to the data (where it exists), and also closer to the values of neighbors. 

There are various choices for how to change the smoothness as a function of iterations. 

E.g. standard methods in numerical analysis include: 
	Gauss-Seidel: η[k_] :=1 /(λ+xs[[k]]). And successive over-relaxation (SOR): η2[k_]:=1.9/(λ+xs[[k]]);

A simulation: Straight line with random missing data points

We look at an example where we have only two data points, both at the boundaries.

Make the data

Consider the problem of interpolating a set of points with missing data, marked by an indicator function with the following notation: 
(Subscript[w, k])= xs[[k]], y* = data, y=f.

We'll assume the true model is that f = y = j, where j=1 to size.  data is a function of the sampling process on  f = j

size=32;xs=Table[0,{i,1,size}];xs[[1]]=1;xs[[size]]=1;data=Table[N[j] xs[[j]],{j,1,size}];g3=ListPlot[Table[N[j],{j,1,size}],Joined->True,PlotStyle->{RGBColor[0,0.5,0]}];g2=ListPlot[data,Joined->False,PlotStyle->{Opacity[0.35],RGBColor[0.75,0.,0],PointSize[Large]}];

The green line shows the a straight line connecting the data points. The red dots on the abscissa mark the points where data are missing.

Show[{g2,g3},ImageSize->Medium]

The update rule is linear, so we can represent in terms of two matrices, Tm and Sm such that the gradient of the energy is equal to: 
Tm . f - Sm . data.  As before, Sm is the filter to exclude non-data points specified by the indicator function (Subscript[w, k])= xs[[k]] . Tm  expresses the "smoothness" constraint.

Sm = DiagonalMatrix[xs];
Tm = Table[0,{i,1,size},{j,1,size}];
For[i=1,i<=size,i++,Tm[[i,i]] = 2];
Tm[[1,1]]=1;Tm[[size,size]]=1; (*Adjust for the boundaries*)
For[i=1,i<size,i++, Tm[[i+1,i]] = -1];
For[i=1,i<size,i++, Tm[[i,i+1]] = -1];

Run gradient descent

Clear[f,d,λ]
(λ*Tm.Array[f,size]-Sm.((Array[d,size])-Array[f,size]))//MatrixForm ;

Clear[Tf,f1,j];
dt = 1; λ=2;
Tf[f1_] := f1 - dt*(1/(λ+xs))*(Tm.f1 - λ*Sm.(data-f1));

We will initialize the state vector to zero, and then run the network for iter iterations:

f0 = Table[0,{i,1,size}];
(*f0 = Table[RandomReal[{0,30}],{i,1,size}];*)
result=f0;
f=f0;
iter=25;

Now update the result after each iteration, and plot the interpolated function with each iteration step.

result=NestList[Tf,f,30*iter];

Manipulate[
g1=ListPlot[result[[j]],Joined->False,AspectRatio->Automatic,PlotRange->{{0,size},{-1,size+1}},ImageSize->Small],{j,1,Dimensions[result][[1]],1}]

Interpolation using Belief Propagation

Same interpolation problem, but now using belief propagation

Example is taken from Yair Weiss. (Weiss, 1999). The data variables are on top, and the hidden  variables to be estimated at the bottom of the graph.

Probabilistic generative model

data[[i]]=(SuperStar[y])[i]=xs[[i]] y[[i]] +dnoise, dnoise~N[0,(Subscript[σ, D])]

y[[i+1]]=y[[i]]+znoise, znoise~N[0,(Subscript[σ, R])]

The first term is the "data formation" model, i.e. how the data is directly influenced by the interaction of the underlying influences or causes: 

	y* is determined by an underlying hidden “y” which can’t be directly measured. But we assume we can measure y*, which is 	determined by sampling some values of y and adding noise. 
	
	The second term reflects our prior assumptions about the smoothness of y, i.e. nearby y's are correlated, and in fact assumed identical except for some added noise. So with no noise the prior reflects the assumption that lines are horizontal--all y's are the same. This is sometimes called a “soft constraint”, because it is a tendency--we don’t insist that it be satisfied exactly.

Summary of the message passing rules

We'd like to know the distribution of the random variables at each node i, conditioned on all the data: I.e. we want the posterior

	p((Subscript[y, i])=u | all the data) 

If we could find this, we'd be able to: 1) say what the most probable value of y is at each node, and 2) give a measure of confidence.

The appendix derives a solution, based on belief propagation, which sequentially updates local conditional distributions for the estimates of the mean of Y at each point, as well as estimating its standard deviation.

Let’s go over the rules of message passing. We follow Weiss, and make a (hopefully not too confusing) notation change to avoid the square superscripts, using the notation substitution: 

	(Subscript[σ, D]^2)→(Subscript[σ, D]), (Subscript[σ, R]^2)→(Subscript[σ, R]).

So  (Subscript[σ, D]) represents the variability or uncertainty in the data. If  (Subscript[σ, D]) is small, we trust the data (SuperStar[Y]).  (Subscript[σ, R])  represents our prior belief
in the “roughness” of the curve--the Y’s--to be estimated.  If (Subscript[σ, R]) is big, then we tolerate more deviations from straightness, while small values would bias the estimates towards smooth, straight lines. In other words, if the variance of the difference between neighboring estimates is small, then they would tend to have nearly the same values.

Let’s look at the messages passed for our example:

The Appendix derives a rule that tells us how to update the mean and variance parameters of the previous node, going right to left:

μα←(((Subscript[μα, p]) (Subscript[σ, D]^2)+(Subscript[y, p]) (Subscript[σα, p]^2))/((Subscript[σ, D]^2)+(Subscript[σα, p]^2)))=(((((Subscript[μα, p]) (Subscript[σ, D]^2))/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))+(((Subscript[y, p]) (Subscript[σα, p]^2))/((Subscript[σα, p]^2) (Subscript[σ, D]^2))))/(((Subscript[σ, D]^2)/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))+((Subscript[σα, p]^2)/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))))=((((Subscript[μα, p])/(Subscript[σα, p]^2))+((Subscript[y, p])/( (Subscript[σ, D]^2))))/((1/(Subscript[σα, p]^2))+(1/(Subscript[σα, p]^2))))

The update rule for the variance is:

((σα)^2)←(Subscript[σ, R]^2)+(1/((1/(Subscript[σ, D]^2))+(1/(Subscript[σα, p]^2))))

And by symmetry, we have a similar pair of update rules for the messages μβ and ((σβ)^2) passed from left to right (see Appendix).

Recall that sometimes we have data and sometimes we don’t. So replace:
			(Subscript[y, p])→(Subscript[w, i-1])(Subscript[y, i-1]^( * ))

Side note: The underlying operation in which estimates get combined weighted by uncertainty,  is similar to how we updated the mean in the Bayesian learning example. This is a recurring theme. In your assignment, you show that the maximum a posteriori estimate for cue combination is:

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

where (Subscript[r, i]) = 1/(Subscript[σ, i]^2), and (Subscript[μ, 1]) and (Subscript[μ, 2]) are the estimates of the means obtained separately. How this could be done neurally has inspired recent work in neural population codes, mentioned below. Combining information weighted by reliability underlies a number of models of visual and multimodal cue integration in psychophysics over the past decade.

A simulation: Belief propagation for interpolation with missing data

Here’s an implementation of the Bayesian belief updates applied to the above interpolation problem.  We change some of the notation for the program, with equivalent symbols:

(Subscript[y, p])→=(Subscript[w, i-1])(Subscript[y, i-1]^( * )) = xs[i-1]data[i-1]

y→ yfit

size=32;xs=Table[0,{i,1,size}];xs[[1]]=1;xs[[size]]=1;data=Table[N[j] xs[[j]],{j,1,size}];g3bp=ListPlot[Table[N[j],{j,1,size}],Joined->True,PlotStyle->{RGBColor[0,0.5,0]}];g2bp=ListPlot[data,Joined->False,PlotStyle->{Opacity[0.35],RGBColor[0.75,0.,0],PointSize[Large]}];

The green line shows the a straight line connecting the data points. The red dots on the abscissa mark the points where data are missing.

Show[{g2bp,g3bp},ImageSize->Medium]

Initialization

size=32;
μ0=1; 
μα=1;σα=100000; (*large uncertainty *)
μβ=1;σβ=100000;(*large*)
σR=4.0;σD=1.0;
µ=Table[μ0,{i,1,size}];
σ=Table[σα,{i,1,size}];
µα=Table[μ0,{i,1,size}];
σα=Table[σα,{i,1,size}];
µβ=Table[μ0,{i,1,size}];
σβ=Table[σβ,{i,1,size}];
iter=0;
i=1;j=size;

The code below implements the above iterative equations, taking care near the boundaries. The plot shows the estimates of (Subscript[y, i])=  µ, and the error bars show ±(Subscript[σ, i]).

Belief Propagation Routine

yfit=Table[{0,0},{i1,1,size}];
g1b=ErrorListPlot[{yfit}];
Dynamic[Show[{g1b,g2bp,g3bp,Graphics[{Text["Iteration="<>ToString[iter],{(size/2),size}]}]},PlotRange->{-50,50},Axes->{False,True}]]

Execute the next cell to run 31 iterations. The display is slowed down so that you can see the progression of the updates in the above graph.

Do[
Pause[.5];
µ[[i]]=((((xs[[i]] data[[i]])/(σD))+((µα[[i]])/(σα[[i]]))+((1. µβ[[i]])/(σβ[[i]])))/(((xs[[i]])/(σD))+(1/(σα[[i]]))+(1/(σβ[[i]]))));
σ[[i]]=(1./(((xs[[i]])/(σD))+(1/(σα[[i]]))+(1/(σβ[[i]]))));
µ[[j]]=((((xs[[j]] data[[j]])/(σD))+((µα[[j]])/(σα[[j]]))+((1. µβ[[j]])/(σβ[[j]])))/(((xs[[j]])/(σD))+(1/(σα[[j]]))+(1/(σβ[[j]]))));
σ[[j]]=(1./(((xs[[j]])/(σD))+(1/(σα[[j]]))+(1/(σβ[[j]]))));
nextj=j-1;
µα[[nextj]]=((((xs[[j]] data[[j]])/(σD))+((1. µα[[j]])/(σα[[j]])))/(((xs[[j]])/(σD))+(1/(σα[[j]]))));
σα[[nextj]]=σR+(1./(((xs[[j]])/(σD))+(1/(σα[[j]]))));
nexti=i+1;
µβ[[nexti]]=((((xs[[i]] data[[i]])/(σD))+((1. µβ[[i]])/(σβ[[i]])))/(((xs[[i]])/(σD))+(1/(σβ[[i]]))));
σβ[[nexti]]=σR+(1./(((xs[[i]])/(σD))+(1/(σβ[[i]]))));
j--;
i++;
iter++;
yfit=Table[{µ[[i1]],σ[[i1]]},{i1,1,size}];
g1b=ErrorListPlot[{yfit}];
,{size-1}];

Relation to neural networks

We’ve seen two very different ways of estimating states through iterative updating. It is easy to see how the first, derived from gradient descent, is related to updating in a traditional neural network. The second, belief propagation, gives us a very different view of how information might be represented and updated in a network.  For example, in the second case applied to interpolation, the nodes represent probability distributions over depths (summarized by the mean and variance) not just estimates of depth. A critical new feature is the explicit representation of uncertainty, i.e. in the standard deviations of the node values. Later in the course, we will discuss whether the brain might have explicit representations of such uncertainties, how these might be represented in populations of neurons, and how the brain might do computations on these. Does the brain do belief propagation, and if so how?

For a preview of the general question of whether the brain represents and computes on distributions, see: Ma et al. (2006, 2012), Knill & Pouget (2004), Zemel & Pouget (1998).

Discussion: relation to behavior?

Yair Weiss and colleagues showed how the probablistic combination of local motion cues and prior motion assumptions could explain a number of visual illusions: Weiss, Y., Simoncelli, E. P., & Adelson, E. H. (2002). Motion illusions as optimal percepts. Nature Neuroscience.

But their conclusion (and probably all other similar applications to perceptual behavior) didn’t depend on the algorithm.

Can you think of an application and test of belief propagation to perception?

More examples and exercises

Generate sine wave with missing data

size=64;xs=Table[RandomInteger[1],{i,1,size}];data=Table[N[Sin[((2 π j)/20)] xs[[j]]],{j,1,size}];g3b=ListPlot[Table[N[Sin[((2 π j)/20)]],{j,1,size}],Joined->True,PlotStyle->{RGBColor[0,0.5,0]}];g2b=ListPlot[data,Joined->False,PlotStyle->{RGBColor[0.75,0.,0]}];

Initialize

μ0=1; 
μα=1;σα=100000; (*large uncertainty *)
μβ=1;σβ=100000;(*large*)
σR=.5;σD=.1;
µ=Table[μ0,{i,1,size}];
σ=Table[σα,{i,1,size}];
µα=Table[μ0,{i,1,size}];
σα=Table[σα,{i,1,size}];
µβ=Table[μ0,{i,1,size}];
σβ=Table[σβ,{i,1,size}];
iter=0;
i=1;j=size;

yfit=Table[{0,0},{i1,1,size}];
g1bb=ErrorListPlot[{yfit}];
Dynamic[Show[{g1bb,g2b,g3b},PlotRange->{-2,2},Axes->{False,True}]]

SINE WAVE DEMO: Belief Propagation Routine

Do[
Pause[0.2];
µ[[i]]=((((xs[[i]] data[[i]])/(σD))+((µα[[i]])/(σα[[i]]))+((1. µβ[[i]])/(σβ[[i]])))/(((xs[[i]])/(σD))+(1/(σα[[i]]))+(1/(σβ[[i]]))));
σ[[i]]=(1./(((xs[[i]])/(σD))+(1/(σα[[i]]))+(1/(σβ[[i]]))));
µ[[j]]=((((xs[[j]] data[[j]])/(σD))+((µα[[j]])/(σα[[j]]))+((1. µβ[[j]])/(σβ[[j]])))/(((xs[[j]])/(σD))+(1/(σα[[j]]))+(1/(σβ[[j]]))));
σ[[j]]=(1./(((xs[[j]])/(σD))+(1/(σα[[j]]))+(1/(σβ[[j]]))));
nextj=j-1;
µα[[nextj]]=((((xs[[j]] data[[j]])/(σD))+((1. µα[[j]])/(σα[[j]])))/(((xs[[j]])/(σD))+(1/(σα[[j]]))));
σα[[nextj]]=σR+(1./(((xs[[j]])/(σD))+(1/(σα[[j]]))));
nexti=i+1;
µβ[[nexti]]=((((xs[[i]] data[[i]])/(σD))+((1. µβ[[i]])/(σβ[[i]])))/(((xs[[i]])/(σD))+(1/(σβ[[i]]))));
σβ[[nexti]]=σR+(1./(((xs[[i]])/(σD))+(1/(σβ[[i]]))));
j--;
i++;
iter++;
yfit=Table[{µ[[i1]],σ[[i1]]},{i1,1,size}];
g1bb=ErrorListPlot[{yfit}];
,{size-1}]

Final projects

Outline due November 21.

Appendix

Derivation of the update rules, called “message passing”

Updating the mean and variance given the data at point i, and current beliefs about mean and variance before and after i

Let p((Subscript[Y, i])=u | all the data) be normally distributed: NormalDistribution[μi,σi].

Consider the ith unit. The posterior p((Subscript[Y, i])=u | all the data) =

p((Subscript[Y, i])=u| all the data) ∝ p((Subscript[Y, i])=u | data before i) p(data at i | (Subscript[Y, i])=u) p((Subscript[Y, i])=u | data after  i)

(“before” and “after” means to the left and right of i, respectively.)
 
 Suppose that p((Subscript[Y, i])=u | data before i) is also gaussian:  
						p((Subscript[Y, i])=u | data before i) = (Subscript[α, i])[u] ∼ NormalDistribution[μα, σα]
						
and so is probability conditioned on the data after i:

						p((Subscript[Y, i])=u | data after  i)= (Subscript[β, i])[u] ∼ NormalDistribution[μβ, σβ]

And the noise model for the data: 
						p(data at i | (Subscript[Y, i])=u) = (Subscript[L, i])[u]∼ NormalDistribution[yp,(Subscript[σ, D])]
						yp=data[[i]]

So in terms of these functions, the posterior probability of the ith unit taking on the value u can be expressed as proportional to a product of the three factors:

p((Subscript[Y, i])=u | all the data) ∝ (Subscript[α, i])[u]*(Subscript[L, i])[u]*(Subscript[β, i])[u]

ClearAll[μα,σα, μβ,σβ,yp,σD];

αudist=NormalDistribution[μα,σα];
α[u]=PDF[αudist,u];

Ddist=NormalDistribution[yp,σD];
L[u]=PDF[Ddist,u];

βudist=NormalDistribution[μβ,σβ];
β[u]=PDF[βudist,u];

α[u]*L[u]*β[u]

This just another gaussian distribution on (Subscript[Y, i])=u. What is its mean and variance? Finding the root enables us to complete the square to see what the numerator looks like. In particular, what the mode (=mean for gaussian) is.

Solve[-D[-((((u-μα))^2)/(2 ((σα)^2)))-((((u-μβ))^2)/(2 ((σβ)^2)))-((((u-(Subscript[y, p])))^2)/(2 (Subscript[σ, D]^2))),u]==0,u]

This suggests that if we had estimates of μα, μβ, ((σα)^2), ((σβ)^2) and the data, we could update the mean of node i using:

u ←  ((((μα)/((σα)^2))+((μβ)/((σβ)^2))+((Subscript[y, p])/(Subscript[σ, D]^2)))/((1/((σα)^2))+(1/((σβ)^2))+(1/(Subscript[σ, D]^2))))

Similarly, the update rule for the variance is:

(σ^2)← (1/((σα)^2))+(1/((σβ)^2))+(1/(Subscript[σ, D]^2))

How do we get μα, μβ , σα, σβ?

We express the probability of the ith unit taking on the value u in terms of the values of the neighbor  before, conditioning on what is known (the observed measurements), and marginalizing over what isn't (the previous "hidden" node value, v, at the i-1th location).

We have three terms to worry about that depend on nodes in the neighborhood preceding i:

α[u]=(∫_(-∞)^∞)(Subscript[α, p])[v]*S[u]*L[v]dv∝(∫_(-∞)^∞)(e^(-((((v-(Subscript[y, p])))^2)/(2 (Subscript[σ, D]^2)))-((((u-v))^2)/(2 (Subscript[σ, R]^2)))-((((v-(Subscript[μα, p])))^2)/(2 (Subscript[σα, p]^2)))))dv

(Subscript[α, p])=(Subscript[α, i-1]). S[u] is our smoothing term, or transition probability: S[u]=p(u|v). 
L[v] is the likelihood of v given the data previous at the previous  node.

Rdist=NormalDistribution[v,(Subscript[σ, R])];
S[u]=PDF[Rdist,u];

αvdist=NormalDistribution[(Subscript[μα, p]),(Subscript[σα, p])];
(Subscript[α, p])[v]=PDF[αvdist,v];

Lp[v]=PDF[Ddist,v];

Integrate[(Subscript[α, p])[v]*S[u]*Lp[v],{v,-Infinity,Infinity}]

Some uninspired Mathematica manipulations

Let’s find an expression for the mode of the above calculated expression for α[u]

D[-((((u-(Subscript[μα, p])))^2) (Subscript[σ, D]^2)+(Subscript[μα, p]^2) (Subscript[σ, R]^2)+(u^2) (Subscript[σα, p]^2)+(Subscript[y, p]^2) ((Subscript[σ, R]^2)+(Subscript[σα, p]^2))-2 (Subscript[y, p]) ((Subscript[μα, p]) (Subscript[σ, R]^2)+u (Subscript[σα, p]^2)))/(2 ((Subscript[σ, R]^2) (Subscript[σα, p]^2)+(Subscript[σ, D]^2) ((Subscript[σ, R]^2)+(Subscript[σα, p]^2)))),u]

Solve[-%==0,u]

Simplify[((((Subscript[μα, p]) (Subscript[σ, D]^2))/((Subscript[σ, R]^2) (Subscript[σα, p]^2)+(Subscript[σ, D]^2) ((Subscript[σ, R]^2)+(Subscript[σα, p]^2))))+(((Subscript[y, p]) (Subscript[σα, p]^2))/((Subscript[σ, R]^2) (Subscript[σα, p]^2)+(Subscript[σ, D]^2) ((Subscript[σ, R]^2)+(Subscript[σα, p]^2)))))/((Subscript[σ, D]^2)*(Subscript[σα, p]^2))]

Simplify[(((Subscript[σ, D]^2)/((Subscript[σ, R]^2) (Subscript[σα, p]^2)+(Subscript[σ, D]^2) ((Subscript[σ, R]^2)+(Subscript[σα, p]^2))))+((Subscript[σα, p]^2)/((Subscript[σ, R]^2) (Subscript[σα, p]^2)+(Subscript[σ, D]^2) ((Subscript[σ, R]^2)+(Subscript[σα, p]^2)))))/((Subscript[σ, D]^2)*(Subscript[σα, p]^2))]

((((Subscript[μα, p]) (Subscript[σ, D]^2)+(Subscript[y, p]) (Subscript[σα, p]^2))/((Subscript[σ, D]^2) (Subscript[σ, R]^2) (Subscript[σα, p]^4)+(Subscript[σ, D]^4) (Subscript[σα, p]^2) ((Subscript[σ, R]^2)+(Subscript[σα, p]^2)))))/((((Subscript[σ, D]^2)+(Subscript[σα, p]^2))/((Subscript[σ, D]^2) (Subscript[σ, R]^2) (Subscript[σα, p]^4)+(Subscript[σ, D]^4) (Subscript[σα, p]^2) ((Subscript[σ, R]^2)+(Subscript[σα, p]^2)))))

We now have a rule that tells us how to update the α(u)=p(yi=u|data before i), in terms of the mean and variance parameters of the previous node:

μα←(((Subscript[μα, p]) (Subscript[σ, D]^2)+(Subscript[y, p]) (Subscript[σα, p]^2))/((Subscript[σ, D]^2)+(Subscript[σα, p]^2)))=(((((Subscript[μα, p]) (Subscript[σ, D]^2))/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))+(((Subscript[y, p]) (Subscript[σα, p]^2))/((Subscript[σα, p]^2) (Subscript[σ, D]^2))))/(((Subscript[σ, D]^2)/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))+((Subscript[σα, p]^2)/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))))=((((Subscript[μα, p])/(Subscript[σα, p]^2))+((Subscript[y, p])/( (Subscript[σ, D]^2))))/((1/(Subscript[σα, p]^2))+(1/(Subscript[σα, p]^2))))

The update rule for the variance is:

((σα)^2)←(Subscript[σ, R]^2)+(1/((1/(Subscript[σ, D]^2))+(1/(Subscript[σα, p]^2))))

A similar derivation gives us the rules for μβ, ((σβ)^2)

μβ ← ((((Subscript[μβ, a])/(Subscript[σβ, a]^2))+((Subscript[y, a])/( (Subscript[σ, D]^2))))/((1/(Subscript[σβ, a]^2))+(1/(Subscript[σα, a]^2))))

((σβ)^2)←(Subscript[σ, R]^2)+(1/((1/(Subscript[σ, D]^2))+(1/(Subscript[σβ, a]^2))))

Where the subscript index  p (for "previous", i.e. unit i-1) is replaced by a (for "after", i.e. unit i+1).

Recall that sometimes we have data and sometimes we don't. So replace:

(Subscript[y, p])→xs[i-1]data[i-1]=(Subscript[w, i-1])(Subscript[y, i-1]^( * ))

And similarly for (Subscript[y, a]).

The ratio, (((((Subscript[σ, D])/(Subscript[σ, R]))))^2)plays the role of λ above. If (Subscript[σ, D]^2)>>(Subscript[σ, R]^2), there is greater smoothing. If (Subscript[σ, D]^2)<<(Subscript[σ, R]^2), there is more fidelity to the data. (Recall (SuperStar[y])->data.(Subscript[w, k])->xs[[k]]). But now we have a principled way of assigning the relative amount of smoothing.

References

Applebaum, D. (1996). Probability and Information . Cambridge, UK: Cambridge University Press.
Frey, B. J. (1998). Graphical Models for Machine Learning and Digital Communication. Cambridge, Massachusetts: MIT Press.
Jepson, A., & Black, M. J. (1993). Mixture models for optical flow computation. Paper presented at the Proc. IEEE Conf. Comput. Vsion Pattern Recog., New York.
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., & Madarasmi, S. (1995). The Visual Perception of Surfaces, their Properties, and Relationships. DIMACS Series in Discrete Mathematics and Theoretical Computer Science, 19, 373-389.
Knill, D. C., & Pouget, A. (2004). The Bayesian brain: the role of uncertainty in neural coding and computation. Trends in Neurosciences, 27(12), 712–719. doi:10.1016/j.tins.2004.10.007
Ma, W. J., Beck, J. M., Latham, P. E., & Pouget, A. (2006). Bayesian inference with probabilistic population codes. Nature Neuroscience, 9(11), 1432–1438. doi:10.1038/nn1790
Ma, W. J. (2012). Organizing probabilistic models of perception. Trends in Cognitive Sciences, 16(10), 511–518. doi:10.1016/j.tics.2012.08.010
Madarasmi, S., Kersten, D., & Pong, T.-C. (1993). The computation of stereo disparity for transparent and for opaque surfaces. In C. L. Giles & S. J. Hanson & J. D. Cowan (Eds.), Advances in Neural Information Processing Systems 5. San Mateo, CA: Morgan Kaufmann Publishers.
Nakayama, K., & Shimojo, S. (1992). Experiencing and perceiving visual surfaces. Science, 257(5075), 1357–1363. doi:10.1126/science.1529336
Pearl, Judea. (1997) Probabilistic Reasoning in Intelligent Systems : Networks of Plausible Inference. (amazon.com link)
Ripley, B. D. (1996). Pattern Reco gnition and Neural Networks. Cambridge, UK: Cambridge University Press.
Weiss Y. (1999) Bayesian Belief Propagation for Image Understanding submitted to SCTV 1999. (gzipped postscript 297K) 
Weiss, Y. (1997). Smoothness in Layers: Motion segmentation using nonparametric mixture estimation. Paper presented at the Proceedings of IEEE conference on Computer Vision and Pattern Recognition.
Weiss, Y., Simoncelli, E. P., & Adelson, E. H. (2002). Motion illusions as optimal percepts. Nature Neuroscience, 5(6), 598–604. doi:10.1038/nn858
Yuille, A., Coughlan J., Kersten D.(1998) (pdf)
Zemel, R. S. & Pouget, A. (1998). Probabilistic interpretation of population codes. Neural Computation, 10(2), 403–430.

For notes on Graphical Models, see:http://www.cs.berkeley.edu/~murphyk/Bayes/bayes.html