• Here we specify a prior over functions specified by a particular Gaussian process which favours smooth functions
• We assume a 0 mean for every x , since we don't have any prior knowledge
• We then add two data points and only consider functions passing through these two points, ending up with the posterior in 1.1 b

• is of observations
• = input vector of dimension
• = scalar output/target
• = design matrix with all input vectors as column vectors

Bayesian analysis of the standard linear regression model with Gaussian noise:

1

• Often a bias is included, but this can easily be implemented by adding a extra element to the input vector , whose value is always
• The noise assumption ( to be normally distributed) together with the model direcly gives us the likelihood, ergo the probability density of the observations given the parameters.

Because of (assumed) independence between samples:

2

In the Bayesian formalism we need to specify a prior over the parameters, expressing our beliefs about the parameters before considering the observations. We simply give a zero mean Gaussian with covariance matrix on the weights

3

Inference in a Bayesian linear model is based on the posterior distribution over the weights:

4

where the normalization constant or marginal likelihood is given by

5

If we can obtain without explicitly calculating this integral, this is what they call "integrating out the weights".

For now we ignore the marginal distribution in the denominator, and then we have

6

where . Notice how this has the form of a Gaussian posterior with mean and covariance matrix :

7

where .

The meaning of the log prior and the MAP estimate substantially differs between Bayesian and non-Bayesian settings:

• Bayesian MAP estimate plays no special role
• Non-Bayesian negative prior is sometimes used as penalty term, and MAP point is known as the penalized maxium likelihood estimate of the weights. In fact, this is what is known as ridge regression, because of the quadratic term from the log prior.

To make predictions for some input , we average over all the parameter values, weighted by their posterior probability. Thus the predictive distribution for , is given by averaging (or rather, taking the expectation for over the weights ) the output of all possible linear models w.r.t. the Gaussian posterior

8

And even the predictive distribution is Gaussian! With a mean given by the posterior mean of the weights from eq. (2.8) multiplied by the new input. All hail Gauss!

Also note how the predictive variance is a quadratic of the new input, and so the prediction will be more uncertain as the magnitude of test input grows.

• Can overcome linearity by projecting features into higher-dimensional space allowing us to implement polynomial regression
• If projections are independent of parameters => model is still linear in the parameters and therefore analytically solvable/tractable!
• Data which is not linearly separable in the original data space may become linearly separable in a high dimensional feature space.

GP is defined as collection of rvs a consistency requirement. This simply means that if the GP e.g. specifies ,

The defintion does not exclude GPs with finite index sets (which would just be Gaussian distributions).

Example of GP can be obtained from our Bayesian linear regression model with prior .

11

Ergo, and are jointly Gaussian with zero mean and covariance given by .

• We know
• = vector of training outputs
• = vector of test outputs
• Join distribution of training and test outputs according to the prior is then

12

• If training points and test points, the matrix = covariaces evaluated at all pairs of training and test points, and similarily for , etc.
• Restrict this joint prior distribution to contain only those functions which agree with the observed data points
  • A lot of functions mate not efficient

  • In probabilistic terms this operation is simple: equivalent of conditioning on the joint Gaussian prior distribution on the observations to give

    13
• Now we can sample values for (corresponding to test inputs ) from the joint posterior distribution by evaluating the mean and covariance matrix from the eqn. above.

• Both forward and backward substitution requires operations, when is of size .
• The computation of the Cholesky factor is considered numerically extremely stable and takes time .

• = input space
• = feature space
• 
• = matrix with columns for all training inputs

17

where is of length . This is analogous to the standard linear model except everywhere is substituted for .

First we import the stuff that needs to be imported

    %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

    from __future__ import division

After having some issues, I found this article from which I used some of the code, so thanks a lot to Chris Fonnesbeck!

We start out by the following specification for our Gaussian Process:

18

def exponential_cov(x, y, theta):
    return theta[0] * np.exp(-0.5 * theta[1] * np.subtract.outer(x, y) ** 2)

Then to compute the conditional probability for some new input conditioned on the previously send inputs with outputs, we use the following:

19

where, again, remember that and represents the test and training outputs, i.e in the case of training (already seen) inputs:

20

def conditional(X_star, X, y, theta):
    B = exponential_cov(X_star, X, theta)
    C = exponential_cov(X, X, theta)
    A = exponential_cov(X_star, X_star, theta)

    mu = np.linalg.inv(C).dot(B.T).T.dot(y)
    sigma = A - B.dot(np.linalg.inv(C).dot(B.T))

    return mu.squeeze(), sigma.squeeze()  # simply unwraps any 1-D arrays, so (1, 100) => (100,)

For now we set :

theta = (1, 1)
sigma_theta = exponential_cov(0, 0, theta)
xpts = np.arange(-3, 3, step=0.01)
plt.errorbar(xpts, np.zeros(len(xpts)), yerr=sigma_theta, capsize=0)
    plt.ylim(-3, 3)

    # some sampled priors
    f1_prior = np.random.multivariate_normal(np.zeros(len(xpts)), exponential_cov(xpts, xpts, theta))
f2_prior = np.random.multivariate_normal(np.zeros(len(xpts)), exponential_cov(xpts, xpts, theta))
f3_prior = np.random.multivariate_normal(np.zeros(len(xpts)), exponential_cov(xpts, xpts, theta))
plt.plot(xpts, f1_prior)
plt.plot(xpts, f2_prior)
plt.plot(xpts, f3_prior)

We select an arbitrary starting point, say , and simply sample from a normal distribution with unit covariance. This is just for us to have a "seen" sample to play with.

x = [1.]
y = [np.random.normal(scale=sigma_theta)]
np.array(y)

We now update our confidence band, given the point we just sampled, by using the covariance function to generate new point-wise confidence intervals, conditional on the value seen above.

sigma_1 = exponential_cov(x, x, theta)  # compute covariance for new point

def predict_noiseless(x, data, kernel, theta, sigma, f):
    K = [kernel(x, y, theta) for y in data]           # covariance between test input and training input
    sigma_inv = np.linalg.inv(sigma)                      # inverting the covariance matrix for training input
    K_dot_sigma_inv = np.dot(K, sigma_inv)
    y_pred = K_dot_sigma_inv.dot(f)                   # compute the predicted mean 
    sigma_new = kernel(x, x, theta) - K_dot_sigma_inv.dot(K)  # covariance of new point

    return y_pred, sigma_new

x_pred = np.linspace(-3, 3, 2000)
predictions = [predict_noiseless(i, x, exponential_cov, theta, sigma_1, y) for i in x_pred]

y_pred, sigmas = np.transpose(predictions)
plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0, alpha=0.6)
plt.plot(x, y, "ro")
plt.plot(x_pred, y_pred, c="c")
plt.ylim(-3, 3)

Notice that here we are predicting a mean and covariance at each test point, but for fun's sake we can sample some functions too!

21

That is; we create a normal distribution of dimension equal to the number of test inputs, draw from this distribution, then the sample at index i corresponds to , yah see? It's weird, yes, but allows the indexed sample points to be any real value!

    # same as above
plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0, alpha=0.3, label="$\sigma_*$")
plt.plot(x, y, "ro")
plt.plot(x_pred, y_pred, c="black", label="$\mu_*$")
plt.ylim(-3, 3)
    plt.xlim(-3, 3)

# Some sampled functions
mu_star, K_star = conditional(x_pred, x, y, theta)       # x_pred conditioned on training input (x, y)
f1 = np.random.multivariate_normal(mu_star, K_star)
f2 = np.random.multivariate_normal(mu_star, K_star)
f3 = np.random.multivariate_normal(mu_star, K_star)
plt.plot(x_pred, f1, label="$f_*^{(1)}$")
plt.plot(x_pred, f2, label="$f_*^{(2)}$")
plt.plot(x_pred, f3, label="$f_*^{(3)}$")
    plt.legend()