Specifying loc and scale on beta prior requires bounds, gives error

Questions version agnostic
1

I’m building a model where the likelihood is gamma and the alpha and beta priors have a LaPlace and beta distribution, respectively.

I need to specify the beta distribution with the loc and scale parameters that are specified in scipy, in addition to ‘a’ and ‘b’.

I found a post that allows me to specify the beta distribution with those additional parameters. However, when I specify as is, I get the following error:

SamplingError : Initial evaluation of model at starting point failed! Starting values: {‘alpha_lowerbound__’: array(0.), ‘beta_testing_interval__’: array(nan)}

So I attempted to do something similar to what I do with the LaPlace, placing a bound on the custom distribution:

import numpy as np
import pymc3 as pm

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

        transform = pm.distributions.transforms.interval(loc, loc + scale)

        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)

with pm.Model() as model1:
    bounded_laplace = pm.Bound(pm.Laplace, lower=0.0)
    alpha = bounded_laplace('alpha', mu = 51704326, b = 767593)

    bounded_beta = pm.Bound(SSBeta, lower=0.0)
    beta_param = bounded_beta(name = 'beta_testing',
                    alpha = 1.578, 
                    beta = 0.855, 
                    loc = 0.0004,
                    scale = 7.338e-5)
    pred = pm.Gamma("test_var", alpha = alpha, beta = beta_param, 
    observed=np.array([1.72e6, 1.58e6, 1.64e6, 1.59e6]))

    trace = pm.sample(400)

This is the error I get:

TypeError : pymc3.distributions.continuous.Beta.init() got multiple values for keyword argument ‘transform’

It does appear that when I change the scale to 1 or above, the model runs.

2

Can we see the model throws the initial evaluation failure?

3

Sure, it’s like this:

with pm.Model() as model1:
    bounded_laplace = pm.Bound(pm.Laplace, lower=0.0)
    alpha = bounded_laplace('alpha', mu = 51704326, b = 767593)

    beta_param = SSBeta(name = 'beta_testing',
                    alpha = 1.578, 
                    beta = 0.855, 
                    loc = 0.0004,
                    scale = 7.338e-5)
    pred = pm.Gamma("test_var", alpha = alpha, beta = beta_param, 
    observed=np.array([1.72e6, 1.58e6, 1.64e6, 1.59e6]))

    trace = pm.sample(400)
1 Like
4

Sorry, I thought that you had an original model that threw the initial point failure error and then tried to bound the beta distribution and got “multiple values for keyword argument ‘transform” error. Is that code the first or the second?

5

Just try them both. The code in the original message with the bounds gives the TypeError, the second one is the original one I had tried that gives the SamplingError.

6

And the solution describe here didn’t work? Or was it unsatisfactory for some other reason? It seems far simpler.

7

That’s exactly the solution that I’m using - I had included the link in my post. Try the reproducible code and you will see the issue.

8

I tried the code and can reproduce the error. But the solution at that link shifts the beta-distributed parameter using a pm.Deterministic. Did you try that?

9

Tried it and it doesn’t error out. I assume that mu is the loc parameter.

However, I am skeptical about the proposed formula:

scale * b + mu

It doesn’t seem right. Also, sampling from the prior predictive gives outlandish results, so something is off.

10

Many of the scipy specifications are a bit quirky.

Here is what I get in scipy:

loc, scale = 10,4
a = b = 1
x = np.linspace(scipy.stats.beta.ppf(0.01, a, b, loc=loc, scale=scale),
                scipy.stats.beta.ppf(0.99, a, b, loc=loc, scale=scale), 100)
rv = scipy.stats.beta(a, b, loc=loc, scale=scale)
fig, ax = plt.subplots(1, 1)
ax.plot(x, rv.pdf(x), 'k-', lw=2)
plt.show()

Note the x-axis.

scipy

And here are some prior predictive samples in pymc:

with pm.Model():
    loc, scale = 10,4
    b = pm.Beta('b', alpha=1, beta=1)
    shifted_b = pm.Deterministic('shifted_b', (scale * b) + loc)
    idata = pm.sample(10000, return_inferencedata=True)
    az.plot_trace(idata)
    plt.show()

b (untransformed) is on top, shifted_b is on the bottom.

pymc
pymc1920×974 144 KB

1 Like