Essays, Posts & Presentations

Embracing Uncertainty: Better Model Selection with Bayesian Linear Regression

Sjoerd Smit

Author

Sjoerd Smit

Title

Embracing Uncertainty: Better Model Selection with Bayesian Linear Regression

Description

Application and interpretation of Bayesian regression in the Wolfram Language compared to regular data-fitting algorithms. Plus a link to download the BayesianLinearRegression function used.

Embracing Uncertainty: Better Model Selection with Bayesian Linear Regression

August 22, 2019

Sjoerd Smit, Technical Consultant, European Sales

Readers who follow the Mathematica Stack Exchange (which I highly recommend to any Wolfram Language user) may have seen this post recently, in which I showed a function I wrote to make Bayesian linear regression easy to do. After finishing that function, I have been playing around with it to get a better feel of what it can do, and how it compares against regular fitting algorithms such as those used by

Fit

. In this blog post, I don’t want to focus too much on the underlying technicalities (check out my previous blog post to learn more about Bayesian neural network regression); rather, I will show you some of the practical applications and interpretations of Bayesian regression, and share some of the surprising results you can get from it.

Getting Started The easiest way to get my BayesianLinearRegression function is through my submission to the Wolfram Function Repository. To use this version of the function with the code examples in this blog, you can evaluate the following line, which creates a shortcut for the resource function:

In[]:=

BayesianLinearRegression = ResourceFunction["BayesianLinearRegression"]

Out[]=

[◼]	BayesianLinearRegression

You can also visit the GitHub repository and follow the installation instructions to load the BayesianInference package with the following:

In[]:=

<<BayesianInference`

Alternatively, you can get the standalone source file of BayesianLinearRegression and evaluate it to make the function available to you, though the function regressionPlot1D I use later will not work in the following examples if you don’t get the full BayesianInference package. Its definition is in the BayesianVisualisations.wl file.

Back to Basics I want to do something that will be very familiar to most people with some background in data fitting: polynomial regression. I could have picked something more complicated, but it turns out that the Bayesian take on fitting data adds quite a bit of depth and new possibilities even to something as simple as fitting polynomials, making it an ideal demonstration example.

I picked the example data here specifically because it hints at a straight trend, while still leaving some doubt:

In[]:=

data=

;plot=ListPlot[data,PlotStyleRed]

Out[]=

So let’s start off by doing the first sensible thing to do: fit a line through it. Concretely, I will fit the model

ya+bx+ϵ

, where

ϵ

NormalDistribution

[0,σ]. BayesianLinearRegression uses the same syntax as

LinearModelFit

, but it returns an association containing all relevant information about the fit:

fit=N@BayesianLinearRegression[data,{1,x.},x.];

Keys[fit]

Out[]=

{LogEvidence,PriorParameters,PosteriorParameters,Posterior,Prior,Basis,IndependentVariables}

I will focus most on explaining the posterior distributions and the log-evidence. The log-evidence is a number that indicates how well the model fits with the data:

In[]:=

fit["LogEvidence"]

Out[]=

-30.0072

Is this value good or bad? We can’t tell yet: it only means something when compared to the log-evidences of other fits of the same data. I will come back to model comparison in the section on the Bayesian Occam’s razor, but let’s first take a look in more detail at the linear fit we just computed. The fit association has a key called "Posterior" that contains a number of useful probability distributions:

In[]:=

KeyValueMap[List,Short/@fit["Posterior"]]//TableForm

Out[]//TableForm=

RegressionCoefficientDistribution	MultivariateTDistribution[{-0.248878,0.714688},{{21,21},1},20.01]
ErrorDistribution	InverseGammaDistribution[10.005,3.09154]
PredictiveDistribution	StudentTDistribution-0.248878+0.714688x.,0.555877 1 ,20.01
UnderlyingValueDistribution	StudentTDistribution-0.248878+0.714688x.,0.555877 1 ,20.01

In Bayesian inference, the word “posterior” refers to a state of knowledge after having observed the data (as opposed to “prior,” which refers to the state where you do not know the data yet). The posterior distribution of the regression parameters tells us how well the parameters

and

are constrained by the data. It’s best visualized with a

ContourPlot

In[]:=

With[{coefficientDist=fit["Posterior","RegressionCoefficientDistribution"]},ContourPlot[Evaluate[PDF[coefficientDist,{a.,b.}]],{a.,-1,1},{b.,0,1},PlotRange{0,All},PlotPoints20,FrameLabel{"a","b"},ImageSize400,PlotLabel"Posterior PDF of regression coefficients"]]

Out[]=

The slant of the ellipses in the contour plot shows that there seems to be some positive correlation between

and

. You can calculate exactly how much correlation there is between them with

Correlation

In[]:=

Correlation[fit["Posterior","RegressionCoefficientDistribution"]]//MatrixForm

Out[]//MatrixForm=



1.	0.202556
0.202556	1.



At first glance, this may seem strange. Why are

and

correlated with each other? One way to look at this is to consider how the fit changes when you force one of the two coefficients to change. For example, you can fix

and then try and find the best value of

that fits the data with

FindFit

. The following animation illustrates this thought experiment:

In[]:=

With{dat=data,p=plot},AnimateWith{fit=a+b.x./.FindFit[dat,a+b.x.,b.,x.]},ShowPlotfit,{x.,-2.5,2.5},PlotLabely.==fit,PlotRange{-2.5,2.5},PlotStyleDirective[Black,Dashed],p,{{a,0.},-1,1},TrackedSymbols{a},AnimationDirectionForwardBackward

Out[]=

Intuitively, there is correlation between

and

because the center of mass of the points is slightly to the bottom left of the origin:

In[]:=

Mean[data]

Out[]=

{-0.259375,-0.434375}

This means that when

increases,

also needs to increase so that it can still catch the center of mass of the cloud. Here I’ve moved all points even further down and left to amplify the effect:

Out[]=

Another way to think of the posterior distribution over

and

is to think of the fit as consisting of infinitely many lines, drawn at random across the page and with each line weighted according to how well it fits the data. Calculating that weight is essentially what Bayesian inference is about. In the following plot, I used

RandomVariate

to draw many lines from the posterior. By making the opacity of each individual line quite low, you can see the regions where the lines tend to concentrate:

In[]:=

lines=With[{x1=0,x2=1},Apply[InfiniteLine[{{x1,#1},{x2,#2}}]&,RandomVariate[fit["Posterior","RegressionCoefficientDistribution"],2000].{{1,1},{x1,x2}},{1}]];Show[plot,Graphics[{Black,Opacity[40/Length[lines]],Thickness[0.002],lines}],PlotRange{{-2.5,2.5},All}]

Out[]=

The distribution of this cloud of lines is what I like to call the “underlying value distribution,” which is one of the distributions returned by the fit. To plot it, it’s convenient to use

InverseCDF

to calculate quantiles of the distribution. In the next example, I plotted the 5-, 50- and 95-percent quantiles, meaning that you’d expect 90% of the lines to fall within the shaded areas:

In[]:=

With[{valueDist=fit["Posterior","UnderlyingValueDistribution"],bands=Quantity[{95,50,5},"Percent"]},Show[Plot[Evaluate@InverseCDF[valueDist,bands],{x.,-5,5},Filling{1{2},3{2}},PlotLegendsbands,Prolog{Black,Opacity[15/Length[lines]],Thickness[0.002],lines}],plot,PlotRangeAll,PlotLabel"Posterior distribution of underlying values"]]

Out[]=

	95 %
	50 %
	5 %

BayesianLinearRegression also estimates the standard deviation

of the error term

, and just like the regression coefficients

and

, it has its own distribution. Variance

follows an

InverseGammaDistribution

, so to get the distribution of

I use

TransformedDistribution

In[]:=

Quiet@PlotEvaluate@PDFTransformedDistribution

,v.fit["Posterior","ErrorDistribution"],x.,{x.,0,1.5},PlotRange{0,All},FillingAxis,FrameTrue,FrameLabel{σ,"PDF"},PlotLabel

fit["Posterior","ErrorDistribution"]

Out[]=

So in short, the posterior error terms are distributed as

ϵNormalDistribution[0,σ]

, where



InverseGammaDistribution[10.005,3.09154]. This distribution we can calculate with

ParameterMixtureDistribution

In[]:=

epsDist=ParameterMixtureDistributionNormalDistribution0,

,v.InverseGammaDistribution[10.005`,3.091536049663807`];Plot[PDF[epsDist,x.],{x.,-2,2},FillingAxis,PlotRange{0,All},FrameTrue,PlotLabelϵepsDist]

Out[]=

Uncertain Uncertainties Just to make the point clear: there is uncertainty about the uncertainty of the model, which for me was an idea that was a bit difficult to get my head around initially. To make things even more complicated, all of the lines I mentioned previously also have distributions of error bars that correlate with

and

. So if you want to make a prediction from this model, you need to consider an infinite number of trend lines that each carry an infinite number of error bars. That’s a lot of uncertainty to consider!

Thinking about regression problems this way makes it clear why Bayesian inference can be a daunting task that involves lots of complicated integrals. For linear regression, though, we’re fortunate enough that it’s possible to do all of these integrals symbolically and plow our way through the infinities. This is one of the reasons why scientists and engineers like linear mathematics so much: everything remains tractable, with nice-looking formulas.

The distribution that tells you where to expect to find future data (taking all of the aforementioned uncertainties into consideration) is called the posterior predictive distribution. For the data I just fitted, it looks like this:

In[]:=

fit["Posterior","PredictiveDistribution"]

Out[]=

StudentTDistribution-0.248878+0.714688x.,0.555877

1.05211+0.0164952x.+0.031814

,20.01

Now I will plot this distribution using the function regressionPlot1D from the Git repository, which is just shorthand for the underlying values plot I showed earlier. I also included the distribution you get when you make a “point estimate” of the predictive model. This means that you take the best values of

and

from the posterior and plot

ya+bx+ϵ

using those values as if they were completely certain. Point estimates can be useful in making your results more manageable since probability distributions can be difficult to work with, but reducing your distributions to single numbers always discards information. Here I used the

Mean

(expected value) to reduce the distributions to single values, but you could also use other measures of centrality like the median or mode:

In[]:=

WithpredictiveDist=fit["Posterior","PredictiveDistribution"],pointEst=NormalDistributionMean@fit["Posterior","RegressionCoefficientDistribution"].{1,x.},MeanTransformedDistribution

,v.fit["Posterior","ErrorDistribution"],bands={0.95,0.5,0.05},Show[regressionPlot1D[predictiveDist,{x.,-5,5}],regressionPlot1D[pointEst,{x.,-5,5},PlotStyleDashed,PlotLegends{"Point estimate"}],plot,PlotRangeAll,PlotLabel"Posterior predictive distribution"]

Out[]=

	95 %
	50 %
	5 %

Point estimate

So that covers fitting a straight line through the points. However, the shape of the data suggests that maybe a quadratic fit is more appropriate (I cooked up the example just for this purpose, of course), so this is a good moment to demonstrate how Bayesian model comparison works. I am going to fit polynomials up to degree 5 to the data, with each of these polynomials representing a model that could explain the data. From this micro-universe of possible explanations of the data, I want to select the best one.

The Bayesian Occam's Razor In Bayesian inference, the models have probability distributions in the same way that the regression coefficients

and

have distributions. The posterior probability of a model depends on two factors:

◼

A quantity commonly referred to as evidence or marginal likelihood. This measures how well a model fits the data while taking into account the uncertainties in the regression parameters such as

and

. BayesianLinearRegression reports this quantity as the "LogEvidence": the higher the log-evidence, the better the model. The evidence automatically accounts for model complexity due to an effect that is sometimes referred to as the Bayesian Occam’s razor (see, for example, chapter 28 in this book by David MacKay).

◼

The prior probability of the model.

Usually you’ll only want to consider models that, before seeing the data, you consider approximately equally likely. However, sometimes you do have strong prior information available that the data was generated by a certain model, and you’d only accept other explanations if the evidence overwhelmingly disproves that model in favor of another.

For example, if the data shown previously came from a noisy (and possibly biased) measurement of current through a standard resistor as a function of voltage, I would have little doubt that Ohm’s law applies here and that fitting a straight line is the correct way to interpret the data. It’s simply more likely that a little noise accidentally made the data look quadratic rather than that Ohm's law suddenly stopped working in my experiment. In this case, I would assign a prior probability very close to 1 to Ohm’s law and divide the remaining sliver of probability among competing models I’m willing to consider. This is in accordance with a principle called Cromwell’s rule, which was coined by Dennis Lindley:

“Leave a little probability for the Moon being made of green cheese; it can be as small as 1 in a million, but have it there since otherwise an army of astronauts returning with samples of the said cheese will leave you unmoved.”

So let’s see how the fits of different polynomials look and what their evidence is:

In[]:=

models=AssociationMapBayesianLinearRegressiondata,

Range[0,#]

,x.&,Range[0,5];Multicolumn[KeyValueMap[Show[regressionPlot1D[#2["Posterior","PredictiveDistribution"],{x.,-5,5}],plot,PlotRangeAll,PlotLabelStringForm["Degree: `1`\nLog evidence: `2`",#1,#2["LogEvidence"]]]&,models],2]

Out[]=

	95 %
	50 %
	5 %

	95 %
	50 %
	5 %

	95 %
	50 %
	5 %

	95 %
	50 %
	5 %

	95 %
	50 %
	5 %

	95 %
	50 %
	5 %

As you can see, the first- and second-degree fits are in close competition for being the best model. If we assume that the prior probabilities for all models are the same, we can calculate the relative probabilities of the fits by exponentiating the log-evidences and normalizing them. For this operation, we can borrow

SoftmaxLayer

from the neural network framework:

calculateRelativeProbabilities[models_]:=Module[{probabilities,temp=SortBy[-#LogEvidence&]@models},probabilities=SoftmaxLayer[]@Replace[temp[[All,"LogEvidence"]],assoc_?AssociationQValues[assoc]];temp[[All,"Probability"]]=probabilities;temp[[All,"AccumulatedProbability"]]=Accumulate[probabilities];temp];

models=calculateRelativeProbabilities[models];Dataset[models][All,{"LogEvidence","Probability","AccumulatedProbability"}]

Out[]=

	LogEvidence	Probability	AccumulatedProbability
1	-30.0072	0.538841	0.538841
2	-30.1774	0.454513	0.993354
3	-34.4292	0.00647144	0.999825
4	-38.7037	0.0000900772	0.999915
0	-38.787	0.0000828723	0.999998
5	-42.6154	0.00000180207	1.0

In this table, I sorted the models by probabilities and then accumulated them so you can easily see how much total probability you’ve covered going from top to bottom. Keep in mind that these probabilities only mean anything in the micro-universe we’re currently looking at: the moment I would start adding more models, all of these probabilities could change. In Bayesian inference, you can only compare models you’re willing to formulate in the first place: there is no universal measure that tells you if the fit of the data is actually good or not in some completely objective sense. All you can really hope for is to find the best model(s) from the group you’re considering because the set of models you’re not considering is impossible to compare against.

Interlude: Comparing against LinearModelFit

At this point, it’s also good to compare the fits by BayesianLinearRegression against those by LinearModelFit, which also presents several goodness-of-fit measures:

standardFits=AssociationMapLinearModelFitdata,

Range[0,#]

,x.,ConfidenceLevel0.90(*Toagreewithourearlierchoiceof95%and5%predictionbands*)&,Range[0,5];

With[{keys={"BIC","AIC","AICc","AdjustedRSquared","RSquared"}},AssociationThread[keys,#[keys]]&/@standardFits//Dataset]

Out[]=

	BIC	AIC	AICc	AdjustedRSquared	RSquared
0	64.8858	62.8943	63.6002	0.0	0.0
1	42.3227	39.3355	40.8355	0.707106	0.722521
2