## Estimate actual coverage probabilities of nominal 95% Bonferroni
## simultaneous CI's for all coefficients in a regression fit (except
## the intercept).

## Sample size.
n = 100

## Rank of the design matrix.
p = 4

## Number of simulation replications.
nrep = 1e3

## Construct a design matrix.
r = 0.2 ## Must be non-negative.
v = sqrt(r/(1-r))
X = array(rnorm(n*p), c(n,p))
X = (X + v*outer(rnorm(n), array(1,p)))/sqrt(1+v^2)

## Count the Bonferroni interval sets that simultaneously cover.
ncb = 0

## Simulation replications.
for (r in 1:nrep)
{
  ## True regression coefficients (population intercept is zero).
  beta = rnorm(p)

  ## Simulate the response data.
  Y = X %*% beta + rnorm(n)

  ## Fit the model.
  mm = lm(Y ~ X)
  sm = summary(mm)

  ## Residuals and error variance estimate.
  sigma2_hat = sm$sigma^2

  ## (X' * X)^-1
  MI = sm$cov.unscaled

  ## Bonferroni t-quantile for 95% nominal coverage.
  qb = qt(1-0.05/(2*p), n-p-1)

  ## Slope estimates (dropping the intercept which we aren't analyzing).
  beta_hat = mm$coeff[2:(p+1)]

  ## Check whether the Bonferroni intervals for all coefficients cover.

  ## The CI half widths.
  c = qb*sqrt(sigma2_hat*diag(MI)[2:(p+1)])

  ## True it at least one element does not cover.
  ncb = ncb + 1*(all(beta > beta_hat-c) & all(beta < beta_hat+c))
}

## The estimated coverage probability.
bcp = ncb / nrep