linear_GP_test.py

import numpy as np

n = 40
d = 20

sigma = .05
sigma2 = sigma**2

X = np.random.normal(size=(n, d))
Y = X.sum(axis=1, keepdims=True)
Y = Y + np.random.normal(size=Y.shape) * sigma

S = np.identity(d) + np.random.normal(size=(d, d)) * 0.0
S = 0.5 * (S + S.T)

Sinv = np.linalg.inv(S)

A = X.T.dot(X) / sigma2 + Sinv
Ainv = np.linalg.inv(A)

w_star = Ainv.dot(X.T).dot(Y) / sigma2

x = np.random.normal(size=(d, 1))
y = x.sum() + np.random.normal() * sigma

B = A + x.dot(x.T) / sigma2
Binv = np.linalg.inv(B)

print("Ratio Var(Book) / Var(My):", (1 - x.T.dot(Binv).dot(x) / sigma2) / sigma2 * (x.T.dot(Ainv).dot(x)).squeeze())

n_mc = 1000000
w_mc = np.random.multivariate_normal(mean=w_star.squeeze(), cov=Ainv, size=n_mc)
y = np.random.normal(size=n_mc) * sigma + w_mc.dot(x).squeeze()
print("Ratio Var(MC) / Var(My):", (y.var() * (1 - x.T.dot(Binv).dot(x) / sigma2) / sigma2).squeeze())