Bayesian inference in the Wolfram Language

Using an implementation of the nested sampling algorithm

dr. ir. Sjoerd Smit
Technical consultant a Wolfram Research Europe Ltd.

WTC 2017

Brief introduction to Bayesian inference

Bayes’ rule

Bayes’ rule (or theorem) is a fundamental mathematical result that allows you to do calculations with conditional probabilities. A conditional probability is written as:

\[ProbabilityPr][X\[Conditioned]Y] (pronounced: "the probability of X given Y")

It is calculated by using the multiplication rule for probabilities:

\[ProbabilityPr][X&&Y]== \[ProbabilityPr][X] × \[ProbabilityPr][Y\[Conditioned]X]

or, after switching X and Y around (since X && Y is just the same as Y && X):

\[ProbabilityPr][X&&Y]== \[ProbabilityPr][Y] × \[ProbabilityPr][X\[Conditioned]Y]

(Note: these formulas are ALWAYS valid. Not just when X and Y are independent)

You get Baye’s rule by simply eliminating the probability for X && Y:

\[ProbabilityPr][X\[Conditioned]Y]==((\[ProbabilityPr][Y\[Conditioned]X]\[ProbabilityPr][X])/(\[ProbabilityPr][Y]))

In other words: Baye’s rule reverses the conditioning of probabilities

Heads and tails: a simple example of inference

Suppose you have a limited amount of data from a random binary phenomenon and you want to say something about how biased the process might be:

binaryData
Counts[binaryData]

We will assume that each event has a probability p of giving True and 1-p of giving False. From a Bayesian perspective, the parameter p is an unknown, so it has a probability distribution.

Thus, the question is: what is the probability distribution for p given the fact that we have seen these data?

\[ProbabilityPr][p\[Conditioned]binaryData, Model]==((\[ProbabilityPr][binaryData\[Conditioned]p, Model]\[ProbabilityPr][p\[Conditioned]Model])/(\[ProbabilityPr][binaryData\[Conditioned]Model]))

\[ProbabilityPr][binaryData\[Conditioned]p, Model] is called the likelihood and can be calculated directly from the BinomialDistribution:

likelihood=Likelihood[
BinomialDistribution[Length[binaryData],\[FormalP]],
{Count[binaryData, True]}
]

The factor \[ProbabilityPr][p\[Conditioned]Model] is called the prior. Let’s assume that (without seeing the data) p \[Distributed] UniformDistribution[{0, 1}] to express that we have no idea what value p has.

prior = PDF[UniformDistribution[{0,1}],\[FormalP]]

Finally, the normalizing factor \[ProbabilityPr][coinData\[Conditioned]Model] is called the evidence and its main use is to compare different statistical descriptions of the data: the model with the highest evidence is generally preferred. It is calculated by evaluating the integral:

evidence == ∫\[ProbabilityPr][binaryData\[Conditioned]p, Model]\[ProbabilityPr][p\[Conditioned]Model]dp == ∫likelihood[p]prior[p]dp

We can plug all this into Bayes’ rule to get the posterior:

evidence = Integrate[likelihood*prior,{\[FormalP],0,1}]

posterior = ProbabilityDistribution[likelihood*prior/evidence,{\[FormalP],0,1}]

Plot[
Evaluate[{
PDF[UniformDistribution[{0,1}],p],
PDF[posterior,p]
}],
{p,0,1},
PlotRange->All,Filling->Axis,AxesLabel->{"(Subscript[p, True])","PDF[(Subscript[p, True])]"},PlotLegends->{"Prior", "Posterior"}
]

As far as Bayesian inference is concerned, this PDF represents the full extent of your knowledge of p. You can use the posterior to answer different questions about your state of knowledge of p:

Mean[posterior]

NumberLinePlot[
{Interval@InverseCDF[posterior,{0.025,0.975}], InverseCDF[posterior,0.5]},
PlotRange->{0,1},PlotLabel->Style["95% credible interval and median","Text"]] 

Probability[\[FormalP]< 0.5,\[FormalP] \[Distributed]posterior]

The nested sampling algorithm

The difficult part of doing inference

As we saw, in order to use Bayes’ rule we need to evaluate the integral:

evidence == ∫\[ProbabilityPr][Data\[Conditioned](UnderBar[θ]), Model]\[ProbabilityPr][(UnderBar[θ])\[Conditioned]Model]d(UnderBar[θ]) == ∫likelihood[(UnderBar[θ])]prior[(UnderBar[θ])]d(UnderBar[θ])

and for most practical problems this is a difficult N-dimensional beast that you cannot do reliably with ordinary analytical or numerical methods.

Solution 1: avoid computing the evidence (e.g., MCMC methods)

Many Bayesian methods avoid the computation of the evidence by using a self-normalizing algorithm. E.g., the Markov Chain Monte Carlo (MCMC) method allows you to draw samples from the posterior distribution without knowing the normalizing factor.
The Metropolis-Hastings algorithm is the most famous example of this and many different MCMC methods have been derived from it.

Solution 2: sample the prior in a smart way to estimate the evidence

How to do nested sampling

BlockRandom[
SeedRandom[7];
dist=CauchyDistribution[0,1];
testdata=RandomVariate[dist,10];
nLiveSamples = 10;
result = nestedSampling[
LogLikelihood[CauchyDistribution[#1, #2],testdata]&,
{"LocationParameter", "ScaleParameter"},
{
{loc,-10,10},
{scale,0.01,5}
},
"SamplePoolSize" -> nLiveSamples,
"MonteCarloSteps"->300,
"InitialMonteCarloStepSize" -> 0.1
];
];
NumberLinePlot[Tooltip[#,#]&/@testdata, PlotLabel->Style["10 Samples from CauchyDistribution[0, 1]","Text"], ImageSize->800, PlotStyle->PointSize[0.01]]
Dataset[result]

ContourPlot[
LogLikelihood[CauchyDistribution[loc,scale],testdata],
{loc,-10,10},{scale,0.01,5},
FrameLabel->(Style[#,"Text",16]&)/@{"location parameter","scale parameter"},ImageSize->Large,
PlotLabel->Style["LogLikelihood[loc, scale]","Text",16]
]

nestedSamplingManipulate[result,testdata]

posteriorBubbleChart[
result,
{loc,scale},
1,(* fraction of posterior mass to plot *)
ScalingFunctions->{None,"Log",None}]

posteriorMarginalPDFDensityPlot2D[result,{loc,scale},ScalingFunctions->{None,"Log",None},PlotPoints->50]

AssociationThread[result["Parameters"],Mean@result["EmpiricalPosteriorDistribution"]]
covarianceMatrixPlot[result,PlotLabel->Style["Covariance matrix of posterior distribution", "Text"]]

LogPlot[
Evaluate[{
PDF[predictiveDistribution[result],x],
PDF[EstimatedDistribution[testdata,result["GeneratingDistribution"][\[FormalX],\[FormalY]],ParameterEstimator->{"MaximumLikelihood", Method -> "FindMaximum"}], x],
PDF[dist,x]
}],
{x,-25,25},
PlotRange->All,
PlotLegends->{"Posterior predictive PDF", "PDF based on maximum likelihood estimate", "Original distribution"},
ImageSize->600,
Filling->Axis
]

Model comparison

Fitting Normal data with a SkewNormal distribution

DynamicModule[{stdSkewNormalDistribution},
Manipulate[
Column[{
Style["Comparison of SkewNormalDistribution with 0 mean and unit standard deviation","Text",16],
Plot[
Evaluate[{
PDF[NormalDistribution[0,1],x],
PDF[stdSkewNormalDistribution[α],x]
}],
{x,-3,3},
Filling->Axis,
PlotLegends->{NormalDistribution[0,1], stdSkewNormalDistribution[α]},
ImageSize->600
]
}],
{{α,0,"Skew paramater"},-5,5,0.1, Appearance->"Open"}
],
Initialization:>(stdSkewNormalDistribution[a_] :=N@SkewNormalDistribution[-(((Sqrt[2]) a)/((Sqrt[1+(a^2)]) (Sqrt[-2+(2/(1+(a^2)))+π]))), ((Sqrt[π])/(Sqrt[-2+(2/(1+(a^2)))+π])),a])
]

dist=NormalDistribution[0,1];
testdata=RandomVariate[dist,200];
result=<||>;
Histogram[testdata,50]

dist=SkewNormalDistribution[-1.1583129633811582`,1.5302577956464851`,3.`];
testdata=RandomVariate[dist,200];
result=<||>;
Histogram[testdata,50]

result["Normal"] = parallelNestedSampling[
LogLikelihood[NormalDistribution[##],testdata]&,
{"LocationParameter", "ScaleParameter"},
{{loc,-20,20},
{scale,0.01,100}},
"SamplePoolSize"->25
];
PrintTemporary["Normal model done..."];
result["SkewNormal"] = parallelNestedSampling[
LogLikelihood[SkewNormalDistribution[##],testdata]&,
{"LocationParameter", "ScaleParameter","LocationParameter"},
{{loc,-20,20},
{scale,0.01,100},
{skew, -20,20}}
];
Dataset[result[[All,{"LogEvidence"}]]]
Dataset[result[[All,{"ParameterValues"}]]]

FindDistribution[testdata]

posteriorBubbleChart[
result["Normal"],
{loc,scale},
0.99,
ScalingFunctions->{None,"Log",None}]

posteriorBubbleChart[
result["SkewNormal"],
{loc,scale},
0.99
]

posteriorBubbleChart[
result["SkewNormal"],
{loc,scale,skew},
0.99
]

skewManipulate[testdata]

Plot[
Evaluate[{
PDF[predictiveDistribution[result[["Normal"]],NormalDistribution,0.99],x],
PDF[predictiveDistribution[result[["SkewNormal"]],SkewNormalDistribution,0.99],x],
PDF[dist,x]
}],
{x,-4,4},
PlotRange->All,
Filling->Axis,
PlotLegends->{"Posterior PDF (Normal)", "Posterior PDF (SkewNormal)", "Original distribution"},
PlotRange->All,
ImageSize->Large
]

Bayesian polynomial regression

BlockRandom[
SeedRandom[1];
testdata={#,(#+2 RandomReal[])}&/@RandomReal[{-1,1},20];
ListPlot[testdata]
]

result=
Association@Table[
order-> With[{vars={
Sequence@@Table[{Symbol["x"<>ToString[i]],-100,100},{i,0,order}],
{s0,0.01,100}
}
},
gaussianProcessNestedSampling[
Rule@@@testdata,
Function[0.],
Function[s0],
Function[Evaluate[Total@Table[vars[[i,1]] ((#)^(i-1)),{i,1,Length[vars]-1}]]],
Flatten@{ConstantArray["LocationParameter",Length[vars]-1], "ScaleParameter"},
vars,
"Likelihood" -> "Custom",
"MonteCarloSteps"->300,
"SamplePoolSize"->25,
"ParallelRuns" ->4
]
],
{order,{1,2,3}}];
Dataset[result[[All,{"LogEvidence"}]]]

predictions = predictFromGaussianProcess[#, Range[-2,2,0.2],0.99]&/@result;
Column[Table[
gaussianProcessPredictionPlot[result[[i]], predictions[[i]], "DistributionPercentiles" -> {0.01,0.5,0.99},ImageSize->400,PlotLegends->{"(1^st) percentile","Median","(99^th) percentile"}],
{i,Length[result]}
]]

References:

E. T. Jaynes, Probability Theory - The Logic of Science, Cambridge University Press (2003)

John Skilling, Nested sampling for general Bayesian computation, Bayesian Anal. Volume 1, Number 4 (2006), 833-859, DOI:10.1214/06-BA127

David J.C. MacKay, Information Theory, Inference, and Learning Algorithms, Cambridge University Press (2003). Freely available at http://www.inference.org.uk/mackay/itila/

Github repository:

Hyperlink[Style["Bayesian Inference in Wolfram Language","Text"],"https://github.com/ssmit1986/BayesianInference"]

Initialization

$HistoryLength=10;
<<BayesianStatistics`;
binaryData={True,True,True,False,True,True,True,True,False,True,True,False};

nestedSamplingManipulate[result_,testdata_]:=DynamicModule[{
samples=Values/@Values@Part[
SortBy[#LogLikelihood&]@result["Samples"],
All,{"Point","LogLikelihood"}],
subset,
nLiveSamples=result["SamplePoolSize"],
dist = result["GeneratingDistribution"]
},
Manipulate[
subset=TakeDrop[Take[samples,UpTo[n]],-nLiveSamples];
ContourPlot[
LogLikelihood[dist[loc,scale],testdata],
{loc,-10,10},{scale,0.01,5},
FrameLabel->(Style[#,"Text",16]&)/@{"location parameter","scale parameter"},ImageSize->Large,
Epilog->{Dynamic[
{
{GrayLevel[0.7],PointSize[0.0075],Point[samples[[All,1]]]},
{Black,PointSize[0.01],Point[Most@subset[[2,All,1]]],PointSize[0.025],Point[Last@subset[[2,All,1]]]},
{Red,PointSize[0.01], Point[Most@subset[[1,All,1]]],PointSize[0.025],Point[Last@subset[[1,All,1]]]}

}]},
Contours->subset[[2,All,2]],
PlotLegends->Grid[{{Graphics[{Red,Disk[]},ImageSize->10],"Alive"},{Graphics[{Black,Disk[]},ImageSize->10],"Dead"},
{Graphics[{GrayLevel[0.7],Disk[]},ImageSize->10],"Unborn"}},BaseStyle->"Text",Alignment->Left]
],
{n,nLiveSamples+1,Min[Length[samples],nLiveSamples+20],1, Appearance->"Open"}
]
];

skewManipulate[testdata_]:=DynamicModule[{
dist, scale,skew
},
Manipulate[
dist=EstimatedDistribution[testdata,SkewNormalDistribution[loc,scale,skew]];
Show[
Histogram[testdata,50,"PDF",PlotLabel->dist,Epilog->{InfiniteLine[{{loc,0},{loc,1}}]}],
Plot[Evaluate@PDF[dist,x],{x,-3,3},Filling->Axis],
PlotRange->All,
ImageSize->500
],
{{loc,0,"Location parameter"},-2,2}
]
];