Using a GP as Covariance Matrix for an MvNormal

I think it is working as expected. Here’s an explanation.

  • gp.prior is putting a N-dimensional re-parameterized multivariate normal prior where N are number of samples (not features) in your Dist_mat array. This means SIGMA in your model is going to be N dimensional sample from a multivariate normal. Here, N is 151.

I think you want to fit a model on Dist_mat to get a distribution over B_ observed values. You can do that in the following way.

class MyMean(pm.gp.Mean):
    def __init__(self, a, b):
        self.a = a
        self.b = b

    def __call__(self, X):
        return self.a + self.b[0]*M + self.b[1]*G


with pm.Model() as m:
    etasq = 1. + pm.HalfNormal('etasq', 0.25)
    rhosq = 3. + pm.HalfNormal('rhosq', 0.25)
    a = pm.Normal("a", 0., 1.)
    b = pm.Normal("b", 0., 0.5, shape=2)

    cov = etasq * pm.gp.cov.ExpQuad(input_dim=1, ls=rhosq)
    mu = MyMean(a, b)
    gp = pm.gp.Marginal(mean_func=mu, cov_func=cov)
    B = gp.prior('B', X=Dist_mat**2, y=B_)
    trace = pm.sample(1000)

I don’t know what M and G are in your model but probably you want to make the mean function a function of your inputs X which is Dist_mat in your model.

1 Like