# Bayesian Linear Regression Predictive Distribution

This is a homework assignment. I need a hint in the right direction and a couple explanations as things are not very clear to me. What is the predictive distribution? How do I draw from it? The notation is quite confusing for me as I'm missing the big picture right now.

Task: Generate new data from $$f(x)=sin(2\pi x)$$ where $$x\in[0,1]$$ and add Gaussian noise with a standard deviation of $$\beta^{-1/2}=0.3$$ to the data points generated. Explore data sets of various size. Consider a model of consisting of 1 constant function $$\phi_0$$ and 9 Gaussian radial basis functions with identical width $$s$$ and means $$\mu_j$$ equally distributed between $$0$$ and $$1$$. Explore different widths of for the basis functions. Plot an analog of Bishop page 157 figure 1 by computing the predictive distribution given the following information:

Consider a target variable $$t$$ given by a deterministic function

$$y(\mathbf{x},\mathbf{w})=\mathbf{w}^T\phi(\mathbf{x})$$

depending on input $$\mathbf{x}$$ and parameters $$\mathbf{w}$$ with additive gaussian noise so that $$t=y(\mathbf{x}, \mathbf{w}) + \epsilon$$, where $$\epsilon$$ is a zero mean Gaussian random variable with precision parameter $$\beta$$.

This can also be written as

$$p(t|\mathbf{x},\mathbf{w}, \beta)=\mathcal{N}(t|y(\mathbf{x},\mathbf{w}),\beta^{-1}).$$

Choosing a Gaussian prior

$$p(\mathbf{w}|\alpha)=\mathcal{N}(\mathbf{w}|\mathbf{0},\alpha^{-1}\mathbf{I})$$

the predictive distribution is also Gaussian and given by

$$p(t|\mathbf{x},\mathbf{t},\mathbf{X},\alpha,\beta)=\mathcal{N}(t|\mathbf{m}_N^T\phi(\mathbf{x}),\sigma_N^2(\mathbf{x}))$$

with mean

$$\mathbf{m}_N=\beta\mathbf{S}_N\mathbf{\Phi}^T\mathbf{t}$$

and variance

$$\sigma_N^2(\mathbf{x})=\frac{1}{\beta}+\phi(\mathbf{x})^T\mathbf{S}_N\phi(\mathbf{x}),$$

where the Matrix $$\mathbf{S}_N$$ is defined as

$$\mathbf{S}_N^{-1}=\alpha\mathbf{I}+\beta\mathbf{\Phi}^T\mathbf{\Phi},$$

the vector of basis functions is given as

$$\phi(\mathbf{x}_n)=(\phi_0(\mathbf{x}_n),\phi_1(\mathbf{x}_n),\ldots,\phi_{M-1}(\mathbf{x}_n))^T,$$

the matrix of basis functions given by

$$\mathbf{\Phi}=\left( \begin{array}{rrrr} \phi_0(\mathbf{x}_1) & \phi_1(\mathbf{x}_1) & \cdots & \phi_{M-1}(\mathbf{x}_1) \\ \phi_0(\mathbf{x}_2) & \phi_1(\mathbf{x}_2) & \cdots & \phi_{M-1}(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_0(\mathbf{x}_N) & \phi_1(\mathbf{x}_N) & \cdots & \phi_{M-1}(\mathbf{x}_N) \\ \end{array}\right)$$

with the vectors of input training data $$\mathbf{X}=\{\mathbf{x_1}, \mathbf{x_2},\ldots,\mathbf{x_N}\}$$ and corresponding output training values $$\mathbf{t}=\{t_1,\ldots,t_n\}$$ and the value t to be predicted for a new input $$\mathbf{x}$$.

Code: This is what I've got so far. As asked, I'm drawing random samples from the sin function and adding noise. Also, I've already defined the radial function. What should I do next?

from math import sin, pi, exp, e
from random import random, gauss
import matplotlib.pyplot as plt
import numpy as np
from typing import List

AMOUNT_RBA = 10 # amount of radial basis functions
ALPHA = 1 # what is this?
BETA = 1 # what is this?

def gen_data_point() -> int:
rand = random()
noise = gauss(0, 0.3**2)
return rand, sin(2 * pi * rand) + noise

def radial_basis(j: int, x: int, amount=9) -> int:
s = 1  # parameter to play with
mu_j = (j-1)/(amount-1)  # spreads the mean evenly over [0,1]
numerator = (x-mu_j)**2
denumerator = 2*(s**2)
phi = e**(-numerator/denumerator)
return phi

return 1 # constant factor

def vector_basis_funtions(x: int):
phis = [radial_basis(i, x) for i in range(1, AMOUNT_RBA)]
basis_func_vector = np.matrix([phi_0] + phis)
return basis_func_vector.T
# print(np.concatenate((basis_func_vector, basis_func_vector)).T)

def matrix_basis_functions(xs: int):
# initalize matrix with the needed dimensions
n, m = len(xs), AMOUNT_RBA
Phi = np.zeros([n, m])
# add row to matrix for every x
for i, x in enumerate(xs):
Phi[i] = vector_basis_funtions(x).T
return Phi

def S_N(xs: List[int]):
identity_matrix = np.identity(AMOUNT_RBA)
Phi = matrix_basis_functions(xs)
S_N = ALPHA*identity_matrix + BETA*Phi.T*Phi
return np.linalg.inv(S_N)

def var_N(x: int, xs: List[int]):
precision = 1/BETA
var = vector_basis_funtions(x).T * S_N(xs) * vector_basis_funtions(x)
return precision + var

def mean_N():
pass

# print(S_N(list(range(10))))
xs = list(range(10))
print(var_N(xs,xs))


Ok, so this was fun. Luckily I have the book on hand and so if you follow some equations of Section 3.3.1 then (I think) it becomes quite easy. Let's write a class to do all the stuff we need.


import numpy as np
import matplotlib.pyplot as plt

def make_design(x, s=0.1, p = 9):

design = np.zeros(shape = (x.size, p))

design[:,0] = 1

for i in np.arange(1,p):
design[:,i] = np.exp(-(x-i/10)**2/(2*s))

return design

class BayesianLinearModel():

def __init__(self, alpha=.010, beta=(1/.3)**2):

self.alpha = alpha
self.beta = beta

def fit(self,x,t):

Phi = make_design(x)
r,c = Phi.shape
self.SN = np.linalg.pinv(self.alpha*np.eye(c) + self.beta*Phi.T@Phi)
self.mN = self.beta*self.SN@Phi.T@t

return self

def predict(self, x, return_predictive = True):

Phi = make_design(x)
preds = Phi@self.mN
cov_mat = 1/self.beta + Phi@self.SN@Phi.T
if return_predictive:
sig = np.sqrt(np.diag(cov_mat))
return preds , sig
else:
return preds


Everything is quite easy now.

true_func = lambda x: np.sin(2*np.pi*x)
x = np.random.rand(25)

y = true_func(x) + np.random.normal(0, 0.3, size = x.size)
X = np.linspace(0,1,1001)

model = BayesianLinearModel().fit(x,y)
preds,lims = model.predict(X, return_predictive = True)
fig, ax = plt.subplots(dpi=120)

plt.fill_between(X, preds - 2*lims, preds + 2*lims, color = 'red', alpha = 0.25)
plt.plot(X, preds, color ='red', label = 'prediction')
plt.plot(X, true_func(X) , color = 'green', label = 'truth')
plt.scatter(x,y, marker =  'o', alpha = 0.5, edgecolor = 'k', facecolor ='None')

plt.legend()


No guarantees this is correct, but it certainly looks correct. EDIT:

when fitting the model we used training data to build the design matrix and to build SN, mN. When predicting we use the new data to build another design matrix and use it to predict. How/when do I know when to use which x?

You only use the training x (in my code x) once, and you use it in the fit call. You only have to fit once. When you want to predict on new data (in my code, X) then you pass that to the predict method.

Could elaborate on the drawing from the pred. distr.? From what I understand in ur code, we only use the computed mean to draw and only use the variance to obtain the uncertainty? I thought I would have to plug those values into a normal distr.

You are technically right. We can draw from the predictive distribution which induces a distribution over curves. I've not implemented that here because I'm lazy. The reason I've made the predictive distribution as I have is because our choice of prior and likelihood make it so our posterior is Gaussian. I can thus summarize the predictive distribution as mean +/- two standard deviations (which is exactly what I've shown in the shaded region).

• Thank you! Your code is quite concise, I still need to learn what's out there in the numpy world.. Question, why is your design matrix of dimension x.size x 25? Shouldn't it be x.size x basis_func_amount(=10)? Jan 7, 2020 at 15:06
• @Marc whoops, that was from me fooling around with the implementation. You can easily change the number of basis functions with the p argument in the design matrix function Jan 7, 2020 at 15:23
• Ah ok. I figured, wasnt sure though. Couple questions: 1) when fitting the model we used training data to build the design matrix and to build SN, mN. When predicting we use the new data to build another design matrix and use it to predict. How/when do I know when to use which x? 2) Could elaborate on the drawing from the pred. distr.? From what I understand in ur code, we only use the computed mean to draw and only use the variance to obtain the uncertainty? I thought I would have to plug those values into a normal distr. (as the definition suggests), but not sure what would happen next. Jan 7, 2020 at 15:30
• Thanks Demetri, you helped a lot! Jan 8, 2020 at 14:20

I will be working backwards from the plot you are supposed to generate.

The plot on page 157 of Bishops book on pattern recognition is a plot with the values of the input $$x$$ on the 1st axis and $$t$$ on the 2nd axis.

The model for the data generating proces is

$$t = sin(2\pi x) + \epsilon,$$

with $$x \in [0,1]$$ and $$\epsilon \sim \mathcal N(0,\sigma^2)$$ with $$\sigma^2 = 1/\beta$$ hence $$\beta$$ being the precision and $$\sigma^2$$ the variance. The plot itself contains the true model so you need a function $$f(x) = sin(2\pi x)$$ and a grid of point in the interval $$[0,1]$$ to plot the true function. This can be done in R in the following manner

G <- 100
f <- function(x)
{
out <- sin(2*pi*x)
return(out)
}

x_grid <- seq(0,1,length.out=G)
t_star <- f(x_grid)
plot(x_grid,t_star,type="l",col="green",ylim=c(-1.5,1.5),ylab=t,xlab=x)


Secondly you need to plot the mean of the predictive distribution which is a normal distribution with some mean $$\mu(x)$$ and some variance $$\sigma_N^2(x)$$ both being functions of the $$x$$ hence for each value of $$x$$ there is a predictive distribution $$p(t\lvert x) = \mathcal N(\mu(x),\sigma_N^2(x))$$. Once you have the functions $$\mu(x)$$ and $$\sigma^2(x)$$ it is easy to plot $$\mu(x) = \mathbf m_N^\top \phi(x)$$ for different $$x$$ to get the read line in the Bishop plot page 157 and the light read area is probably a 95% credibility interval so find 0.025 and 0.975 quantiles of the distribution predictive distribution $$p(t\lvert x) = \mathcal N(\mu(x),\sigma_N^2(x))$$ for different values of x and plot them.

So changing the workflow no longer working backwards here are the steps you need to take:

(1) Generate $$N$$ observations of the input variable $$x$$ on interval $$[0,1]$$ can be done drawing random uniform $$\mathcal U(0,1)$$ and store in vector $$\mathbf X = (x_1,...,x_N)^\top$$.

(3) Set values for $$\beta$$ and $$\alpha$$.

(4) Generate $$\mathbf t = (t_1,...,t_N)^\top$$ by the formula $$t_i = sin(2\pi x_i) + \epsilon_i$$ where $$\epsilon_i \sim \mathcal N(0,1/\beta)$$.

You have now succesfully generated a dataset $$\{t_i,x_i\}_{i=1}^N$$ according to the datagenerating process described.

Next step is to find the predictive distribution which you know is normal so you need to find mean and variance. Since however you need the mean and variance for different values of $$x$$ the mean and variance are functions of $$x$$. So you need to wirte a function that for a given x computes the mean and variance

(4) Calculate the $$\mathbf m_N$$ this can be done by solving $$\alpha \mathbf I + \beta \Phi^\top\Phi = \beta \Phi^\top \mathbf t$$

(5) Calculate $$S_N^{-1}$$ to be used in the variance formula $$\sigma^2_N(x)$$ for different values of $$x$$

.... to be continued if you have more questions.

• I haven't read your answer yet as I wanna try a little more by myself, but thanks so far. I updated my code above, I hope you can read Python and NumPy. What I have done so far is implement all the functions (or matrices, which kind of act like functions) and see if the dimensions that come out match. Does this look like the right approach? I'm definitely lacking an understanding of the terminology and what I'm supposed to do here. Jan 6, 2020 at 22:43
• The answer is schematic and incomplete not a full blown solution. I do not use python but you seems to be on the right track from how I can read the code. Jan 6, 2020 at 22:51
• Ok, I read your answer and I think I now know what to do. I think my approach is pretty much what you have described. Really, everything depends on $X$ and the single $x$, the functions/matrices are given so its just plugging in I guess. Thinking about it, how do I obtain the read area, which I think refers to the uncertainty? Jan 6, 2020 at 22:53
• Calculate quantiles of the predictive distribution as explained in the answer. Ask for further details if needed. Jan 6, 2020 at 22:55
• Ah. That makes sense. I will get back to you if I get stuck. Thank you, you helped a lot! Jan 6, 2020 at 22:56