Calculating marginal likelihood of GP Model

Questions
1

Hi, I am using the marginal likelihood implementation of PyMC3 GP. How can I calculate the marginal likelihood for model comparison? My understanding is that model.logp(y) would give me the log posterior.

A bit unrelated, I tried calculating MLE for the GP model using flat priors and pm.find_MAP() but it doesn’t work because of inf or Nan values. The initial idea was to calculate BIC but as that cannot be done I want to calculate marginal likelihood for model comparison.

I was seeing bridge sampling notebook @junpenglao but then read The Harmonic Mean of the Likelihood: Worst Monte Carlo Method Ever. I know bridge sampling wouldn’t work here. As far as BIC is concerned, BIC GPy, here by Neil comments that theoretically it’s invalid for GPs as they assume correlation among parameters. After a lot of reading, I have got confused.

Any suggestion to compute BIC or Marginal likelihood would be appreciated!

Best Regards,
Rohan

2

Neal is right. And in most cases, when you have a complex model, computing marginal likelihood is extremely difficult. I would not trust anything and advise you to choose other methods to do model comparison (loo or waic).

3

The nan in find_map() might related to the bad initial value, try changing the start kwarg of the find_map() call.

4

Hi,

Thanks for your comments on metrics.

Regarding the MLE:

lengthscale_true = 1
mean = pm.gp.mean.Zero()
cov = pm.gp.cov.ExpQuad(1,lengthscale_true)

n = 200
X_true = np.linspace(0,5,n)[:,None]
y_true = np.random.multivariate_normal(mean(X_true).eval(),cov(true).eval() + 1e-8*np.eye(n), 1).flatten()

noiseobs = 0.5
X_obs = X_true
y_obs = y_true + noiseobs * np.random.randn(n)

with pm.Model() as standardmodel:
noise_obs = pm.Flat(‘noise_obs’)
lengthscale = pm.Flat(‘lengthscale’)

covariance = pm.gp.cov.ExpQuad(1,lengthscale)
mean = pm.gp.mean.Zero()

gp_obs = pm.gp.Marginal(mean_func=mean, cov_func=covariance)

yobs = gp_obs.marginal_likelihood("yobs", X=X_obs, y=y_obs, noise=noise_obs)
start={'noise_obs': 0.5, 'lengthscale': 1}
mp = pm.find_MAP(start=start)

The error occurs on yobs line. The error is ValueError: array must not contain infs or NaNs. It can be reporoduced using the above code. I went through the trace but couldn’t find the problem. It might be happening due to Flat priors.

5

The length scale and noise only takes positive value - you should use HalfFlat instead.

    noise_obs = pm.HalfFlat('noise_obs')
    lengthscale = pm.HalfFlat('lengthscale')
1 Like
6

So silly to miss that on my part. Sorry for asking so many questions!

Just to confirm the output of logp by MAP is logp of the model. So the optimizer is maximizing log the probablity of the model.

logp = -220.68, ||grad|| = 149.66: 100%|██████████| 11/11 [00:00<00:00, 97.39it/s]

mp

{'noise_obs_log__': array(-0.7523434),
 'lengthscale_log__': array(-0.14632332),
 'noise_obs': array(0.47126091),
 'lengthscale': array(0.86387835)}

So, how do we get the log MLE ?

yobs.logp(mp)

Thanks a lot!

7

You can do model.logp(mp) to get the model logp (up to a constant) at the point of MLE.