Getting beta distribution subclass with loc and scale parameters to work

Questions version agnostic
1

I’m trying to get a modification of a beta distribution that includes the loc and scale parameters to work.
The subclass is described here.

As a test, I’m just specifying weak priors.

import numpy as np
import pymc3 as pm

data = np.array([1.726567e+06, 1.589836e+06, 1.643981e+06, 1.584314e+06])
data_norm = data = (data - 1.5e6) / (2e6 - 1.5e6)

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.interval(loc, loc + scale)
        #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, point=None, size=None):
        return super().random()*self.scale + self.loc

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

# Specifying weak priors first to test SSBeta
with pm.Model() as model1:
    
    a = pm.Uniform("alpha", 0,1)
    b = pm.Uniform("beta", 0,1)
    loc = pm.Uniform("loc", 0,1)
    scale =  pm.Uniform("scale", 0,1)

    pred = SSBeta("pred", alpha = a, beta = b, loc = loc, scale = scale, observed = data)
    #pred = SSBeta("pred", alpha = a, beta = b, loc = loc, scale = scale, observed = data_norm)

    trace = pm.sample(400)

Whether I try it with the raw or normalized data, I get the same error

SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'alpha_interval__': array(0.), 'beta_interval__': array(0.), 'loc_interval__': array(0.), 'scale_interval__': array(0.)}

Initial evaluation results:
alpha_interval__   -1.39
beta_interval__    -1.39
loc_interval__     -1.39
scale_interval__   -1.39
pred                -inf
Name: Log-probability of test_point, dtype: float64

I’ve also tried it with (1,2) as the Uniform distribution parameters.

Note that I have also tried to just normalize my prior data and then model it with floc=0 and fscale=1, but I’m seeing some issues with that too, so wanted to try this approach.

2

I am not sure about the specific error message, but the priors over alpha and beta are pretty extreme. They imply a beta distribution living somewhere between a uniform (beta(1,1)) and a Haldane (beta(0,0)). The prior over loc is also a bit odd in that you expect to shift your beta distribution between 0 and 1 units in the positive direction, but this means that an unshifted beta (loc=0) is nearly impossible. Trying to use the raw data with these priors would make these problems even worse.

This seems to work for me:

with pm.Model() as model:

    a = pm.Uniform("alpha", 1, 2)
    b = pm.Uniform("beta", 1, 2)
    loc = pm.Uniform("loc", -1, 1)
    scale =  pm.Uniform("scale", 0, 2)

    pred = SSBeta("pred",
                  alpha = a,
                  beta = b,
                  loc = loc,
                  scale = scale,
                  observed = data_norm)

    trace = pm.sample(init='adapt_diag')
3

Thanks. Very strange - I launched exactly what you have and mine is just stuck with the progress bar at 0%. It has been running for 17 minutes.

4

Very strange. What version of pymc are you using?

[Edit:] Just for comparison:

Sampling 2 chains for 1_000 tune and 400 draw iterations (2_000 + 800 draws total) took 3 seconds.s]
5

3.11.5