Posterior mean with MCMC

Let's say we have a posterior distribution: $$\pi(\theta_1, \theta_2, \theta_3 | \bf{y})$$

and that we've run an MCMC algorithm to approximate this distribution.

I know that there is a Markov chain version of the CLT that states: $$\sqrt{n}\left(\frac{1}{n} \sum\limits_{i=1}^{n} X_i - E(X) \right) \xrightarrow[]{d} \mathcal{N}(0,\sigma^2)$$ where $$n$$ is the sample size, and $$X \sim \pi$$.

Normally with a joint pdf, you could integrate out the other parameters and then take the expectation to get the posterior means, but what does that look like with my MCMC samples? Does it follow from this CLT? Is it also true that a good posterior mean estimate for $$E[\theta_j]$$ is just $$\frac{1}{n} \sum\limits_{i=1}^{n} \theta_{ij}$$ (where $$\theta_{ij}$$ is the $$i$$th sample of the $$j$$th parameter)?

$$\pi(\theta_1,\theta_2,\theta_3 \lvert y)$$
Because it is a random sample $$\{ (\theta_{i1},\theta_{i2},\theta_{i3})\}_{i=1}^n$$ anything that holds for random samples apply. This means that Law of Large Numbers will apply such that a consistent estimator of $$\mathbb E[\theta_j ]$$ is the average $$(1/n) \sum_i \theta_{ij}$$. The CLT can be used to justify asymptotic distribution and make tests if you should desire. And you can approximate integrals, for example $$\mathbb E[h(\theta_2)]$$ can be approximated
$$\mathbb E[h(\theta_2)] = \int_{\theta_1} \int_{\theta_2}\int_{\theta_3} h(\theta_2) \pi(\theta_1,\theta_2,\theta_3 \lvert y)\approx \frac{1}{n} \sum_i h(\theta_{i2})$$ which is just Law of Large Numbers again.