# How to estimate variance in Bayesian matrix factorization using Gibbs samples?

I have implemented a Gibbs sampler for Bayesian Matrix Factorization /Completion of matrix $$R = (r_{ij})$$ which is $$(N, M)$$ dimensional and $$p(r_{ij} | \textbf{u}_i, \textbf{v}_j) = N(r_{ij}|\textbf{u}_i^T \textbf{v}_j, \beta^{-1})$$ or simply $$r_{ij} = \textbf{u}_i^T \textbf{v}_j + \epsilon_{ij}$$ where $$\epsilon_{ij} \sim N(\epsilon_{ij}|0, \beta^{-1})$$. (Exactly like linear regression)

I calculated the two conditional posteriors and got the 'S' samples as $$\{U^s, V^s\}_{s=1}^S$$ where $$U^{s} =\{ \textbf{u}_i^{\;s} \}_{i=1}^N \text{ and } V^{s} =\{ \textbf{v}_j^{\;s} \}_{j=1}^M$$

Now, I want to approximate the mean and variance of $$r_{ij}$$ using these samples.

In my notes, while discussing MCMC, we only discussed how we can do Monte Carlo averaging and I was able to calculate

$$E(r_{ij}) = \frac1{S} \sum_{i=1}^S(\textbf{u}_i^{\;s})^T \textbf{v}_j^{\;s}$$

However, I don't have any idea how we can approximate the variance of $$r_{ij}$$ using these samples.