I’m happy to announce the release of BayesianLinearRegression.jl. Initial work on this package began in 2015, but there was some significant bitrot. Thanks to RelationalAI for supporting this development!

I’ve started registration for this package, but it’s currently in the waiting period. It should be registered within the next few days.

If you’re new to Bayesian modeling, you can think of this as ridge regression with automatically-tuned regularization.

For all the details, see Chapter 3 of Pattern Recognition and Machine Learning by Chris Bishop.

Say we have something like

```
julia> X = randn(100,10);
julia> β = randn(5);
julia> y = X[:,1:5] * β + randn(100);
```

Then we can say

```
julia> m = BayesianLinReg(X,y);
julia> fit!(m)
BayesianLinReg model
Log evidence: -156.14
Prior scale: 0.57
Noise scale: 0.97
Coefficients:
1: 1.54 ± 0.09
2: 0.46 ± 0.10
3: 0.55 ± 0.09
4: 0.16 ± 0.11
5: -0.43 ± 0.10
6: 0.03 ± 0.10
7: 0.03 ± 0.10
8: 0.18 ± 0.11
9: -0.07 ± 0.11
10: 0.04 ± 0.10
```

This shows a few things:

- The
*log evidence*, which can be used for model selection, - The
*prior scale*, a “reasonable” prior for each weight, - The
*noise scale*, which you can think of as \hat{\sigma}, and - The
*coefficients*, or posterior weights, with standard errors.

Note that thanks to the clever tricks of `Measurements.jl`

, the posterior weights include correlation information:

```
julia> w = posteriorWeights(m);
julia> w[1].der
Measurements.Derivatives{Float64} with 10 entries:
(0.0, 1.0, 0x000000000000000a) => -0.00411411
(0.0, 1.0, 0x0000000000000008) => 0.0106236
(0.0, 1.0, 0x0000000000000006) => -0.00479298
(0.0, 1.0, 0x0000000000000004) => -0.00320726
(0.0, 1.0, 0x0000000000000002) => 0.0920992
(0.0, 1.0, 0x0000000000000003) => 0.00134676
(0.0, 1.0, 0x0000000000000005) => -0.00552799
(0.0, 1.0, 0x0000000000000007) => -0.00393502
(0.0, 1.0, 0x0000000000000009) => 0.000518892
(0.0, 1.0, 0x000000000000000b) => -0.00409676
```

We can also change the *active weights* (any inactive are set to zero):

```
julia> m.active = [1:5;];
julia> update!(m)
BayesianLinReg model
Log evidence: -149.04
Prior scale: 1.00
Noise scale: 1.00
Coefficients:
1: 1.56 ± 0.10
2: 0.48 ± 0.10
3: 0.54 ± 0.09
4: 0.16 ± 0.10
5: -0.42 ± 0.11
julia> fit!(m)
BayesianLinReg model
Log evidence: -148.69
Prior scale: 0.80
Noise scale: 0.96
Coefficients:
1: 1.56 ± 0.09
2: 0.48 ± 0.10
3: 0.54 ± 0.08
4: 0.16 ± 0.10
5: -0.42 ± 0.10
```

So you could, for example, build some variable selection methods around this.

For prediction, we have

```
julia> predict(m,X[1:10,m.active])
10-element Array{Measurements.Measurement{Float64},1}:
-0.62 ± 0.99
-2.18 ± 0.97
1.21 ± 0.99
0.33 ± 0.99
-2.58 ± 0.99
-1.94 ± 0.98
1.4 ± 1.0
-0.05 ± 0.97
-3.51 ± 0.99
-0.97 ± 0.98
```

You can also do this without observation noise:

```
julia> predict(m,X[1:10,m.active]; noise=false)
10-element Array{Measurements.Measurement{Float64},1}:
-0.62 ± 0.22
-2.18 ± 0.14
1.21 ± 0.24
0.33 ± 0.24
-2.58 ± 0.23
-1.94 ± 0.2
1.36 ± 0.29
-0.05 ± 0.14
-3.51 ± 0.24
-0.97 ± 0.17
```

Or ignore uncertainty altogether:

```
julia> predict(m,X[1:10,m.active]; uncertainty=false)
10-element Array{Float64,1}:
-0.6198367846446444
-2.1822748489031185
1.2097200082593464
0.3324144009375185
-2.5818035785619786
-1.9366523399339064
1.3611242250751907
-0.0458769531369336
-3.5130276372775886
-0.9710314161681255
```

Finally, when you `fit!`

, there’s an option to pass a `callback`

. This is expected to be a function taking a `BayesianLinReg`

as its argument. One useful thing a callback can do is set `m.done = true`

if some convergence criteria are met. The default is `callback = stopAtIteration(10)`

, which does what it says on the tin