Calculating the cumulative sum of a modeled random count variable

Hello,

It seems like you might be able to just use the .cumsum() method, together with the posterior samples you have already sampled in post, to save yourself a lot of work. For example:

(post.posterior_predictive['count']
             .cumsum(dim='count_dim_0')
             .stack(sample = ['chain', 'draw'])
             .quantile([0.05, 0.5, 0.95], dim='sample')
             .plot
             .line(x='count_dim_0'));

Here, I picked out the predicted counts from the posterior predictive, took the cumulative sum over the “count_dim_0” dimension (which is the time dimension, it’s useful to rename these using coords!), stacked all the samples from all chains into a single dimension called “sample”, computed quantiles, then plotted. Here are the results:

I don’t think it necessary to look at every possible combination of trajectories as you do, because your samples are all independent in time anyway. If you need more trajectories, rather than building them by permuting existing trajectories, it is easier to just draw more samples from posterior. In this case you can do either, but be aware that if you were to add some kind of time-series dynamics to your model, it would not be legit to jump across trajectories like this.

You also don’t need to manually build the posterior by sampling mu from a kde then sampling a poisson, you already have all this using post = pm.sample_posterior_predictive(trace).