A gamma GLM models the expected value of a gamma distribution as a linear function, but the gamma function itself is non-linear, so you end up with that curving in the data (due to the exponential link function). I guess what you want is a fixed mean with gamma noise, which is different. You could try something like this:
import pymc as pm
import arviz as az
def gamma_noise_regression(mu, alpha, beta, size=None):
return mu + pm.Gamma.dist(alpha=alpha, beta=beta, size=size)
with pm.Model() as model:
intercept = pm.Normal('intercept', sigma=10)
slope = pm.Normal('slope', sigma=10)
mu = intercept + slope * x
alpha = pm.Exponential('alpha', 1)
beta = pm.Exponential('beta', 1)
y_obs = pm.CustomDist('y_obs', mu, alpha, beta, dist=gamma_noise_regression, observed=y)
idata = pm.sample(init='jitter+adapt_diag_grad')
idata = pm.sample_posterior_predictive(idata, extend_inferencedata=True)
I’m using the powers of PyMC’s automatic logp inference to turn mu into the location parameter of a generalized gamma, and using that to model the data. I had to add some normal noise to mu when I generated the data, otherwise it was divergence city (this is common in models with noiseless components). I got this as my output:
