3

I'm trying to create a beta RV that behaves like the one in scipy. Namely, it takes an alpha, beta, loc, and scale. I followed the discussion on creating a shifted gamma, and came up with something that looks like this:

import pymc3 as pm

class SSBeta(pm.Beta):
    def __init__(self, alpha, beta, loc, scale, *args, **kwargs):

        # Not sure what this is doing but it won't work without it. 
        transform = pm.distributions.transforms.lowerbound(loc)

        super().__init__(alpha=alpha, beta=beta, *args, **kwargs, transform=transform)

        self.scale = scale
        self.loc = loc
        self.mean += loc

    def random(self):
        return super().random()*self.scale + self.loc

    def logp(self, x):
        return super().logp((x - self.loc)/self.scale)

So I have two questions:

  • Is this implementation correct (random, logp)?
  • What's the purpose of the transform at the top of the class? I can't find anything useful in the docs and the code didn't help a ton.
Taku
  • 31,927
  • 11
  • 74
  • 85
Mike
  • 251
  • 1
  • 12
  • 1
    That implementation looks fine to me - you should not need a transform (which sometimes makes sampling easier - for example, allowing you to sample from all the reals instead of just the positives). Is something going wrong with it? I have posted an alternative as an answer. – colcarroll Feb 20 '18 at 19:26
  • When calling it with price = SSBeta("price", 2, 5, 20, 10) I get "ValueError: Bad initial energy: nan. The model might be misspecified." when I remove that transform line. (Using: joblib-0.11 pymc3-3.3 theano-1.0.1 tqdm-4.19.5) – Mike Feb 22 '18 at 01:44
  • 1
    Oh interesting - yes, `Beta` uses the `logodds` transform, which assumes all proposals will be in [0, 1] (and so throws the error, because price must be between 20 and 30 in your example). I would use instead `transform = pm.distributions.transforms.interval(loc, loc + scale)`. – colcarroll Feb 22 '18 at 16:13

1 Answers1

1

An alternative is using pm.Deterministic. This does not allow you to pass observed data, but might be what you want anyways?

mu, scale = 2, 4
with pm.Model():
    b = pm.Beta('b', alpha=3, beta=8)
    shifted_b = pm.Deterministic('shifted_b', scale * b + mu)

enter image description here

colcarroll
  • 3,632
  • 17
  • 25
  • That's an interesting idea but, yes, not being able to pass observed data would be limiting. (Reference: https://github.com/pymc-devs/pymc3/issues/1971) – Mike Feb 22 '18 at 01:27