Let $\beta$ be parameters in a regression ($\mu = X\beta$), $\hat{\beta}$ the estimated parameters, and $V$ the estimated covariance matrix.

I want to simulate from the posterior distribution (assumed to be multivariate normal).

What is the best way to obtain (very approximate, since the posterior distribution is approximate) confidence (?credible?) intervals?

Currently, I just simulate from the posterior a million times, and then I take the 2.5 and 97.5 quantiles on either side.


A million samples sounds like overkill. If you think about computing the quantile of a particular value (rather than the value of a particular quantile, i.e. estimating that Q(100.7) = 0.975 rather than V(0.975)=100.7), this will be determined by binomial errors. If you have 100 samples and try to get the 97.5th quantile, that will fall between observations 97 and 98 (so, not very precise): the standard errors will be on the order of $\sqrt{p(1-p)/n} \approx 0.016$. Every time you increase the sample size by 100 you'll get 10-fold higher precision. If I take a sample of 1000 and I compute the 97.5th percentile I'd expect it to be good to about $\pm 1$ percentile (2*SE = 0.0099).

My rule of thumb for interval estimates of a 1-dimensional summary statistic would be 1000 samples (maybe 10,000 if I were being obsessive), more if I need a multi-dimensional summary.

Unfortunately I don't have a reference to offer, perhaps someone else will.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.