Computational Statistical Inference: Monte Carlo Methods and the Bootstrap

Modern computing makes it possible to approximate sampling distributions, standard errors, confidence intervals, and p-values through repeated simulation.

Two foundational computational methods are:

Monte Carlo simulation
The bootstrap

Both use repeated random computation, but they generate data differently and answer different questions.

Classical Analytic Inference

For a population mean \(\mu\), the sample mean is\[ \bar{X} = \frac{1}{n} \sum_{i=1}^{n}X_i. \]

Under independent sampling,\[ SE(\bar{X}) = \frac{\sigma}{\sqrt{n}}. \]

If the sample mean has an approximately normal sampling distribution and \(\sigma\) is known, a confidence interval takes the form\[ \bar{X} \pm z^* \frac{\sigma}{\sqrt{n}}. \]

When \(\sigma\) is unknown and estimated using \(s\), the conventional interval is\[ \boxed{ \bar{X} \pm t^*_{n-1} \frac{s}{\sqrt{n}} }. \]

This analytic approach is convenient because the standard error and reference distribution have known formulas.

General Parameters and Estimators

Suppose the parameter of interest is\[ \theta \]

and its estimator is\[ \hat{\theta}. \]

The parameter could be:

A mean
A median
A percentile
A correlation
A regression coefficient
A ratio
A difference between groups
A prediction error
A nonlinear function of several parameters

The sampling distribution of \(\hat{\theta}\) describes how the estimator varies across repeated samples from the population.

Its standard error is\[ \boxed{ SE(\hat{\theta}) = \sqrt{ \operatorname{Var}(\hat{\theta}) } }. \]

When Analytic Formulas Become Difficult

For some estimators, there is no simple expression analogous to\[ SE(\bar{X})=\frac{\sigma}{\sqrt{n}}. \]

Difficulties can arise when:

The estimator is nonlinear.
The sampling distribution is skewed.
The statistic involves optimization.
The estimator is a ratio.
The parameter lies near a boundary.
The data contain censoring or truncation.
The model is hierarchical.
The sample is clustered or dependent.
The estimator is selected adaptively.
The sample size is too small for a normal approximation.

In such situations, a confidence interval of the form\[ \hat{\theta}\pm z^*SE(\hat{\theta}) \]

may be inaccurate even if the standard error can be estimated.

Why the Standard Normal Formula Can Fail

A normal-theory interval assumes that the standardized estimator is approximately normal:\[ \frac{ \hat{\theta}-\theta }{ SE(\hat{\theta}) } \approx N(0,1). \]

This can fail for several reasons.

Skewness

The sampling distribution may have a long tail on one side.

Bias

The estimator may not be centered at the parameter:\[ E(\hat{\theta})\neq\theta. \]

Boundaries

A proportion must remain between 0 and 1, but a symmetric normal interval can extend outside those limits.

Heavy Tails

Extreme sampling outcomes may occur more frequently than the normal model predicts.

Discreteness

A count statistic may take only a small number of possible values.

Instability

A statistic such as a ratio can change dramatically when its denominator is near zero.

Computational methods allow the sampling distribution’s shape to be approximated rather than automatically replaced with a normal curve.

The Simulation Principle

Suppose we can repeatedly generate datasets under an appropriate model.

For simulated dataset \(b\), compute\[ \hat{\theta}^{(b)}. \]

After \(B\) repetitions, we obtain\[ \hat{\theta}^{(1)}, \hat{\theta}^{(2)}, \ldots, \hat{\theta}^{(B)}. \]

These simulated estimates approximate the sampling distribution of \(\hat{\theta}\).

From them, we can estimate:

Expected value
Bias
Standard error
Quantiles
Confidence-interval coverage
Type I error
Power
Prediction error

Monte Carlo Simulation

Monte Carlo simulation generates data from a specified probability model.

Suppose\[ X\sim F_\eta, \]

where \(F_\eta\) is a probability distribution with known or selected parameter values \(\eta\).

A Monte Carlo study proceeds as follows:

Generate a dataset from \(F_\eta\).
Calculate \(\hat{\theta}\).
Repeat independently \(B\) times.
Examine the distribution of the resulting estimates.

The Monte Carlo estimate of the expected estimator value is\[ \overline{\hat{\theta}} = \frac{1}{B} \sum_{b=1}^{B} \hat{\theta}^{(b)}. \]

The estimated standard error is\[ \boxed{ \widehat{SE}_{MC} = \sqrt{ \frac{ \sum_{b=1}^{B} \left( \hat{\theta}^{(b)} – \overline{\hat{\theta}} \right)^2 }{ B-1 } } }. \]

What Monte Carlo Simulation Requires

Monte Carlo simulation requires a model from which data can be generated.

For example, we might assume\[ X_i\sim N(\mu,\sigma^2) \]

or\[ X_i\sim\operatorname{Bernoulli}(p). \]

The method can then show how an estimator behaves when that model is true.

Its conclusions are model-dependent. If the assumed data-generating process is unrealistic, the simulation may precisely describe the wrong situation.

The Bootstrap

The bootstrap is designed for situations in which the population distribution is unknown but an observed sample is available.

Suppose the data are\[ X_1,X_2,\ldots,X_n. \]

The nonparametric bootstrap treats the empirical distribution of these observations as an approximation to the unknown population distribution.

Each observed value receives probability\[ \frac{1}{n}. \]

A bootstrap sample is created by drawing \(n\) observations with replacement from the original sample.

Bootstrap Algorithm

For bootstrap repetition \(b\):

Draw\[ X_1^{*(b)},\ldots,X_n^{*(b)} \]with replacement from the observed sample.
Calculate the bootstrap estimate\[ \hat{\theta}^{*(b)}. \]
Repeat for\[ b=1,\ldots,B. \]

The collection\[ \hat{\theta}^{*(1)}, \ldots, \hat{\theta}^{*(B)} \]

forms the bootstrap distribution.

The bootstrap estimate of the standard error is\[ \boxed{ \widehat{SE}_{\text{boot}} = SD( \hat{\theta}^{*(1)}, \ldots, \hat{\theta}^{*(B)} ) }. \]

Why Sampling with Replacement Matters

A bootstrap sample has the same size as the original sample, but because sampling is performed with replacement:

Some original observations appear multiple times.
Some appear once.
Some do not appear at all.

This reproduces sample-to-sample variation relative to the empirical distribution.

Sampling without replacement would merely rearrange the original observations. Any statistic unaffected by order would remain unchanged, producing no bootstrap variability.

Monte Carlo and Bootstrap Compared

Feature	Monte Carlo simulation	Nonparametric bootstrap
Source of generated data	Specified probability model	Observed empirical distribution
True parameter known?	Often known in a simulation study	Usually unknown
Primary goal	Study behavior under assumed scenarios	Estimate sampling uncertainty from observed data
Main dependence	Quality of the simulation model	Representativeness of the sample
Typical output	Bias, power, coverage, expected error	Standard error, confidence interval, bias estimate

Both use repeated random computation, but they should not be treated as identical methods.

Parametric Bootstrap

A parametric bootstrap combines aspects of both approaches.

First, fit a probability model:\[ F_{\hat{\eta}} \]

to the observed data.

Then:

Generate samples from \(F_{\hat{\eta}}\).
Refit the model to each sample.
Calculate the statistic.
Examine the resulting distribution.

The parametric bootstrap can be more efficient than the nonparametric bootstrap when the fitted model is correct, but it is more vulnerable to model misspecification.

Estimating a Standard Error

Suppose no convenient formula exists for \(SE(\hat{\theta})\).

The bootstrap approximates it with\[ \widehat{SE}_{\text{boot}} = \sqrt{ \frac{ \sum_{b=1}^{B} \left( \hat{\theta}^{*(b)} – \overline{\hat{\theta}^*} \right)^2 }{ B-1 } }, \]

where\[ \overline{\hat{\theta}^*} = \frac{1}{B} \sum_{b=1}^{B} \hat{\theta}^{*(b)}. \]

This can be used in a normal-style interval:\[ \hat{\theta} \pm z^* \widehat{SE}_{\text{boot}}. \]

However, if the bootstrap distribution is strongly skewed, a symmetric interval may still be inappropriate.

Bootstrap Percentile Interval

The percentile method uses quantiles of the bootstrap distribution directly.

For a nominal 95% interval:

Sort the bootstrap estimates.
Find the 2.5th percentile.
Find the 97.5th percentile.

The interval is\[ \boxed{ \left( q^*_{0.025}, q^*_{0.975} \right) }. \]

This method allows asymmetric intervals, but it does not automatically correct bias or other distortions.

Basic Bootstrap Interval

The basic interval reflects the bootstrap quantiles around the original estimate:\[ \boxed{ \left( 2\hat{\theta}-q^*_{0.975}, \, 2\hat{\theta}-q^*_{0.025} \right) }. \]

This reverses the bootstrap error distribution around \(\hat{\theta}\).

BCa Interval

The bias-corrected and accelerated interval adjusts for:

Bias in the bootstrap distribution
Changes in the estimator’s standard error across parameter values
Asymmetry

BCa intervals often perform better than simple percentile intervals, although they are more computationally involved.

No bootstrap interval is universally best. The method should match the estimator and inferential setting.

Example: Estimating a Median

For a sample mean, the formula\[ s/\sqrt{n} \]

is readily available.

For a sample median, the standard error depends on the shape and density of the population distribution near the median. A simple universal formula is not available without additional assumptions.

A nonparametric bootstrap can:

Resample the observed data.
Calculate a median for each resample.
Use the standard deviation of the bootstrap medians as an estimated standard error.
Use bootstrap quantiles to construct an interval.

This illustrates the main benefit of computational inference: the same resampling logic can be applied to many statistics.

Example: A Ratio

Suppose the parameter is\[ \theta=\frac{\mu_X}{\mu_Y}. \]

The estimator is\[ \hat{\theta} = \frac{\bar{X}}{\bar{Y}}. \]

Its sampling distribution may be skewed, particularly when \(\bar{Y}\) can be close to zero.

A symmetric normal interval may perform poorly. Bootstrap resampling can approximate the distribution of the ratio directly, although ratio parameters near a zero denominator remain intrinsically unstable.

Bootstrap in Regression

For a regression model,\[ Y_i=\beta_0+\beta_1X_i+\varepsilon_i, \]

the bootstrap can estimate uncertainty in:

Regression coefficients
Predicted values
Variable-selection procedures
Nonlinear functions of coefficients
Model-performance measures

Several resampling methods are possible.

Case Resampling

Resample complete observations:\[ (X_i,Y_i). \]

This treats the observed cases as a sample from a joint population.

Residual Resampling

Fit the regression model and resample residuals, then construct new responses:\[ Y_i^* = \hat{Y}_i+\varepsilon_i^*. \]

This relies more strongly on the regression model and commonly assumes similarly distributed errors.

Wild Bootstrap

Randomly perturb residuals in a way designed to accommodate heteroscedasticity.

Cluster Bootstrap

Resample entire clusters when observations are grouped within schools, hospitals, households, or other units.

The resampling scheme must preserve the dependence structure relevant to the original design.

The Bootstrap Does Not Make Assumptions Disappear

The bootstrap replaces some analytic derivations with computation, but it still requires assumptions.

The ordinary nonparametric bootstrap assumes that:

The observed sample is representative.
Sampling units are independent and identically distributed.
The empirical distribution approximates the population adequately.
The statistic behaves sufficiently smoothly.
The sample includes the important features of the population.

If observations are dependent, ordinary individual resampling is generally invalid.

Situations Where the Ordinary Bootstrap Can Fail

The bootstrap may perform poorly for:

Extremely small samples
Rare-event probabilities
Sample maxima and minima
Heavy-tailed distributions
Parameters on boundaries
Nonsmooth estimators
Ratios with near-zero denominators
Dependent observations
Data selected through complex adaptive processes

Modified methods may include:

Block bootstrap
Cluster bootstrap
Stratified bootstrap
Subsampling
\(m\)-out-of-\(n\) bootstrap
Parametric bootstrap
Wild bootstrap

Simulation Is Not Automatically Better

Computational methods can outperform crude normal approximations when they represent the sampling process more accurately.

However, simulation is not automatically superior to an analytic method.

An analytic t-interval for a normally distributed mean may be:

Exact under its assumptions
Faster
Easier to interpret
More accurate than a poorly designed bootstrap

Computational methods are advantageous when they solve a real inferential difficulty, not merely because they use more computation.

Monte Carlo Error

Using a finite number \(B\) of simulations introduces Monte Carlo error.

If a probability \(q\) is estimated from \(B\) simulation indicators, its Monte Carlo standard error is approximately\[ \boxed{ MCSE = \sqrt{ \frac{ \hat q(1-\hat q) }{ B } } }. \]

For example, if\[ \hat q=0.95 \]

and\[ B=10{,}000, \]

then\[ MCSE = \sqrt{ \frac{ 0.95(0.05) }{ 10{,}000 } } \approx 0.00218. \]

Increasing \(B\) reduces numerical simulation error but does not correct model or sampling bias.

Statistical Error Versus Computational Error

Computational inference contains at least two distinct uncertainties:

Statistical Uncertainty

This arises because only one finite sample was observed.

Monte Carlo Uncertainty

This arises because only a finite number of simulated or bootstrap samples were generated.

Increasing the number of bootstrap repetitions reduces Monte Carlo uncertainty. It does not create new information about the population.

Reproducibility

A computational analysis should document:

Random seed
Number of repetitions
Resampling unit
Sample size per repetition
Statistic recalculated
Confidence-interval method
Treatment of missing data
Handling of dependence
Software and version
Monte Carlo error when relevant

The statistic must be recomputed completely in every resample, including model selection or tuning steps that are part of the original procedure.

Key Takeaway

Classical inference often uses\[ \text{estimate} \pm \text{critical value}\times SE. \]

This works well when the standard error has a reliable formula and the estimator has an appropriate reference distribution.

When those conditions fail, repeated simulation can approximate the sampling distribution.

Monte Carlo simulation generates data from a specified probability model.
The nonparametric bootstrap resamples from the observed empirical distribution.
The parametric bootstrap generates data from an estimated probability model.

These methods can estimate standard errors, bias, confidence intervals, power, and coverage for complicated estimators. They do not eliminate the need for valid sampling, appropriate dependence modeling, and sound assumptions. More computation reduces numerical error, but it cannot repair an incorrect data-generating model or a biased sample.