I’m have a 2D dataset containing a linear relationship + gamma noise, like what follows:
import numpy as np length = 200 rng = np.random.default_rng(1337) x = np.linspace(1e-2, 1, num=length) true_regression_line = 5 * x + 4 y = true_regression_line + rng.gamma(0.5, 0.2, size=length)
I can fit an ordinary linear regression with Gaussian noise with the following code:
with pm.Model() as guassian_model: sigma = pm.HalfCauchy("sigma", beta=10) intercept = pm.Normal("intercept", 0, sigma=20) slope = pm.Normal("slope", 0, sigma=20) likelihood = pm.Normal("y", mu=intercept + slope * x, sigma=sigma, observed=y) map_estimate = pm.find_MAP()
resulting in the following fit:
So far so good. I’m now trying to improve the fit by recognising that the noise is Gamma-distributed, and replacing the Gaussian likelihood from the model above. I found the following code snippet that claims to achieve this:
with pm.Model() as gamma_model: a = pm.Normal("a", 0, 0.5) b = pm.Normal("b", 0, 0.5) shape = pm.Uniform("shape", 0, 50) μ = pm.math.exp(a + b * x) likelihood = pm.Gamma("y", alpha=shape, beta=shape / μ, observed=y) map_estimate = pm.find_MAP(progressbar=False)
which results in the following fit, which I find a bit disappointing, although I may be missing the point:
I was expecting a line closer to the true regression line, and suspect that there’s something odd with
gamma_model, as the Gamma noise must be >0 and that doesn’t seem to be the case with the fitted line in the plot. In all honesty, I don’t follow all the details of
gamma_model—is it using the appropriate link function for the Gamma distribution?
Thanks for your time, any tips appreciated!