MCMC convergence, analytic derivations, Monte Carlo error

Question

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.

Answer 1

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.