error estimate in Metropolis-Hastings integration

Question

I am evaluating an integral of the form $$I=\frac{1}{c}\int f(x) \omega(x)dx$$ with normalization $c=\int \omega(x) dx$. Using the Metropolis-Hastings algorithm I generate points $x_1,\dots,x_n$ distributed proportional to $\omega$, and the estimate for the integral is then simply $$I =\langle f \rangle_\omega \approx \frac{1}{n}\sum_{i=1}^n f(x_n)$$ where $\langle \cdot \rangle_\omega$ is the expected value w.r.t. $\omega$. But how do I get an error estimate on this value? The obvious idea is using the sample variance $$(E_I)^2 \approx \frac{Var(f)_\omega}{n}\approx\frac{\langle f^2\rangle-\langle f\rangle^2}{n}$$ But this seems to give wrong results (the computed errors are too small). I think the problem is due to the auto-correlation between the $x_i$ coming from the Markov process. Any idea how to get a beter error estimate?