<p>Just an update for anyone who looked at this. I found the solution reading this paper:</p>

<p><a href="https://arxiv.org/pdf/1711.10834.pdf" rel="nofollow noreferrer">https://arxiv.org/pdf/1711.10834.pdf</a></p>

<p>and the following code:</p>

<pre><code>mean_function = np.dot(np.dot(K_pt ,np.linalg.inv(K_train)), ytrain) 

covariance_function = K_prior - np.dot(np.dot(K_pt ,np.linalg.inv(K_train)) , K_tp) 

f = np.random.multivariate_normal(mean_function[:,0],covariance_function , 100)
</code></pre>

<p>where f is the posterior joint Gaussian from which you sample from</p>
Gaussian Process Posterior (Python)
<p>I have created and sampled a jointly Gaussian prior with mean=0 using the code below:</p>

<pre><code>import numpy as np
import matplotlib.pyplot as plt 
from math import pi 
from scipy.spatial.distance import cdist
import scipy.stats as sts

x_prior = np.linspace(-10,10,101)
x_prior = x_prior.reshape(-1,1)
mu = np.zeros(x_prior.shape)

#defining the Kernel for the covariance function

def sec(a,b, length_scale , sigma) : 
    K = sigma * np.exp(-1/(2*length_scale) * cdist(a,b)**2)
    return K 

#defining the Gaussian Process prior

def GP(a , b, mu , kernel , length_scale, sigma , samples ) :
    f = np.random.multivariate_normal(mu.flatten(), kernel(a ,b , length_scale , sigma ) , samples)
    return f

prior = GP(x_prior ,x_prior, mu , sec , 100, 1 , 5)

plt.figure()
plt.grid()
plt.title('samples from the Gaussian prior')
plt.plot(x_prior , prior.T)
plt.show()
</code></pre>

<p><a href="https://i.stack.imgur.com/w5Yim.png" rel="noreferrer"><img src="https://i.stack.imgur.com/w5Yim.png" alt="GP prior sampling"></a></p>

<p>Then, when adding in some 'observed' data, I wish to compute the posterior over these points but this is where I become stuck. </p>

<p>Here's my code for introducing new data:</p>

<pre><code>x_train = np.array([-10,-8,5,-1,2])
x_train = x_train.reshape(-1,1)
def straight_line(m , x , c):
    y = 5*x + c
    return y
ytrain = straight_line(5 , x_train , 0)
</code></pre>

<p>It's my understanding that you calculate a conditional distribution over the new data given the prior and new x values associated with the observed data.</p>

<p>Do you then wish to update the multivariate prior to become the posterior by performing some sort of change to the mean values to include the new y values? </p>

<p>I have used the following resources to try and attempt this:</p>

<p><a href="http://katbailey.github.io/post/gaussian-processes-for-dummies/" rel="noreferrer">http://katbailey.github.io/post/gaussian-processes-for-dummies/</a>
<a href="https://www.robots.ox.ac.uk/~mebden/reports/GPtutorial.pdf" rel="noreferrer">https://www.robots.ox.ac.uk/~mebden/reports/GPtutorial.pdf</a></p>

<p>but I'm really trying to understand what happens at each stage, and why, so that when I get a posterior (which I can't do) I know exactly how I got there. </p>

<p>Here's some solutions I've been trying to implement but so far no avail:</p>

<pre><code>K_train = sec(x_train , x_train , 1,1)
K_prior = sec(x_prior , x_prior , 1,1)
K_pt =  sec(x_prior , x_train , 1,1)
K_tp = sec(x_train , x_prior ,  1,1)  ## = k_tp transpose
prior = sts.multivariate_normal(mu.flatten(), K_prior) 
#mean_test = np.dot(K_p , np.linalg.inv(K_prior))
mean_function = np.dot(np.dot(K_tp ,np.linalg.inv(K_prior).T) , prior )
covariance_function = K_train - np.dot(np.dot(K_tp ,np.linalg.inv(K_prior).T) , K_pt) 
</code></pre>