Using Gamma distribution parameters from SciPy

Hello,

I like to use input parameters that have been fitted to Gamma distributions with SciPy in a PyMc3 model. The SciPy implementation of the Gamma distribution has three parameters (a, loc and scale, see SciPy Docs for gamma). For an PyMc3 model, I’m limited to either alpha and beta or mu and sigma, according to the docs.

Is there any way to implement a Gamma distribution with the three parameters from SciPy?

I tried this as a work-around, based on a snippet from the getting started notebook:

import pymc3 as pm
import scipy.stats as st
with pm.Model() as model:
    a = 7.43
    loc = 124.57
    scale = 3.15
    pm.DensityDist('gamma', lambda value, a=a, loc=loc, scale=scale: st.gamma.logpdf(value, a, loc, scale), testval=150)

But that just gives me TypeError: Variables do not support boolean operations., which is entirely cryptic to me. Could you please help me with this issue?

Thanks in advance!

I see two possible issues with your code:

  1. I think the custom function used in the DensityDist can only take one argument (value) Nope, this seems fine the way you did it.
  2. You have to explicitly write down the logpdf yourself, so that it works with theano.

This is what I would have tried:

import pymc3 as pm
import scipy.stats as st

def my_gamma(a, loc, scale):
   def logp(value):
       return ... # (explicit equation for lopg using a, loc, scale and value, goes here)

   return logp

with pm.Model() as model:
    a = 7.43
    loc = 124.57
    scale = 3.15
    pm.DensityDist('gamma', my_gamma(a=a, loc=loc, scale=scale), testval=150)

PS: Your error message is probably related to trying to use theano variables inside the scipy function, which is why you have to write the formula yourself.

I am not completely sure, so hopefully someone else will be able to chime in :slight_smile:

Thank you for your reply.

I took your advice and added the @as_op decorator from Theano

import theano.tensor as tt
from theano.compile.ops import as_op
...
def my_gamma(a, loc, scale):
    @as_op(itypes=[tt.dscalar], otypes=[tt.dscalar])
    def logp(value):
        return np.array(st.gamma.logpdf(value, a, loc, scale))
    return logp

That works in general, but is slower due to the wrapping. As you said, a Theano implementation would be better. So far, I had no luck finding the implementation in the SciPy code though and both my mathematics and Theano skills are not enough to do it myself.

I don’t know what is the goal but if it is just parameterization, you could use the ShiftedGamma that was defined here (code below). It extends the PyMC3 Gamma distribution to be able to receive a parameter equivalent to loc for SciPy.

class ShiftedGamma(pm.Gamma):
    def __init__(self, alpha, beta, shift, *args, **kwargs):
        transform = pm.distributions.transforms.lowerbound(shift)
        super().__init__(alpha=alpha, beta=beta, *args, **kwargs,
                         transform=transform)
        self.shift = shift
        self.mean += shift
        self.mode += shift
        
    def random(self):
        return super().random() + self.shift
    
    def logp(self, x):
        return super().logp(x - self.shift)

Then just make sure that you transform scale to beta and you are good to go.

with pm.Model() as model:
    
    alpha = pm.HalfNormal('alpha', 1)
    scale = pm.HalfNormal('scale', 1)
    loc = pm.HalfNormal('loc', 10)
    
    beta = 1/scale
    
    # likelihood
    y_g = ShiftedGamma("y_g", 
                          alpha=alpha, 
                          beta=beta, 
                          shift=loc, observed=y)
    
    trace = pm.sample(
        2000,
        tune=2000
    )

Hope it helps!

3 Likes

Thank you very much, that’s it! I managed to do a similar Implementation for the Beta distribution myself.