Generating Posterior Predictive Samples Over Forecast Period

Questions
1

I have a simple linear model here that takes time as an argument in years from 2012. It also takes a flag argument that allows for the modification of slope and intercept in the first year of the time series. My data has a first-year effect that i’m trying to account for with these terms, and so have modeled it in an interrupted times series approach. It’s also hierarchical with 36 different categories.

What I would like to do is generate posterior predictive samples over the future, calculated the hpd over these samples, and plot them.

I’ve used plot_posterior_predictive_glm() which allows me only to specify a very limiting, single variable function. It doesn’t allow me to include the terms that account for the first year effects. Any guidance on how I can generate posterior predictive samples over a forecast period?

with pm.Model() as boe_multilevel:
    mu_b1 = pm.TruncatedNormal('mu_b1', mu=-0.09, sigma = 1, upper = -0.03)
    a = pm.Normal('a', mu = 8, sigma = 10, shape=cat_len)
    b1 = pm.TruncatedNormal('b1', mu = mu_b1, sigma = 1, upper=0.0, shape = cat_len)
    a2 = pm.Normal('a2', mu = 0, sigma = 10, shape = cat_len)
    b2 = pm.Normal('b2', mu = 0, sigma = 10, shape= cat_len)
    eps = pm.HalfCauchy('eps',0.5, shape=cat_len)
    sigma_y = eps[cat_codes]
    nu = pm.DiscreteUniform('nu',1,10)

    
    boe_est = a[cat_codes] + b1[cat_codes]*train1['time'].values + a2[cat_codes]*train1['capital_year_flag'] + b2[cat_codes]*train1['time']*train1['capital_year_flag'] 
    
    y = pm.StudentT('y', mu = boe_est,nu = nu , sigma = sigma_y, observed=train1['boe'].apply(np.log))

with boe_multilevel:
    boe_trace = pm.sample(4000,cores=4)

ppc_y = pm.sample_posterior_predictive(trace=boe_trace, model=boe_multilevel)
ppc_coefs = pm.sample_posterior_predictive(trace=boe_trace, model=boe_multilevel, var_names = ['a','a2','b1','b2'])

for i in range(36):
    lm = lambda x, sample: np.exp(sample['a'][i] + sample['b1'][i]*x)
    pm.plot_posterior_predictive_glm(boe_trace,lm =lm, samples=200,eval=np.linspace(7,9,36),alpha=0.3)```
2

ppc_y, in the code above, takes the shape of (#samples) x (#observations in y). Is there a way to get samples for any array of inputs of arbitrary length? Similar to how plot_posterior_predictive_glm() does so for a the vector “eval=”

3

Hi Justin,
I think you need to cast your time variable – or any variable over which you want posterior predictive samples – as a pm.Data variable. That way you’ll be able to resample from the model in no time, by just telling PyMC to change the predictors’ values.

Here is a notebook that explains how to do it, and I recently wrote a detailed NB on how to do it for a multinomial regression.

Hope this helps :vulcan_salute:

2 Likes
4

This is super helpful. I think i’m now on the right track, but I’ve somehow broken my code causing my python kernel to crash.

with pm.Model() as boe_multilevel:
    data_time = pm.Data('data_time',train1['time'].values)
    data_boe = pm.Data('data_boe',train1['boe'].values)
    data_codes = pm.Data('data_codes',train1['cat_codes'].values)
    data_wedgeflg = pm.Data('data_wedgeflg', train1['capital_year_flag'].values)
   
    mu_b1 = pm.TruncatedNormal('mu_b1', mu=-0.09, sigma = 1, upper = -0.03)
  
    a = pm.Normal('a', mu = 8, sigma = 10, shape=cat_len)
    b1 = pm.TruncatedNormal('b1', mu = mu_b1, sigma = 1, upper=0.0, shape = cat_len)
    a2 = pm.Normal('a2', mu = 0, sigma = 10, shape = cat_len)
    b2 = pm.Normal('b2', mu = 0, sigma = 10, shape= cat_len)
    eps = pm.HalfCauchy('eps',0.5, shape=cat_len)
    sigma_y = eps[cat_codes]
    nu = pm.DiscreteUniform('nu',1,10,shape=cat_len)[cat_codes]

    
    boe_est = a[cat_codes] + b1[cat_codes]*data_time + a2[cat_codes]*data_wedgeflg + b2[cat_codes]*data_time*data_wedgeflg 
    
    y = pm.StudentT('y', mu = boe_est,nu = nu , sigma = sigma_y, observed=np.log(data_boe))

So, whereas before I was using a pandas dataframe as the data within the model. I have separated them out into individual SharedVariables using pm.Data(). This, however, causes my kernel to crash in jupyter lab. I’m not detecting my mistake. If I can find it, perhaps with your help, my next move will be to update the sharedvariables with pm.set_values() per your notebook examples (super helpful thank you!) before running pm.sample_posterior_predictive()

Does it sound like i’ve got it?

1 Like
5

Hi Justin,
Glad to know it was helpful! Yeah, sounds like a good roadmap to me :slight_smile:
What do you mean by “causes my kernel to crash” though?

6

I call sample and Jupyter Lab crashes as seen below.

Capture2
Capture21309×248 22.3 KB

Capture
Capture1570×890 61.3 KB

7

Well that’s weird! It looks like a Theano-related issue, doesn’t it?

8

Okay. The issue is definitely Theano on my end. I shut everything down. Deleted Theano folder from appdata/local/ and it actually finished sampling. I reran it and it crashed again… But at least I know it’s not the code.

9

What if you do cores=1?

10

When I set cores=1 I don’t get a kernel crash, but instead get “Sampling Error: Bad Initial Energy”

Edit: Neverymind, I figured this error out. I was taking np.log() over data with 0 values causing -infs being passed to observed. I added a small epsilon, and everything is running now.

What do you think is the deal with the multiple cores causing the Theano issues? I read maybe I need to increase the compile time?

11

In this case I suspect it is a memory problem.

12

Ok issues persist.

I’m now trying to set the values for the data I set up before:

with boe_multilevel:
    pm.set_data({'data_time': np.linspace(7,10,36),'data_wedgeflg': np.zeros(36)})
    ppc_y = pm.sample_posterior_predictive(trace=boe_trace)

and I’m getting an dimension mismatch error from theano. I suspect it has to do with the fact that I have multiple time series (36) and each one is given a category code (cat_codes) that I use to specify coefficients for each. I can’t set the cat_codes within the model because it’s passed into my model as an array and not declared under the pm.Data() method. However, when I try to use the pm.Data() method on the category codes and then use the category codes (as below) I get an error because theano is converting the integer codes into floats, and I get an error for passing floats into the indexer.

data_codes = pm.Data('data_codes',train1['cat_codes'].values)

boe_est = a[data_codes] + b1[data_codes]*data_time + a2[data_codes]*data_wedgeflg + b2[data_codes]*data_time*data_wedgeflg

So how should I be handling the category codes? How does posterior_predictive_sample() work when the sample is dependent on a category code to specify which coefficients are needed to generate the posterior pred for y?

EDIT:
Calling pm.intX(pm.Data(‘data_codes’,train1[‘cat_codes’].values)) stopped the error I was getting from trying to pass floats into the indexer. So this code is NOT throwing errors when defining the model.

data_codes = pm.intX(pm.Data('data_codes',train1['cat_codes'].values))

boe_est = a[data_codes] + b1[data_codes]*data_time + a2[data_codes]*data_wedgeflg + b2[data_codes]*data_time*data_wedgeflg
1 Like
13

Just to clarify. Here is the code all together. With the errors:

with pm.Model() as boe_multilevel:
    data_time = pm.Data('data_time',train1['time'].values)
    data_oil = pm.Data('data_oil',np.log(train1['oil_prod'].values))
    data_codes = pm.intX(pm.Data('data_codes',train1['cat_codes']))
    data_wedgeflg = pm.Data('data_wedgeflg', train1['capital_year_flag'].values)
    #data = pm.Data('data', train1[['boe','time','cat_codes','capital_year_flag']].values)
    
    # Hyperpriors
    mu_b1 = pm.TruncatedNormal('mu_b1', mu=-0.09, sigma = 1, upper = -0.03)

    # Priors
    a = pm.Normal('a', mu = 0, sigma = 10, shape=cat_len)
    b1 = pm.TruncatedNormal('b1', mu = mu_b1, sigma = 1, upper=0.0, shape = cat_len)
    a2 = pm.Normal('a2', mu = 0, sigma = 10, shape = cat_len)
    b2 = pm.Normal('b2', mu = 0, sigma = 10, shape= cat_len)
    eps = pm.HalfCauchy('eps',0.5, shape=cat_len)
    sigma_y = eps[cat_codes]
    #nu = pm.DiscreteUniform('nu',1,10,shape=cat_len)[cat_codes]

    
    oil_est = a[data_codes] + b1[data_codes]*data_time + a2[data_codes]*data_wedgeflg + b2[data_codes]*data_time*data_wedgeflg 
    
    y = pm.StudentT('y', mu = oil_est, nu = 5 , sigma = sigma_y, observed=data_oil)

with boe_multilevel:
    boe_trace = pm.sample(4000,cores=2)

with boe_multilevel:
    pm.set_data({'data_time': np.linspace(7,10,36),'data_wedgeflg': np.zeros(36),'data_codes': np.ones(36)})
    ppc_y = pm.sample_posterior_predictive(trace=boe_trace)

ValueError                                Traceback (most recent call last)
~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\shape_utils.py in broadcast_dist_samples_shape(shapes, size)
    165     try:
--> 166         broadcast_shape = shapes_broadcasting(*sp_shapes, raise_exception=True)
    167     except ValueError:

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\shape_utils.py in shapes_broadcasting(raise_exception, *args)
     91                     "Supplied shapes {} do not broadcast together".format(
---> 92                         ", ".join(["{}".format(a) for a in args])
     93                     )

ValueError: Supplied shapes (), (36,), (952,), (952,) do not broadcast together

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py in _draw_value(param, point, givens, size)
    575                 try:
--> 576                     return dist_tmp.random(point=point, size=size)
    577                 except (ValueError, TypeError):

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\continuous.py in random(self, point, size)
   1963                                 dist_shape=self.shape,
-> 1964                                 size=size)
   1965 

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py in generate_samples(generator, *args, **kwargs)
    696             [np.asarray(i).shape for i in inputs] + [_dist_shape],
--> 697             size=size_tup
    698         )

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\shape_utils.py in broadcast_dist_samples_shape(shapes, size)
    169             "Cannot broadcast provided shapes {} given size: {}".format(
--> 170                 ", ".join(["{}".format(s) for s in shapes]), size
    171             )

ValueError: Cannot broadcast provided shapes (), (36,), (952,), (952,) given size: ()

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\shape_utils.py in broadcast_dist_samples_shape(shapes, size)
    165     try:
--> 166         broadcast_shape = shapes_broadcasting(*sp_shapes, raise_exception=True)
    167     except ValueError:

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\shape_utils.py in shapes_broadcasting(raise_exception, *args)
     91                     "Supplied shapes {} do not broadcast together".format(
---> 92                         ", ".join(["{}".format(a) for a in args])
     93                     )

ValueError: Supplied shapes (), (36,), (952,), () do not broadcast together

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-241-4ff2bd752351> in <module>
      1 with boe_multilevel:
      2     pm.set_data({'data_time': np.linspace(7,10,36),'data_wedgeflg': np.zeros(36),'data_codes': np.ones(36)})
----> 3     ppc_y = pm.sample_posterior_predictive(trace=boe_trace)
      4 #ppc_coefs = pm.sample_posterior_predictive(trace=boe_trace, model=boe_multilevel, var_names = ['a','a2','b1','b2'])

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\sampling.py in sample_posterior_predictive(trace, samples, model, vars, var_names, size, random_seed, progressbar)
   1112                 param = trace[idx % len_trace]
   1113 
-> 1114             values = draw_values(vars, point=param, size=size)
   1115             for k, v in zip(vars, values):
   1116                 ppc_trace[k.name].append(v)

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py in draw_values(params, point, size)
    393                                         point=point,
    394                                         givens=temp_givens,
--> 395                                         size=size)
    396                     givens[next_.name] = (next_, value)
    397                     drawn[(next_, size)] = value

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py in _draw_value(param, point, givens, size)
    583                     with _DrawValuesContextBlocker():
    584                         val = np.atleast_1d(dist_tmp.random(point=point,
--> 585                                                             size=None))
    586                     # Sometimes point may change the size of val but not the
    587                     # distribution's shape

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\continuous.py in random(self, point, size)
   1962         return generate_samples(stats.t.rvs, nu, loc=mu, scale=lam**-0.5,
   1963                                 dist_shape=self.shape,
-> 1964                                 size=size)
   1965 
   1966     def logp(self, value):

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\distribution.py in generate_samples(generator, *args, **kwargs)
    695         broadcast_shape = broadcast_dist_samples_shape(
    696             [np.asarray(i).shape for i in inputs] + [_dist_shape],
--> 697             size=size_tup
    698         )
    699         # We do this instead of broadcast_distribution_samples to avoid

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pymc3\distributions\shape_utils.py in broadcast_dist_samples_shape(shapes, size)
    168         raise ValueError(
    169             "Cannot broadcast provided shapes {} given size: {}".format(
--> 170                 ", ".join(["{}".format(s) for s in shapes]), size
    171             )
    172         )

ValueError: Cannot broadcast provided shapes (), (36,), (952,), () given size: ()

EDIT: I FIGURED IT OUT!!!

sigma_y = eps[cat_codes]

should be

sigma_y = eps[data_codes]

1 Like
14

In this video, https://www.youtube.com/watch?v=rZvro4-nFIk&t=1934s, Nicole Carlson says that I need also replace the observed data to an arbitrary vector of the same length.

I’m going to try this and report back.

Edit: Didn’t work :frowning:

Edit2: I think I figured it out. See above.