Sequential sampling from Gaussian conditional not working

Data Science Asked by Jacob Holm on September 28, 2020

I’m trying to sequentially sample from a Gaussian Process prior.

The problem is that the samples eventually converge to zero or diverge to infinity.

I’m using the basic conditionals described e.g. here

Note: the kernel(X,X) function returns the squared exponential kernel with isometric noise.

Here is my code:

n = 32

x_grid = np.linspace(-5,5,n)

x_all = []
y_all = []
for x in x_grid:
    x_all = [x] + x_all
    X = np.array(x_all).reshape(-1, 1)
    # Mean and covariance of the prior
    mu = np.zeros((X.shape), np.float)
    cov = kernel(X, X)
    if len(mu)==1: # first sample is not conditional
        y = np.random.randn()*cov + mu        
        # condition on all previous samples
        u1 = mu[0]
        u2 = mu[1:]
        y2 = np.atleast_2d(np.array(y_all)).T
        C11 = cov[:1,:1] # dependent sample
        C12 = np.atleast_2d(cov[0,1:])
        C21 = np.atleast_2d(cov[1:,0]).T
        C22 = np.atleast_2d(cov[1:, 1:])
        C22_ = la.inv(C22)
        u = u1 +,, (y2 - u2)))
        C22_xC21 =, C21)
        C_minus =, C22_xC21) # this weirdly becomes larger than C!
        C = C11 - C_minus        
        y = u + np.random.randn()*C
    y_all = [y.flatten()[0]] + y_all

Here’s an example with 32 samples, where it collapses:

enter image description here

Here’s an example with 34 samples, where it explodes:

enter image description here

(for this particular kernel, 34 is the number of samples at which (or more) the samples start to diverge.

I’ve gone through this code so many times that I’m going blind – something must be fundamentally wrong with it but I just can’t see it.

Add your own answers!

Related Questions

Hot to use the formula model in t.test

1  Asked on May 14, 2021 by molitoris


How do I plot data in Octave?

1  Asked on May 14, 2021 by mrpythonic


What is the opposite of baseline?

2  Asked on May 12, 2021 by joe_mind


Extended ResNet

0  Asked on May 12, 2021 by coderhk


Ask a Question

Get help from others!

© 2023 All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP