AssertionError with logcdf method for a Beta

Questions
1

Hi,

I’m having a AssertionError issue when using the logcdf method for a Beta, but strangely not for a Normal or a Uniform. From what I understand from other posts on AssertionErrors, this error appears when the incoming object in my deterministic distribution is of the wrong format. But I’m surprised it happens only for Betas. Here is my Code:

# Initialize random number generator
np.random.seed(123)

# True parameter values
alpha, sigma = 1, 1
beta = [1, 2.5]

# Size of dataset
size = 100

# Predictor variable
X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.2

# Simulate outcome variable
data = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma
>
basic_model = pm.Model()

with basic_model:

    a = pm.Normal('a', mu=2, sd=2)
    b = pm.Normal('b', mu=2, sd=2)
    s = pm.HalfNormal('sigma', sd=0.01)

    m1 = np.exp(pm.Beta.dist(alpha=a, beta=b).logcdf(0.1))
  
    Y_obs = pm.Normal('Y_obs1', mu=m1, sd=s, observed=data)


map_estimate = pm.find_MAP(model=basic_model)
print map_estimate

with basic_model:
    # Sample from the posterior
      trace = pm.sample(draws=1000, chains=2, tune=500, discard_tuned_samples=True)

      pm.plot_posterior(trace)

and here is my error log :

Traceback (most recent call last):
File “/Users/Colo/Google Drive/Projects/greghec/learning pymc2.py”, line 47, in
trace = pm.sample(tune=500)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/sampling.py”, line 395, in sample
progressbar=progressbar, **args)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/sampling.py”, line 1515, in init_nuts
step = pm.NUTS(potential=potential, model=model, **kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/step_methods/hmc/nuts.py”, line 154, in init
super(NUTS, self).init(vars, **kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/step_methods/hmc/base_hmc.py”, line 75, in init
vars, blocked=blocked, model=model, dtype=dtype, **theano_kwargs
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/step_methods/arraystep.py”, line 228, in init
vars, dtype=dtype, **theano_kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/model.py”, line 710, in logp_dlogp_function
return ValueGradFunction(self.logpt, grad_vars, extra_vars, **kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/model.py”, line 443, in init
grad = tt.grad(self._cost_joined, self._vars_joined)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 605, in grad
grad_dict, wrt, cost_name)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1371, in _populate_grad_dict
rval = [access_grad_cache(elem) for elem in wrt]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1162, in access_term_cache
new_output_grads)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/scan_module/scan_op.py”, line 2222, in L_op
assert (x[::-1][:-1].tag.test_value.shape[0] == n)
AssertionError

Any ideas why this happens ? Thanks !

2

I think that the problem is that the beta distribution still doesn’t implement its mode, so there is no tag test value set properly. Try giving setting tag.test_value of the beta RV to some value.

3

Another thing you could try is to replace np with theano.tensor. I’m not sure how well np.exp mixes with theano tensors.

4

I works well with floats or arrays, but if I try

with basic_model:

    # Priors for unknown model parameters
    a = pm.HalfNormal('alpha', sd=10)
    b = 1.0
    sigma = pm.HalfNormal('sigma', sd=1)
    a.tag.test_value = pm.HalfNormal('alphat', sd=10)

    m1 = pm.Beta.dist(alpha=a, beta=b).logcdf(0.1)


    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=m1, sd=sigma, observed=Y)

with basic_model:
    # Sample from the posterior
      trace = pm.sample(tune=500)
      pm.plot_posterior(trace)

I get :

Traceback (most recent call last):
File “/Users/Colo/Google Drive/Projects/greghec/learning pymc2.py”, line 48, in
m1 = pm.Beta.dist(alpha=a, beta=b).logp(0.1)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/distributions/distribution.py”, line 52, in dist
dist.init(*args, **kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/distributions/continuous.py”, line 1122, in init
self.mean = self.alpha / (self.alpha + self.beta)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/tensor/var.py”, line 128, in add
return theano.tensor.basic.add(self, other)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gof/op.py”, line 625, in call
storage_map[ins] = [self._get_test_value(ins)]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gof/op.py”, line 562, in _get_test_value
ret = v.type.filter(v.tag.test_value)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/tensor/type.py”, line 87, in filter
'Expected an array-like object, but found a Variable: ’
TypeError: For compute_test_value, one input test value does not have the requested type.

The error when converting the test value to that variable type:
Expected an array-like object, but found a Variable: maybe you are trying to call a function on a (possibly shared) variable instead of a numeric array?

5

This is a different error. You’re setting a's test value to an RV. That is not supported because test values must always be numbers or arrays.

6

I see, thanks. So with this code :

basic_model = pm.Model()

with basic_model:

    # Priors for unknown model parameters
    a = pm.HalfNormal('alpha', sd=10)
    b = 1.0
    sigma = pm.HalfNormal('sigma', sd=1)
    a.tag.test_value = 1

    m1 = pm.Beta.dist(alpha=a, beta=b).logcdf(0.1)


    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=m1, sd=sigma, observed=Y)

with basic_model:
      trace = pm.sample(tune=500)
      pm.plot_posterior(trace)

I get :

Running on PyMC3 v3.6
Auto-assigning NUTS sampler…
Initializing NUTS using jitter+adapt_diag…
Traceback (most recent call last):
File “/Users/Colo/Google Drive/Projects/greghec/learning pymc2.py”, line 62, in
trace = pm.sample(tune=500)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/sampling.py”, line 395, in sample
progressbar=progressbar, **args)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/sampling.py”, line 1515, in init_nuts
step = pm.NUTS(potential=potential, model=model, **kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/step_methods/hmc/nuts.py”, line 154, in init
super(NUTS, self).init(vars, **kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/step_methods/hmc/base_hmc.py”, line 75, in init
vars, blocked=blocked, model=model, dtype=dtype, **theano_kwargs
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/step_methods/arraystep.py”, line 228, in init
vars, dtype=dtype, **theano_kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/model.py”, line 710, in logp_dlogp_function
return ValueGradFunction(self.logpt, grad_vars, extra_vars, **kwargs)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/pymc3/model.py”, line 443, in init
grad = tt.grad(self._cost_joined, self._vars_joined)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 605, in grad
grad_dict, wrt, cost_name)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1371, in _populate_grad_dict
rval = [access_grad_cache(elem) for elem in wrt]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1021, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1326, in access_grad_cache
term = access_term_cache(node)[idx]
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/gradient.py”, line 1162, in access_term_cache
new_output_grads)
File “/Users/Colo/anaconda/lib/python2.7/site-packages/theano/scan_module/scan_op.py”, line 2222, in L_op
assert (x[::-1][:-1].tag.test_value.shape[0] == n)
AssertionError

7

I found some time to take a look at your code. It seems that the logcdf of the Beta distribution uses some scan operations to perform the computations. These seem to fail when either alpha or beta are RVs… I did not manage to find a workaround to the problem, sorry. Maybe @junpenglao can help with this.

8

hmmm, seems the gradient is broken in the scan, but I dont see an easy fix as well - doesn’t seems trivial to fix unfortunately…

9

Thanks for having tried. Any idea if I can do without using

pm.Beta.dist(alpha=a, beta=b).logcdf(0.1)

Could I implement a beta cdf by hand and plug it into a PyMC model ?

10

Yes, if you’re up for the challenge. You could implement the Op and its gradient in theano and then use that.

11

To follow up: have a look at pymc3.distributions.dist_math.incomplete_beta which calls power-series and continued-fraction expansions for the incomplete beta function (beta CDF). As these are power expansions, they use scan; and one or both of these has a broken gradient.

You might unwrap the scan to some term of the approximation that is good enough; and this should deal with the problem.

Alternatively you can implement the gradients directly rather than using automatic differentiation - as power series are easily differentiated. Continued fractions are nastier, but sum-product formulas exist for their derivatives.

2 Likes