n <- 10
sigma_sq <- 4
n_reps <- 10000
n_draws <- n_reps * n
sigma <- sqrt(sigma_sq)
set.seed(571)
draws <- rnorm(n_draws, mean = 0, sd = sigma)
# Each row is one replication's data set.
x <- matrix(draws, nrow = n_reps)
x_bar <- rowMeans(x)
deviations <- x - x_bar
squared_deviations <- deviations^2
var_hat <- rowMeans(squared_deviations)
bias_hat <- mean(var_hat) - sigma_sq
mc_se <- sd(var_hat) / sqrt(n_reps)
c(
simulation_estimate = bias_hat,
mc_standard_error = mc_se,
known_bias = -sigma_sq / n
)
#> simulation_estimate mc_standard_error known_bias
#> -0.3915413 0.0169017 -0.4000000Simulation Studies
Last modified: 2026-11-07 21:34:29 (PST)
Example 7 (Simulating the bias of the divide-by-\(n\) variance estimator) For the data-generating process with \(\mu= 0\), \(\sigma^2= 4\) and \(n = 10\), we simulate the bias of the divide-by-\(n\) estimator \(\hat{\sigma}^2\) of \(\sigma^2\) from the example of a biased variance estimator. That example shows that the bias is exactly \(-\sigma^2/ n = -0.4\), so the simulation can be checked against the known answer.
Each replication draws \(n = 10\) values and computes one estimate of \(\sigma^2\). We run \(M = 10000\) replications, set the seed once at the start, and report the simulation estimate of the bias with its Monte Carlo standard error. The choice \(M = 10000\) gives a Monte Carlo standard error of about \(0.02\), small next to the bias.
import numpy as np
n = 10
sigma_sq = 4
n_reps = 10000
rng = np.random.default_rng(571)
# Each row is one replication's data set.
x = rng.normal(loc=0, scale=np.sqrt(sigma_sq), size=(n_reps, n))
x_bar = x.mean(axis=1, keepdims=True)
var_hat = ((x - x_bar) ** 2).mean(axis=1)
bias_hat = var_hat.mean() - sigma_sq
mc_se = var_hat.std(ddof=1) / np.sqrt(n_reps)
print(f"simulation estimate: {bias_hat:.3f}")
#> simulation estimate: -0.390
print(f"mc standard error: {mc_se:.3f}")
#> mc standard error: 0.017
print(f"known bias: {-sigma_sq / n:.3f}")
#> known bias: -0.400The R simulation estimate of the bias is -0.392, with a Monte Carlo standard error of 0.017. It is within 0.5 Monte Carlo standard errors of the known bias -0.4. The Python run uses a different random number generator, so its estimate differs slightly, but it too agrees with -0.4 up to its Monte Carlo standard error.
Figure 1 shows the sampling distribution of \(\hat{\sigma}^2\) that the replications trace out. Its histogram uses 40 equal-width bins, which put about 250 of the 10^{4} estimates in the average bin. Fewer bins, such as 10, would hide the right skew of the distribution in a few wide bars. More bins, such as 200, would leave about 50 estimates per bin on average, and random bar-to-bar variation would hide the shape. The bin count changes only the picture, not any number reported above.
[Code]
estimates <- tibble::tibble(var_hat = var_hat)
mean_var_hat <- mean(var_hat)
line_values <- c(sigma_sq, mean_var_hat)
line_labels <- c("True variance", "Average estimate")
means <- tibble::tibble(value = line_values, line = line_labels)
estimates |>
ggplot() +
aes(x = var_hat) +
# check-unexplained-choices: allow
geom_histogram(bins = 40, fill = "grey80", color = "grey40") +
geom_vline(
data = means,
mapping = aes(xintercept = value, linetype = line)
) +
labs(x = "Estimate of the variance", y = "Replications", linetype = NULL)
The distribution is right-skewed, and its average sits below the true variance, which is the bias.