Understanding how a term in GP model is used

Questions
1

I’m following the Pymc repo of chapter 14 of McElreath’s Statistical Rethiking.

In cell 61 (Code 14.39), we have the following model and I don’t understand where rhosq is being used:

with pm.Model() as m14_8:
    a = pm.Exponential("a", 1.0)
    b = pm.Exponential("b", 1.0)
    g = pm.Exponential("g", 1.0)

    etasq = pm.Exponential("etasq", 2.0)
    ls_inv = pm.HalfNormal("ls_inv", 2.0)
    rhosq = pm.Deterministic("rhosq", 0.5 * ls_inv**2)

    # Implementation with PyMC's GP module:
    cov = etasq * pm.gp.cov.ExpQuad(input_dim=1, ls_inv=ls_inv)
    gp = pm.gp.Latent(cov_func=cov)
    k = gp.prior("k", X=Dmat)

    lam = (a * P**b / g) * at.exp(k[society])

    T = pm.Poisson("total_tools", lam, observed=total_tools)

    trace_14_8 = pm.sample(4000, tune=2000, target_accept=0.99, random_seed=RANDOM_SEED)

I’m new to Gaussian Processes and have gone through to understand the ExpQuad and Latent syntax. However, I’m not following - rhosq doesn’t appear to be used anywhere after its definition. I see that it’s dependent on ls_inv, and ls_inv is used to define the covariance function.

Note that in the latest Pymc version, you would need to replace from aesara import tensor as at with import pytensor.tensor

1 Like
2

Hi Matsuo, I just ran into the same mystery and spent some time puzzling through it.

You’re right that rhosq isn’t actually doing anything in the model. ls_inv does the whole job. There are a couple of different ways to express the exponential quadratic kernel. McElreath uses one. Pymc uses another. Check the docs for exponential quadratic and compare it to what is in the textbook to see how they compare. Anyway, I think the writer of the GP example wanted to be able to get a posterior on rho squared to benchmark their work to McElreath’s, even though the pymc model only wants ls_inv.

2 Likes
3

Daniel, thanks for looking into this in-depth and clearing this up. Yes, you’re probably right they wanted to compare rhosq posterior to the book’s result.

1 Like
4

Typically with PyMC if you see pm.Deterministic thats a good indication that the modeler wants to look at the variable after sampling in az.summary or some other manner. Sometimes that deterministic is used “in the model”, sometimes its “outside”