Thanks for the response @ricardoV94. I was wondering if I could do something like this.
But then I realized that if I have 1000 posterior samples, I would end up running 4 MCMC chains with say, 500 tunes and 500 draws at least, for every 1000 posterior samples, and would have to verify the convergence for each 1000 predicted datasets. This would be a lot of computational effort I guess. But if you suggest that this option is technically correct, then I would certainly try that.
Or maybe I would just need 1 draw after maybe 500 tune samples. But then how would I know that my chains converged to a stationary distribution and the “1” sample is from the posterior predictive distribution? I don’t know what would be technically correct. I will think more about this, thanks for the suggestion.