How does one sample from a tempered posterior?

Questions v5
1

Hello! I’m trying to sample from a tempered posterior distribution. As I understand, I should be doing this using pm.Potential, but no matter how I arrange things I can’t get the temperature parameter to have any effect on the sampled distribution.

Here’s the code I’m using to test things:

import numpy as np
import pymc as pm
import pytensor.tensor as pt
import arviz as az

np.random.seed(42)
X = np.random.normal(loc=0, scale=1, size=500)

def tempered_gaussian(
        observations,
        beta=1.0,
    ):
    with pm.Model() as model:
        ## Priors
        mu = pm.Normal("mu", mu=0, sigma=10)
        sigma = pm.HalfNormal("sigma", sigma=10)

        likelihood = pm.Normal("likelihood", mu=mu, sigma=sigma, observed=observations)
        tempered_ll = pm.Potential("tempered_ll", beta*likelihood)
                                  
    return model

model = tempered_gaussian(X, beta=1/np.log(X.size))
with model:
    idata = pm.sample(
        draws=10000,
        tune=1000,
        chains=4,
    )

az.plot_trace(idata.posterior)

No matter what I set beta to, it doesn’t change the output distribution. Even e.g. beta=0 doesn’t seem to do anything. What am I doing wrong?

Thanks in advance.

2

The variable you call likelihood is not a likelihood, it’s just another random variable that happens to have associated observed data. It will be used to compute a likelihood term when pymc transforms your forward model into a logp program, but you’re not there yet.

If you want to directly work the the logp, you can do something like:

obs_dist = pm.Normal.dist(mu=mu, sigma=sigma)
tempered_ll = pm.Potential('tempered_ll', beta * pm.logp(obs_dist, observations))

Or you can use CustomDist to make a new distribution with your desired logp, then it will work more like a normal pymc distribution (depending on how you set up the CustomDist – see its docstring for an exhaustive discussion)

Finally, if you’re actually interested in the normal case and this isn’t just a simplified example, you can also just set sigma = sigma / sqrt(beta), since F(x | \mu, \sigma) \propto \exp \left [ -\frac{(x - \mu)^2}{2\sigma^2} \right ] , so (F^\beta) \propto \exp \left [ -\frac{\beta}{2\sigma^2} (x - \mu)^2 \right ] is clearly also a normal distribution with the same mean but a new variance.

2 Likes