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
```