Help with fitting Gamma Distribution

Questions
1

Hi all,
New to pymc3 so thanks in advance for your patience.

I have a Poisson process that gives me the following distribution.
I would like to fit a gamma distribution and use the distributions of the alpha and beta parameters.

I can get the alpha and beta parameters from scipy.stats.gamma.fit but would like the distributions of these parameters.

Any advice on how to respecify this model?

image
image.png746×366 15.3 KB

Using the parameterizations for gamma here:
https://juanitorduz.github.io/intro_pymc3/

My model:

model = pm.Model()

with model:
# Prior distribution for mu.
mu = pm.Gamma('mu', alpha=100, beta=1.0/5)

# Prior distribution for sigma
sigma = pm.Exponential('sigma', 100.0)

# Parametrization for the (alpha) shape parameter.
alpha =  mu**2/sigma**2

# Parametrization for the scale parameter.
beta = mu/sigma**2  

y_obs = pm.Gamma('y_obs', alpha=alpha, beta=beta, observed=df['bw'].values)

trace = pm.sample(draws=2000)

Error:

ValueError: Bad initial energy: inf. The model might be misspecified.
2

May be it is due to the division by sigma which has an exponential prior and during the sampling it may sample 0s (which leads to division by 0 --> inf).

Instead of Exponential try using an prior that behave according to the expected behavior.

1 Like
3

The main problem is that the Gamma distribution can only have positive observed values, that is no observed can be equal to zero. I suggest you just change the values of the observed from 0 to 1e-3. That, along with what @Nadheesh said about the exponential having is most likely value equal to zero, should help you.

4

Thanks! I realize I used a 1e-4 offset in my scipy implementation and forgot to apply that in my model above!

5

Thanks all for the feedback above. for zero values, I added 1e-4.

Here’s the working model:

model = pm.Model()

with model: 

# alpha
alpha = pm.Exponential('alpha', 10)

# beta
beta = pm.Exponential('beta', 100)

g = pm.Gamma('g', alpha=alpha, beta=beta, observed=bw)
    
trace = pm.sample(2000)
1 Like