# This demonstration gives a slightly different view of the two stage sampling
# process demonstrated in gp_minimal.py. See that file for more details. Here we
# reinforce that sampling from the posterior really is just continuing the prior
# sampling process given the values that we've seen. I like this version of the
# demo because there's less linear algebra than in gp_minimal.py, but this
# presentation is less standard.
#
# Iain Murray, November 2016
import numpy as np
import matplotlib.pyplot as plt
## The kernel function (as in gp_minimal.py)
######################################################################
rbf_fn = lambda X1, X2: \
np.exp((np.dot(X1,(2*X2.T))-np.sum(X1*X1,1)[:,None]) - np.sum(X2*X2,1)[None,:])
gauss_kernel_fn = lambda X1, X2, ell, sigma_f: \
sigma_f**2 * rbf_fn(X1/(np.sqrt(2)*ell), X2/(np.sqrt(2)*ell))
k_fn = lambda X1, X2: gauss_kernel_fn(X1, X2, 3.0, 10.0)
## Sampling from the prior
######################################################################
# Pick the input locations that we want to see the function at.
X_train = np.array([2,4,6,8])[:,None] + 0.01
X_test = np.arange(0, 10, 0.02)[:,None]
X_all = np.vstack([X_train, X_test])
N_train = X_train.shape[0]
N_all = X_all.shape[0]
# The joint distribution over function values has zero mean and covariance
# K_all = np.dot(L_all, L_all.T)
K_all = k_fn(X_all, X_all) + 1e-9*np.eye(N_all)
L_all = np.linalg.cholesky(K_all)
# Function values can be sampled with: L_all*nu, where nu = randn(N_all).
# Because L_all is lower-triangular, the first N_train function values depend
# only on the first N_train values of nu. We pick those first:
nu1 = np.random.randn(N_train)
plt.figure(1)
plt.clf()
for ii in range(3):
# Then we consider different samples from the prior that complete those
# first N_train values in different ways:
nu2 = np.random.randn(N_all - N_train)
nu = np.hstack([nu1, nu2])
f_all = np.dot(L_all, nu)
# These x's will fall on top of each other for each loop, as nu1 is shared:
plt.plot(X_train, f_all[:N_train], 'x', markersize=20, markeredgewidth=2)
# But we'll get different completions for different nu2. These are
# samples from the posterior given the 'x' observations.
plt.plot(X_test, f_all[N_train:], '-', linewidth=2)
plt.legend(['train points', 'completions / posterior samples'])
plt.xlabel('x')
plt.ylabel('f')
plt.show()
# Want to see samples from the posterior given noisy observations? You could
# insert the following two lines beneath the definition of K_all:
#noise_var = 1.0
#K_all[:N_train, :N_train] = K_all[:N_train, :N_train] + noise_var*np.eye(N_train)
# You could extend the demo to plot a mean and error bars like in gp_minimal.py
# Of course we don't see the random numbers nu1 directly when we observe data.
# However, they are known: we can solve for nu1 from the observed values.
# I should use specialist triangular solver, but just an illustration:
nu1_from_obs = np.linalg.solve(L_all[:N_train, :N_train], f_all[:N_train])
assert(np.max(np.abs(nu1_from_obs - nu1)) < 1e-9)
# Notice how almost all of the code above is comments, plotting, and tracking
# which data points are which. Little maths is required to sample realizations
# of complex functions given data.