# Law of Total Variance

I trying to experiment with law of total variance in order to empirically recreate theoretical results.
In particular I am interested in verifying that:
$$Var(Y) = E[Var(Y|X)] + Var(E[Y|X])$$ Let's assume I have the following random variables:
$$X \sim Gamma(\alpha, \beta)$$

$$Y \sim Poisson(X)$$

The total variance of Y should be equal to:
$$Var(Y) = E[Var(Y|X)] + Var(E[Y|X]) = E[X] + Var(X) = \alpha*\beta + \alpha*\beta^2$$ This follow from $$E[X] = \alpha*\beta$$, $$Var(X) = \alpha*\beta^2$$, $$E[Y|X] = Var(Y|X) = X$$, which are known results for the Gamma and Poisson distribution.
I tried to recreate such results in R with the following code:

alpha <- 250
beta <- 10

n <- 10000

X <-  rgamma(n, shape = alpha, scale = beta)

Y <- replicate(n, rpois(1, sample(X, 1)))

alpha * beta + alpha * beta^2

[1] 27500

mean(X) + var(X)

[1] 27552.47

var(Y)

[1] 26610.7


In this situation the total variance is 27500 and when I calculate it manually with the previous results it is, in fact, 27500 (or close enough).
My question is, why when I do var(Y) I get a wrong result? Am I defining the variable Y in the right way?

Thank you

Your estimator is quite close to the true variance, but it appears you might need a higher sample size to get as close as you want. (You also forgot to set the seed to get a reproducible analysis.) I'm going to try it with a larger number of samples:

#Set parameters
alpha <- 250
beta <- 10

#Set sample size
n <- 10^7

#Generate random variables
set.seed(1)
X <- rep(0, n)
Y <- rep(0, n)
for (i in 1:n) {
X[i] <- rgamma(1, shape = alpha, scale = beta)
Y[i] <-  rpois(1, lambda = X[i]) }

#Check sample variance
var(Y)
[1] 27522.87

• It would be interesting to see how to decompose the total variance in case I have a distribution that depends on two parameters which also are stochastic in nature. Commented Feb 28, 2022 at 22:39