Bayesian Linear Regression Predictive Distribution

Question

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

def radial_basis_null():
    return 1 # constant factor

def vector_basis_funtions(x: int):
    phi_0 = radial_basis_null()
    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[0],xs))

Answer 1

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*[email protected]@t

        return self

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

        Phi = make_design(x)
        preds = [email protected]
        cov_mat = 1/self.beta + [email protected]@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:

Re: Your questions

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).

Answer 2

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.