set.seed(715)  # retain for exactly same simulation, delete for fresh run
m = 10^6; n = 5; mu = 0; sg = 4
x = rnorm(m*n, mu, sg);  MAT = matrix(x, nrow = m)
v = rowMeans((MAT - mu)^2)  # using 'known' population mean
mean(v);  mean((v-sg^2)^2)
[1] 15.99998    # aprx E(V) = 16
[1] 102.5       # aprs MSE(V) = 102.4

hist(5*v/sg^2, prob=T, br=40, xlab="q", col="skyblue2", main="")
  curve(dchisq(x, 5), add=T, lwd=2, n=1001)