DensityDist conditional random() method

Questions
1

Hi pymc Comunity,
I’d like to do posterior predictive checks on a DensityDist.
I’m implementing a random method for the custom DensityDist distribution, however, the ‘random’ values produced are conditional on other observed data. How can I get this data into the random() method? Maybe something to do with the point input?

I have a simplified version on my code below, where a value (next_epoch) is sampled categorically, conditioned on another value (current_epoch)

with pm.Model() as model:
    BoundedNormal = pm.Bound(pm.Normal, lower=0, upper=1)
    trans_baserate = BoundedNormal('trans_baserate', mu=0.5,
                         sd=1e5, shape=(nstages, nstages))

    def logp(current_epoch, next_epoch):
        current_idx = tt.cast(current_epoch, 'int16')
        next_idx = tt.cast(next_epoch, 'int16')
        cat = pm.Categorical.dist(p=trans_baserate[current_idx, :])
        return cat.logp(next_idx)


    def random(point=None, size=None):
        trans_baserate_ = pm.distributions.multivariate.draw_values([trans_baserate], point=point, size=size)
        def _random(point, trans_baserate, size=None):
            trans_p = trans_baserate[point['current_epoch'], :] #Complains here :frowning: 
            return pm.distributions.dist_math.random_choice(p=trans_p, size=size)  

        return pm.distributions.multivariate.generate_samples(_random,
                                                              point=point,
                                                              trans_baserate=trans_baserate_,
                                                              size=size)

    next_epoch = pm.DensityDist('next_epoch', logp,
                              observed={'current_epoch': data['current_epoch'],
                                                'next_epoch': data['next_epoch']},
                              random=random)

The low response rate here makes me wonder if this is even possible? Maybe shared theano variables could help?

2

Hi, can any of the devs help me here please? Or is this not possible? Conversely, I will raise this on github.

3

There is no FreeRV in your model call ‘current_epoch’ - that’s why it complains.
How about casting current_epoch to theano shared variable and just use it like in the logp: trans_p = trans_baserate[current_epoch, :] ?

4

Amazing!
Seems a little hacky using a variable from the parent scope, but it works :slight_smile:
Further, I had to call .eval() on the shared current_epoch in the random() function, as the random function should be using numpy arrays and not theano tensors.

5

I dont think it is a good practice to call .eval() - try draw_valued instead:

6

Thanks,
I tried that but:

current_epoch_values = draw_values(current_epoch_shared)

throws:


  File "/home/bdyetton/PSleep/src/modeling/sleep_stage_models.py", line 208, in _random
    current_epoch_values = draw_values(current_epoch_shared)
  File "/home/bdyetton/anaconda3/envs/psleep/lib/python3.7/site-packages/pymc3/distributions/distribution.py", line 283, in draw_values
    params = dict(enumerate(params))
  File "/home/bdyetton/anaconda3/envs/psleep/lib/python3.7/site-packages/theano/tensor/var.py", line 628, in __iter__
    for i in xrange(theano.tensor.basic.get_vector_length(self)):
  File "/home/bdyetton/anaconda3/envs/psleep/lib/python3.7/site-packages/theano/tensor/basic.py", line 4828, in get_vector_length
    raise ValueError("length not known: %s" % msg)
ValueError: length not known: <TensorType(int64, vector)> [id A]

Giving a size input to draw values does not help either.

7

I think it needs a list as input, try:

current_epoch_values = draw_values([current_epoch_shared])[0]
8

Ahhh, silly me. Thanks!
+1 for type hints in pymc4 :wink: