# Pointwise Log Likelihood For Gaussian Processes

**URL:** <https://discourse.pymc.io/t/pointwise-log-likelihood-for-gaussian-processes/16423>\
**Category:** v5\
**Created:** [January 21, 2025, 5:09pm UTC](https://discourse.pymc.io/t/pointwise-log-likelihood-for-gaussian-processes/16423 "2025-01-21T17:09:43Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![rob00](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@rob00](https://discourse.pymc.io/u/rob00)\
**Post date:** [January 21, 2025, 5:09pm UTC](https://discourse.pymc.io/t/pointwise-log-likelihood-for-gaussian-processes/16423/1 "2025-01-21T17:09:43Z")

</div>

I am trying to compare different GP models using `az.loo` or `az.waic`. In PyMC these methods work by computing the pointwise log-likelihood with `pm.compute_log_likelihood` function.  
The problem i am facing is that with GP models the compute\_log\_likelihood does not return the values for each element.  
I was expecting the loglikelihood group to have shape (chain, draw, y\_dim\_0) instead it has (chain,draw).  
Is there something i am missing?  
I will leave a reproducible example from the Marginal Likelihood Implementation notebook.

```python
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import scipy as sp

# set the seed
np.random.seed(1)

n = 100 # The number of data points
X = np.linspace(0, 10, n)[:, None] # The inputs to the GP, they must be arranged as a column vector

# Define the true covariance function and its parameters
ell_true = 1.0
eta_true = 3.0
cov_func = eta_true**2 * pm.gp.cov.Matern52(1, ell_true)

# A mean function that is zero everywhere
mean_func = pm.gp.mean.Zero()

# The latent function values are one sample from a multivariate normal
# Note that we have to call `eval()` because PyMC3 built on top of Theano
f_true = np.random.multivariate_normal(
    mean_func(X).eval(), cov_func(X).eval() + 1e-8 * np.eye(n), 1
).flatten()

# The observed data is the latent function plus a small amount of IID Gaussian noise
# The standard deviation of the noise is `sigma`
sigma_true = 2.0
y = f_true + sigma_true * np.random.randn(n)

with pm.Model() as model:
    ell = pm.Gamma("ell", alpha=2, beta=1)
    eta = pm.HalfCauchy("eta", beta=5)

    cov = eta**2 * pm.gp.cov.Matern52(1, ell)
    gp = pm.gp.Marginal(cov_func=cov)

    sigma = pm.HalfCauchy("sigma", beta=5)
    y_ = gp.marginal_likelihood("y", X=X, y=y, sigma=sigma)

with model:
    marginal_post = pm.sample(nuts_sampler="pymc", idata_kwargs={"log_likelihood": True})

az.loo(marginal_post)

```

The loo computation of course provide the following Userwarning:  
‘The point-wise LOO is the same with the sum LOO, please double check the Observed RV in your model to make sure it returns element-wise logp.’

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 21, 2025, 9:56pm UTC](https://discourse.pymc.io/t/pointwise-log-likelihood-for-gaussian-processes/16423/2 "2025-01-21T21:56:51Z")

</div>

CC @bwengals

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [January 23, 2025, 4:41am UTC](https://discourse.pymc.io/t/pointwise-log-likelihood-for-gaussian-processes/16423/3 "2025-01-23T04:41:17Z")

</div>

This would be expected because `pm.gp.Marginal.marginal_likelihood` is just a `pm.MvNormal` under the hood. It’s evaluating your `y` data as if it’s one single data point – one draw from an multivariate normal.

You’ll need to either calculate the log-likelihoods manually, ~~or resample your model using `pm.gp.Latent` with separate `Normal` likelihood.~~ what @jessegrabowski said below

You’ll find that the sum of the log-likelihoods calculated this way will equal the log-likelihood from `pm.gp.Marginal`.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 23, 2025, 4:58am UTC](https://discourse.pymc.io/t/pointwise-log-likelihood-for-gaussian-processes/16423/4 "2025-01-23T04:58:49Z")

</div>

You shouldn’t need to resample the whole model. Instead, write the `pm.gp.Latent` + `pm.Normal` version of the model, then inside that context call `pm.compute_log_likelihood`, passing in the idata from the marginalized version.

---

<div class="post-metadata">

**Author:** ![rob00](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@rob00](https://discourse.pymc.io/u/rob00)\
**Post date:** [February 6, 2025, 5:52pm UTC](https://discourse.pymc.io/t/pointwise-log-likelihood-for-gaussian-processes/16423/5 "2025-02-06T17:52:19Z")

</div>

First of all, thanks for the quick replies.  
I tested some of your suggestions.

1. 

> resample your model using `pm.gp.Latent` with separate `Normal` likelihood. The sum of the log-likelihoods calculated this way will equal the log-likelihood from `pm.gp.Marginal`.

Using the model from my first message i have written its Latent version (which calculate the correct elementwise log likelihood) but the sum of the likelihoods do not match. I have -3.37e+05 for the Latent model and -4.39e+05 for the Marginal model.I have set the random seed as suggested [here](https://discourse.pymc.io/t/how-to-set-a-seed-for-pm-sample/11497/2).  
The implementation of the Latent model is the following:

```auto
rng = np.random.default_rng(42)

with pm.Model() as model:
    ell = pm.Gamma("ell", alpha=2, beta=1)
    eta = pm.HalfCauchy("eta", beta=5)

    cov = eta**2 * pm.gp.cov.Matern52(1, ell)
    gp = pm.gp.Latent(cov_func=cov)

    f = gp.prior("f", X=X)
    
    sigma = pm.HalfCauchy("sigma", beta=5)
    y_ = pm.Normal('y', f, sigma, observed=y)

    latent_post = pm.sample(nuts_sampler="pymc", idata_kwargs={"log_likelihood": True}, random_seed=rng)

# i calculated the sum in this way
print(latent_post.log_likelihood.sum()) 
# compare it to the marginal model
print(marginal_post.log_likelihood.sum())

```

1. 

> inside that context call `pm.compute_log_likelihood`, passing in the idata from the marginalized version.

As suggested i did the following:

```auto
with pm.Model() as model:
    ell = pm.Gamma("ell", alpha=2, beta=1)
    eta = pm.HalfCauchy("eta", beta=5)

    cov = eta**2 * pm.gp.cov.Matern52(1, ell)
    gp = pm.gp.Latent(cov_func=cov)

    f = gp.prior("f", X=X)
    
    sigma = pm.HalfCauchy("sigma", beta=5)
    y_ = pm.Normal('y', f, sigma, observed=y)

    latent_post = pm.compute_log_likelihood(marginal_post)

```

and i have the following error:

```auto
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
/opt/conda/envs/dev/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, names)
   1474 variables[name] = self._variables[name]
   1475 except KeyError:
-> 1476 ref_name, var_name, var = _get_virtual_variable(
   1477 self._variables, name, self.sizes

KeyError: 'f_rotated_'

```

The only method that worked so far is to resample the model using the Latent implementation but it is very slow and not feasible in my application.

Do you have any ideas on how to solve this? My last approach would be calculate it manually. If you could give me an input also on that it would be very helpful.

I need also to ask you another thing related to this topic. I am using `pm.compute_log_likelihood` on a sparse gp. My understanding is that the sparse gp implementation is done using a `pm.Potential` which the `pm.compute_log_likelihood` ignores. Is that correct?  
Then how can i compute it? My guess is that the best way is to compute it manually because writing the Full Latent model and resampling would be too slow.  
In theory, at least from my understanding, this pointwise log likelihood (from the sparse model) would be calculated from the approximation of the likelihood given by FITC or VFE.

Thanks again

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [February 11, 2025, 8:21pm UTC](https://discourse.pymc.io/t/pointwise-log-likelihood-for-gaussian-processes/16423/6 "2025-02-11T20:21:04Z")

</div>

> [@rob00](#):
>
> I need also to ask you another thing related to this topic. I am using `pm.compute_log_likelihood` on a sparse gp. My understanding is that the sparse gp implementation is done using a `pm.Potential` which the `pm.compute_log_likelihood` ignores. Is that correct?  
> Then how can i compute it? My guess is that the best way is to compute it manually because writing the Full Latent model and resampling would be too slow.  
> In theory, at least from my understanding, this pointwise log likelihood (from the sparse model) would be calculated from the approximation of the likelihood given by FITC or VFE.

So I’d be careful here, because VFE is a variational _approximation_, so the likelihood isn’t calculated as part of the model, but a lower bound on it. FITC is also viewed as a likelihood approximation, see [this paper](https://quinonero.net/Publications/tr-2007-124.pdf).

If I were you, I’d check [Section 5.4.2 in GPML](https://gaussianprocess.org/gpml/chapters/RW5.pdf) for the analytic calculation.

For your point 2, that error might be happening because you have `reparameterize = True` (the default) in the call `f = gp.prior("f", X=X)`. Try `reparameterize=False`.