AnswerBun.com

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        
    else:
        # 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 + np.dot(C12, np.dot(C22_, (y2 - u2)))
        C22_xC21 = np.dot(C22_, C21)
        C_minus = np.dot(C12, 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 AnswerBun.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP