Censored Gamma giving reciprocal betas?

Questions version agnostic
1

I used a censored Gamma distribution for my outcome and got what seems to be the reciprocal of beta in the posterior.

I’m not sure if there’s a problem with my understanding of the statistics, understanding of PyMC, a limitation of the sampling, or a bug somewhere. I’m using pymc==5.5.0

I’m aware that the Gamma distribution can be parametrized by the shape, plus the beta (the “rate”) or theta == 1/beta (the “scale”).

Here is the censored code:

import pymc as pm
import matplotlib.pyplot as plt

n = 1000
g_dist = pm.Gamma.dist(alpha=24, beta=6)
y = pm.draw(g_dist, n, random_seed=3)
ceiling = 11  # I had it set to actually censor, but same results w/ 11
y = [x if x < ceiling else ceiling for x in y]

# sum([x == ceiling for x in y])  # count of censored
# _, ax = plt.subplots()
# ax.set_ylim(0,1)
# ax.hist(y, density=True)

with pm.Model() as model:
    # Data
    y_data = pm.ConstantData("y_data", y)

    # Priors
    alpha = pm.Normal("alpha (shape)", mu=22, sigma=3)
    beta = pm.Normal("beta (rate)", mu=6, sigma=2)

    # lik
    Y_obs = pm.Censored("Y_obs", 
                        dist=pm.Gamma.dist(alpha=alpha, beta=beta),
                        lower=None,
                        upper=ceiling,
                        observed=y_data)

with model:
    idata = pm.sample(random_seed=3)

beta_mean = idata.posterior["beta (rate)"].mean().to_numpy()
alpha_mean = idata.posterior["alpha (shape)"].mean().to_numpy()
alpha_mean, beta_mean, 1 / beta_mean

And the non-censored version:

import pymc as pm
import matplotlib.pyplot as plt

n = 1000
g_dist = pm.Gamma.dist(alpha=24, beta=6)
y = pm.draw(g_dist, n, random_seed=3)

with pm.Model() as model:
    # Data
    y_data = pm.ConstantData("y_data", y)

    # Priors
    alpha = pm.Normal("alpha (shape)", mu=22, sigma=2)
    beta = pm.HalfNormal("beta (rate)", sigma=4)

    # lik
    Y_obs = pm.Gamma("Y_obs", alpha=alpha, beta=beta, observed=y_data)

with model:
    idata = pm.sample(random_seed=3)

beta_mean = idata.posterior["beta (rate)"].mean().to_numpy()
alpha_mean = idata.posterior["alpha (shape)"].mean().to_numpy()

alpha_mean, beta_mean
2

Can you test with a more recent version pf PyMC?

3

Oh, duh. Yeah it’s giving the right answer now. :+1: Thanks.

4

Probably some bug long fixed. :slight_smile: