Последовательная выборка из условного многомерного нормального

Question

Я пытаюсь последовательно сэмплировать из предыдущего гауссовского процесса.

Проблема в том, что выборки в конечном итоге сходятся к нулю или расходятся до бесконечности.

Я использую основные условия, описанные, например, здесь

Примечание: функция kernel(X,X) возвращает квадрат экспоненциального ядра с изометрическим шумом.

Вот мой код:

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

Вот пример с 32 семплами, где он сворачивается:

введите описание изображения здесь

Вот пример с 34 образцами, где он взрывается:

введите описание изображения здесь

(для этого конкретного ядра 34 - это количество выборок, при котором (или больше) выборки начинают расходиться.