# Predictions from Latent GP in hierarchical model appear not conditioned on anything

**URL:** <https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309>\
**Category:** Questions\
**Created:** [June 22, 2020, 10:11pm UTC](https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309 "2020-06-22T22:11:34Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![BioGoertz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/biogoertz/32/4489_2.png) [@BioGoertz](https://discourse.pymc.io/u/BioGoertz)\
**Post date:** [June 22, 2020, 10:11pm UTC](https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309/1 "2020-06-22T22:11:34Z")

</div>

I’m using `gp.Latent` as a (top-level) prior in a hierarchical GLM to estimate parameter values in a mechanistic model of timeseries data (experimental instrument readings). The parameter estimates are well-defined and reasonable, but when I go to predict _new_ parameters using `gp.conditional` the GP draws don’t appear conditioned on the estimates. Even supplying the posterior estimates as `given` I just get flat predictions with broad, uniform variance. However, if I then use `gp.Marginal` to regress a GP on the posteriors using the same priors on the hyperparameters, the GP draws fit the data nicely. I’ve put a “toy” example in [this colab notebook](https://colab.research.google.com/drive/1PF3ni-8zU4nYCcP2UPW1aK9kDSHNM7Tm?usp=sharing). (Sorry, it takes ~30-60 minutes to run on colab, it was the simplest I could make it while retaining the structure.)

The “true” parameter values (as estimated without the GP) are very close to the fake ones I generate in the notebook. (Note, despited the appearance of dim2, I don’t want to use a periodic kernel; I don’t expect to observe much beyond the extents of this plot.)

 ![True parameter values](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/0/04b4d3c5458435cb1ffae6b8aea941e106e05812.png)

Fake experimental measurements, ~50 timepoints each:

 ![Reaction Measurements](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/a/ace4a0bdadf02806668ba5ea87384cc30a4712dd.png)

The model looks like this. I have 14 groups (reaction conditions), each of which has 2 replicates. Each group is assumed to have a typical parameter value from which individual reactions can vary slightly. The groups are defined by two components, called dim1 and dim2 below. The parameter `τ` (under an appropriate transform) varies very close to linearly along dim1, but has more complex behavior along dim2, hence the GP. My actual model has four parameters, each estimated with a similar hierarchical GLM structure, albeit without a GP component (for now).

```python
S_obs = np.hstack(data.sig)
grp = data.group.values
n_grps = len(set(grp))
n_pts = data.sig.map(len)
rxn = np.hstack([int(i)*np.ones(n, dtype=int) for i,n in enumerate(n_pts)])
n_rxns = len(set(rxn))
t = np.hstack(data.time)

dim2_lst = np.sort(list(set(dim2)))[:,None]
dim2_idx = [list(dim2_lst).index(v) for v in dim2]
dim2_new = np.linspace(dim2_lst.min()-0.25,dim2_lst.max()+0.25,100)[:,None]
coord1 = theano.shared(dim1)
coord2 = theano.shared(dim2)

with pm.Model() as latent_model:
    # Dim1 behaves nearly linearly
    τ_lgt_β_dim1 = pm.Normal('τ_lgt_β_dim1', mu=0, sigma=0.5)
    
    # Dim2 behavior is more complicated, requiring a GP
    ℓ = pm.InverseGamma('ℓ', mu=1.7, sigma=0.5)
    η = pm.Gamma('η', alpha=2, beta=1)
    cov = η**2 * pm.gp.cov.ExpQuad(1, ℓ)
    gp = pm.gp.Latent(cov_func=cov)
    τ_lgt_β_dim2 = gp.prior('τ_lgt_β_dim2', X=dim2_lst)
    
    ## Group-level hyperpriors
    # Expected value of τ based on dim1 and dim2
    τ_lgt_μ_grp = pm.Deterministic('τ_lgt_μ_grp', τ_lgt_β_dim1*coord1 + τ_lgt_β_dim2[dim2_idx])
    # Non-centered parameterization of actual value of τ for each group
    τ_lgt_μ_grp_σ = pm.Exponential('τ_lgt_μ_grp_σ', lam=1)
    τ_lgt_μ_grp_nc = pm.Normal('τ_lgt_μ_grp_nc', mu=0, sigma=1, shape=n_grps)
    τ_lgt_μ = pm.Deterministic('τ_lgt_μ', τ_lgt_μ_grp_nc*τ_lgt_μ_grp_σ+τ_lgt_μ_grp)

    # Non-centered parameterization of actual value of τ for each reaction
    τ_lgt_σ = pm.Exponential('τ_lgt_σ', lam=1)
    τ_lgt_z_nc = pm.Normal('τ_lgt_z_nc', mu=0, sigma=1, shape=n_rxns)
    τ_lgt_z = pm.Deterministic('τ_lgt_z', τ_lgt_z_nc*τ_lgt_σ+τ_lgt_μ[grp])
    # Un-transform τ back to the natural scale
    τ_lgt = pm.Deterministic('τ_lgt', τ_lgt_z*scale+shift)
    τ = pm.Deterministic('τ', pm.invlogit(τ_lgt))
    
    σ = pm.Exponential('σ', lam=1)
    
    # Simplified mechanistic model
    S = 1 / (1+pm.math.exp(t/τ[rxn]))
    S_lik = pm.StudentT('S_lik', nu=3, mu=S, sigma=σ*np.exp(-7), observed=S_obs)

```

After running `pm.sample()`, I draw predictions from the Latent GP:

```python
with latent_model:
    τ_new = gp.conditional('τ_new', Xnew=dim2_new, given={'X':dim2_lst,
                                                      'y':trace['τ_lgt_β_dim2'].mean(axis=0)[:,None],
                                                      'noise':trace['τ_lgt_β_dim2'].std(axis=0)[:,None],
                                                      'gp':gp})

    gp_draws = pm.sample_posterior_predictive([trace], vars=[τ_new], samples=1000)

```

The parameter estimates are very close to their true values. The estimates of `τ_lgt_β_dim2` along this axis are shown as white circles below, with the respective `τ_lgt_μ` values as blue circles. However, `gp.conditional` is giving what appears to be a GP conditioned on nothing at all.

 ![Latent GP Predictions](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/a/ae0562f2db8b7366104e7fef1022b90afa7afae1.png)

On the other hand, if I directly regress on the posterior estimates of `τ_lgt_β_dim2`, the result is as expected:

```python
with pm.Model() as marginal_model:
    ℓ = pm.InverseGamma('ℓ', mu=1.7, sigma=0.5)
    η = pm.Gamma('η', alpha=2, beta=1)
    
    cov = η**2 * pm.gp.cov.ExpQuad(1, ℓ)
    gp = pm.gp.Marginal(cov_func=cov)

    grp_τ = gp.marginal_likelihood('grp_τ', X=x, y=y, noise=e)

    mp = pm.find_MAP()
    trace_marginal = pm.sample()

with marginal_model:
    τ_new = gp.conditional('τ_new', dim2_new)

    map_samples = pm.sample_posterior_predictive([mp], vars=[τ_new], samples=1000)
    full_samples = pm.sample_posterior_predictive([trace_marginal], vars=[τ_new], samples=1000)

...

```

 ![Marginal GP Predictions](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/c/c425a16aed93abd7f7a0dadc1e3e0c97ab2f8a78.jpeg)

_ **My main question is, how do I get better predictions from the latent GP? Is this difference between `Latent` and `Marginal` implementations expected, or have I made a mistake?** _

As a side note, please let me know if you have suggestions for analyzing this sort of timeseries data with a mechanistic model better. The point-by-point likelihood evaluation seems inefficient (my real dataset ~11000 points, and I expect it to grow by two or three fold) and naive, ignorant of the inherent correlation in points (and therefore likelihoods) resulting from a given set of parameters.

---

<div class="post-metadata">

**Author:** ![AlexAndorra](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/alexandorra/32/9142_2.png) [@AlexAndorra](https://discourse.pymc.io/u/AlexAndorra)\
**Post date:** [June 23, 2020, 10:35am UTC](https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309/2 "2020-06-23T10:35:18Z")

</div>

Hi John,  
Welcome, and thanks for this detailed question!

Using `Latent` and `gp.prior` with a `StudentT` seems relevant to me, as shown in [this example](https://docs.pymc.io/notebooks/GP-Latent.html).  
However, I think @bwengals will be able to help you more quickly than I on this 😉

---

<div class="post-metadata">

**Author:** ![BioGoertz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/biogoertz/32/4489_2.png) [@BioGoertz](https://discourse.pymc.io/u/BioGoertz)\
**Post date:** [June 23, 2020, 11:01am UTC](https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309/3 "2020-06-23T11:01:59Z")

</div>

Thanks for responding quickly! Looking at that example more closely, it seems they first sample draws from the GP hyperprior at the _X_ coordinates in the dataset. That works for me, though the result is just linearly-interpolated retrodictions:

 ![Latent GP Retrodictions](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/3/39b42eb4de7be72dd965e86db58b8d1e2ef49f21.jpeg)

The `conditional` method is where things go haywire.

(As a side note, the second application in that example samples from `conditional`, but then just plots the retrodictions. Maybe that should be more consistent?)

---

<div class="post-metadata">

**Author:** ![BioGoertz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/biogoertz/32/4489_2.png) [@BioGoertz](https://discourse.pymc.io/u/BioGoertz)\
**Post date:** [July 1, 2020, 3:50pm UTC](https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309/4 "2020-07-01T15:50:41Z")

</div>

Apologies for the bump, but does anyone have suggestions here?

---

<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:** [July 3, 2020, 10:01pm UTC](https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309/5 "2020-07-03T22:01:46Z")

</div>

I tried running your collab, made a few tweaks and got:

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/1/19d64ec86f25d128486612e638f0009a33744d0b.jpeg)

Does this look right…? All I changed was:

```
with latent_model:
    τ_new = gp.conditional('τ_new', Xnew=dim2_new)
    gp_draws = pm.sample_posterior_predictive(trace, vars=[τ_new])
```

---

<div class="post-metadata">

**Author:** ![BioGoertz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/biogoertz/32/4489_2.png) [@BioGoertz](https://discourse.pymc.io/u/BioGoertz)\
**Post date:** [July 6, 2020, 7:55pm UTC](https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309/6 "2020-07-06T19:55:59Z")

</div>

Yup, that worked!

 ![Latent GP Predictions](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/f/f8133c790461f4519549195632381b910d41ba0d.jpeg)

All it took was a pair of square brackets… I guess I had looked at the examples for using MAP estimates first, in which you need to pass the MAP result as a single-element list. Anyways, thanks a lot for taking the time to dig through and find that!

---

<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:** [July 6, 2020, 9:15pm UTC](https://discourse.pymc.io/t/predictions-from-latent-gp-in-hierarchical-model-appear-not-conditioned-on-anything/5309/7 "2020-07-06T21:15:27Z")

</div>

Great! I’m not sure why exactly the brackets messed things up… I’ll double check that. I also noticed that the docstring on the `condition` method for `Latent` were a bit wrong (they look incorrectly copy/pasted from `Marginal`).

The conditioning would be on `gp`, or `X`, and `f`. Not on `y` or `noise`, since those certainly aren’t part of the `Latent` GP model. Sorry for the confusion!