search

helske / KFAS

KFAS: R Package for Exponential Family State Space Models
64 stars 17 forks source link

Returning Importance Samples in predict.SSModel #72

Closed bking124 closed 1 year ago

bking124 commented 2 years ago

When working with non-Gaussian models such as Poisson, I would like to access draws from the forecast distribution. It seems that such a sample is being drawn in the background by predict.SSModel when nsim is set to a positive number. However, the function only returns the point predictions and errors. Is there any way to access those draws? Would I need a separate call to importanceSSM?

helske commented 2 years ago

You're right, the predict method only returns the summaries computed from the importance sampling, so you would need to call importanceSSM directly for obtaining the (weighted) samples.

I think the reason why I didn't provide an option to return the actual samples was to keep the function consistent between nsim>0, nsim=0 and linear-gaussian models, but of course, another predict function for this could be useful (which would always sample, also in case of linear-gaussian models).

bking124 commented 2 years ago

Appreciate the response and explanation! Another predict method that always samples would definitely be useful.

For now, I think I've figured out how to appropriately augment the data so that calling importanceSSM will output forecasts, but I'm not 100% sure how to combine the samples appropriately using the weights. Is it just a matter of multiplying the appropriate matrices from the function output? Do I need to scale the weights matrix so that it sums to 1 first?

helske commented 2 years ago

For computing summary statistics other than ones provided by predict you can compute weighted averages of arbitrary functions, see e.g. Section 3 of the KFAS paper. If you scale the weights first so that they sum to 1 then the weighted average turns to the weighted sum.

For doing something with individual draws (e.g. visualizations), probably the easiest way is to resample (with replacement) the draws using the importance weights, so you end up with a non-weighted sample which is easier to handle.

Btw, as importance sampling doesn't scale well with respect to the number of time points, you could also use our bssm package for particle filtering based sampling of the predictions after converting the model to proper form with as_bssm function.