tom-hc-park / MSc-RA-Bayesian-evidence-synthesis

Research Project at M.Sc. Statistics at University of British Columbia
4 stars 0 forks source link

a recommended package for posterior predictive samples #7

Closed tom-hc-park closed 4 years ago

tom-hc-park commented 4 years ago

Dear Prof. Paul,

I hope you are doing great.

I would like to ask your opinion/experience about how to get posterior predictive samples in JAGS.

I am worrying about which package should I use for getting posterior predictive samples (U_t, X_a). I've been only using 'rjags' and 'coda' packages and I could not find a function to generate posterior predictive samples from both of the packages. So I googled how to do this in JAGS, and found the following two packages mentioning about the functionality:

  1. jagsUI https://cran.r-project.org/web/packages/jagsUI/index.html
  2. morse https://cran.csiro.au/web/packages/morse/morse.pdf

So I guess I should use either of the two packages to get the posterior predictive samples... ? Or I would like to know how did you obtain posterior predictive samples in general..

Thanks, Tom

paulgstf commented 4 years ago

Good question @tom-hc-park . I don't have any experience with either package you mention. But a different option that might be attractive is to extract the MCMC draws from say JAGS to R and then write your own script to sample from the posterior predictive.

Small/simple example. Say the posterior distribution of the q unknown parameters theta is numerically represented by an m by q matrix, theta.sim, i.e., have m MCMC draws from the posterior.

And say the distribution for a scalar Y_new given theta is Poisson( h(theta)).

Then

rpois(m, apply(theta, 1, h))

would give m draws from the predictive dist of Y_new | Y_obs (because literally each of the m elements arises from first drawing a theta from the posterior and then drawing a Y_new given that theta).

Hope this helps.

tom-hc-park commented 4 years ago

Thank you for your reply @paulgstf. I think I understand the idea!

From the example, I guess you meant rpois(m, apply(theta.sim , 1, h)) instead of rpois(m, apply(theta, 1, h)).

I will implement this to the Rmarkdown file.

Have a great day.

paulgstf commented 4 years ago

Yes, that's right.