Collectives™ on Stack Overflow

Find centralized, trusted content and collaborate around the technologies you use most.

Teams

Q&A for work

Connect and share knowledge within a single location that is structured and easy to search.

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.  
                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, 2018 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, 2018 at 1:44
                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, 2018 at 16:13
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)
                That's an interesting idea but, yes, not being able to pass observed data would be limiting. (Reference: github.com/pymc-devs/pymc3/issues/1971)
– Mike
                Feb 22, 2018 at 1:27
        
Thanks for contributing an answer to Stack Overflow!
Please be sure to answer the question. Provide details and share your research!
But avoid …
Asking for help, clarification, or responding to other answers.
Making statements based on opinion; back them up with references or personal experience.
To learn more, see our tips on writing great answers.