set.seed(2019)
m = 10^6;  n = 5;  mu = 50;  sg = 7
x = rnorm(m*n, mu, sg)
MAT = matrix(x, nrow=m)      # each row a sample of             n
v1 = apply(MAT, 1, var)      # uses sample mean
v2 = rowSums((MAT - mu)^2)/n # uses population mean
q1 = (n-1)*v1 /sg^2
q2 = n*v2 / sg^2

par(mfrow=c(1,2))
hist(q1, prob=T, br=50, col="skyblue2", ylim=c(0,.2))
 curve(dchisq(x, n-1), add=T, lwd=2)
 curve(dchisq(x, n), add=T, col="red", lwd=2, lty="dotted")
hist(q2, prob=T, br=50, col="skyblue2")
  curve(dchisq(x, n), add=T, lwd=2)
par(mfrow=c(1,1))