2
$\begingroup$

I'm trying to figure out some convergence statements on an MCMC example.

The setup is: I'm generating data samples as observations from a (known) deterministic parameter, say $s$ (using a forward equation and a random pertubation term). The data are input to my MCMC and I try to get estimates of the parameter $s$. I know it's not sensible to use MCMC to estimate a known deterministic parameter but this is just a toy example. So I assume I don't know much about $s$, give it an uninformative prior and set up a likelihood according to my problem.

Then I run my MCMC and I want to quantify how good my solutions are. So the assessements I want to make are along the lines of "I have this and that assumptions on $s$, I would expect e.g. the posterior mean be of an error magnitude of [whatever]... and my simulations meet these expectations yes or no..." How do I go about this?

Is it possible to make such statements without an analytical MCMC target/posterior? (i.e. in my example I know the "true" distribution of $s$ (deterministic in this case), but I'm not putting this into the MCMC algorithm - the posterior for $s$ I get from likelihood and prior are something I don't have an analytical expression for) So far I figured that for error assessement of the posterior mean one could examine

$| s - \int_{\mathbb{R}} s \ p(s | data) \ ds | \leq \ldots$

but I'm not getting anywhere here if I don't know the posterior p(s | data)...

I would appreciate any hints.

Edit: Thanks for the help! My original question was obviously a little vague and also contains more than one question. In think the numerical error by MCMC is answered very clearly by Greenparkers post. I have more questions on the model error of the inverse problem itself but I think that's a different issue and so I'll mark this question answered and maybe make new ones for the other stuff.

$\endgroup$
  • $\begingroup$ Why not sample from the posterior distribution using MCMC and average them? $n^{-1}\sum_{i=1}^{n} X_{i} \to \mathbb{E}\left[X\right]$ if $n \to \infty$ almost surely by strong law of large numbers. So the integral in your error assessment could be computed using $n^{-1}\sum_{i=1}^{n}X_{i}$. $\endgroup$ – Daeyoung Lim Mar 15 '16 at 11:06
  • $\begingroup$ Thanks. Your're right... I can compute the intergral this way. Then, so I figured out, there are known results on (Markov Chain) Monte Carlo standard error I can use for quantifying the error $| \mathbb{E}(X) - \frac{1}{n} \sum_{i=1}^{n} X_i$ |. So that would give the numerical error I do by MCMC. Anyway I'm still lacking something for the analysis I want - along the lines of model error rather than numerical error. Maybe I update the question in a few days... $\endgroup$ – jenna Mar 17 '16 at 14:16
2
$\begingroup$

Once you obtain samples from the posterior $X_1, X_2, \dots, X_n$ and calculate the sample mean, $$s_n = \dfrac{1}{n} \sum_{i=1}^{n} X_i \to E(s), $$

your question now is, how much error have I made in estimating $E(s)$ (assuming $E(s)$ perfectly summarizes $s$). Thus you want to know the Monte Carlo error $$s_n - E(s). $$

Since you don't know $E(s)$, you cannot find that quantity exactly. However, if a Markov chain CLT exists,

$$\sqrt{n}(s_n - E(s)) \overset{d}{\to}N(0, \sigma^2), $$

where $\sigma^2 > 0$ is the asymptotic variance of the Monte Carlo error. Thus, in order to understand the Monte Carlo error, you want to learn $\sigma^2$. $\sigma^2$ is tricky because it is not the posterior variance of $s$, and includes the correlation in the Markov chain. There are consistent estimators of $\sigma^2$ available using batch means or spectral variance methods. See

http://arxiv.org/abs/math/0703746

http://arxiv.org/abs/0811.1729v5

http://arxiv.org/abs/math/0601446

$\sigma^2$ can be estimated using the R package mcmcse. Below I provide example code for a simple AR(1) Markov chain. The mcse function finds the estimatate of $\sigma/\sqrt{n}$.

> library(mcmcse)
> n <- 10^5
> x <- numeric(length = n)
> x[1] <- 0
> 
> # Ar(1) Markov chain
> for(i in 2:n)
+ {
+   x[i] <- .4*x[i-1]+ rnorm(1)
+ }
> 
> # Batch means estimate of sigma/sqrt(n) and posterior mean
> (st_error <- mcse(x, method = "bm"))
$est
[1] 0.0001526641

$se
[1] 0.004943979

> 
> # asymptotic variance is not equal to variance of s
> sd(x)/sqrt(n)
[1] 0.00343409

Thus, in this example code you will report the Monte Carlo standard error to be .00494, and that summarizes your Monte Carlo error. That is all you can say about the Monte Carlo error.

If $s$ is of dimension larger than $1$ so that $E(s)$ is a vector, than the same idea can be repeated by looking at the multivariate Markov chain CLT $$\sqrt{n}(s_n - E(s)) \overset{d}{\to} N(0, \Sigma), $$

where now $\Sigma$ is the asymptotic covariance structure of the Monte Carlo error and can be estimated using batch means methods. See http://arxiv.org/abs/1512.07713. For implementation, you can estimate $\Sigma$ using the mcse.multi function in the same package.

$\endgroup$

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.