Simulation, Bootstrap and Randomization Methods

Introduction

Sampling distributions are central to statistical quality management because they describe how an estimator varies over repeated samples. In practice, however, exact sampling distributions are often unavailable: the statistic may be nonlinear, the population may be non-normal, the sample may be small, or the estimator may involve a complex production rule. Simulation and resampling methods provide flexible computational alternatives.

Monte Carlo simulation generates repeated samples from a specified probability model. Bootstrap methods approximate repeated sampling by resampling observed data, either from an estimated parametric model or directly from the empirical distribution. Randomization and permutation methods generate reference distributions under a null hypothesis by rearranging labels or signs. Cross-validation uses repeated data splitting to estimate predictive performance and its variability.

These methods are not substitutes for statistical reasoning. Their validity depends on the data-generating model, independence assumptions, exchangeability, the resampling design, and faithful implementation. Reproducibility therefore requires documented code, fixed random seeds, explicit preprocessing rules, and preservation of the original data.

Learning objectives

After completing this chapter, the reader should be able to:

1. Explain why computational sampling distributions are useful.
2. Construct Monte Carlo, parametric-bootstrap, nonparametric-bootstrap, permutation, and cross-validation procedures.
3. Estimate bias, standard error, confidence intervals, p-values, and predictive performance from simulated distributions.
4. Select a resampling method appropriate to the research design.
5. Implement reproducible analyses and diagnose simulation error.
6. Evaluate limitations, governance risks, and managerial implications in quality-improvement studies.

1. The computational view of a sampling distribution

Let X₁, …, Xₙ denote observations and let T(X₁, …, Xₙ) be a statistic. Its sampling distribution is the probability distribution of T over repeated samples from the target population. If an analytical distribution is unavailable, generate B simulated or resampled statistics T₁, …, Tᵦ and use their empirical distribution.

The empirical cumulative distribution function is

where I(·) equals 1 when its condition is true and 0 otherwise. The approximation improves as B increases, although it can still be systematically wrong if the assumed data-generating mechanism is inappropriate.

For a simulated distribution, the estimated mean and standard deviation are

The first measures the centre of the computational distribution; the second estimates Monte Carlo variability of the statistic under the generating procedure. These are distinct from uncertainty caused by observing only one sample.

A useful approximation to the Monte Carlo standard error of a simulated mean is

For a simulated proportion p̂, the corresponding approximation is

These formulas assume independent simulation replicates. For Markov-chain simulation or other dependent algorithms, an effective sample size is required instead of B.

2. Monte Carlo simulation

2.1 Principle

Monte Carlo simulation samples from a fully specified probability model. Suppose a quality characteristic Y is assumed to follow a normal distribution with mean μ and standard deviation σ. A simulation replicate consists of generating n observations from that distribution and calculating the statistic of interest, such as the sample mean, defect proportion, or capability index.

The generic algorithm is:

1. Specify the population model and parameter values.
2. Set a documented random seed.
3. Generate a sample of size n.
4. Calculate the statistic T.
5. Repeat steps 3–4 B times.
6. Summarize and plot the resulting statistics.
7. Compare simulation results with analytical expectations or a second implementation.

For example, to study the distribution of a sample mean under a normal process, generate Yᵢ ~ N(μ, σ²), calculate Ȳ for each replicate, and inspect the simulated means. Under independent observations, the theoretical standard error is σ/√n. Agreement between the simulated standard deviation and this value is a useful validation check.

2.2 Simulation design choices

The model must represent the operational context. Important decisions include:

• Distribution: normal, binomial, Poisson, lognormal, gamma, or an empirically estimated distribution.
• Dependence: serial correlation, clustering, repeated measurements, or batch effects.
• Heterogeneity: changes in variance or mean across shifts, machines, suppliers, or operators.
• Finite population: sampling without replacement rather than independent sampling.
• Rare events: sufficiently large B to observe tail behaviour.

A simulation is a model-based experiment, not evidence that the model is true. Sensitivity analysis should vary plausible parameter values and structural assumptions.

2.3 Worked example: simulated defect rate

Assume a stable process has defect probability p = 0.02 and inspection samples contain n = 100 units. The number of defects D follows a binomial distribution:

The sample defect rate is p̂ = D/n. Its theoretical mean and standard deviation are

Thus,

A Monte Carlo study with B = 100,000 replicates should produce a mean near 0.02 and a standard deviation near 0.0140. Because the distribution is discrete and bounded between 0 and 1, a normal approximation may be poor in the extreme tails; the simulated distribution retains the correct discreteness.

3. Bootstrap methods

The bootstrap treats the observed sample as an approximation to the population and mimics repeated sampling from it. Let the observed data be x₁, …, xₙ and let the statistic be T̂ = T(x₁, …, xₙ). A bootstrap sample contains n observations drawn with replacement from the observed data. Applying T to B such samples gives T₁, …, Tᵦ.

The bootstrap standard error is

and the estimated bootstrap bias is

The method assumes that the empirical distribution is a reasonable approximation to the population and that the resampling unit reflects the study design.

3.1 Nonparametric bootstrap

The ordinary nonparametric bootstrap resamples individual observations with replacement. It is attractive when the population distribution is unknown and the statistic is complicated. It is commonly used for means, medians, correlations, regression coefficients, reliability measures, and process-performance indices.

Its basic algorithm is:

1. Retain the original sample and calculate T̂.
2. Draw n observations with replacement from the sample.
3. Calculate T*.
4. Repeat B times.
5. Use the bootstrap distribution for standard errors, bias assessment, or confidence intervals.

A percentile interval uses empirical quantiles:

where qᵣ is the r quantile of the bootstrap statistics. For a 95% interval, α = 0.05 and the endpoints are the 2.5th and 97.5th percentiles.

The percentile interval is intuitive but may perform poorly with strong bias or skewness. Alternative intervals include the basic interval, studentized interval, and bias-corrected and accelerated (BCa) interval. BCa intervals account, approximately, for bias and acceleration but require additional calculations and careful software implementation.

3.2 Parametric bootstrap

The parametric bootstrap assumes a family of distributions. First estimate parameters from the observed data, then generate bootstrap samples from the fitted model. For a normal model, estimate μ̂ and σ̂, generate Y*ᵢ ~ N(μ̂, σ̂²), and recompute the statistic.

Parametric bootstrap is often more efficient than nonparametric bootstrap when the model is credible, especially for tail probabilities, small samples, or statistics whose sampling behaviour depends strongly on a fitted distribution. It is vulnerable to model misspecification. Goodness-of-fit diagnostics and sensitivity analyses are therefore essential.

3.3 Structured bootstrap

Independent resampling is invalid when observations are dependent or clustered. Appropriate alternatives include:

• Cluster bootstrap: resample entire hospitals, plants, suppliers, or production lots.
• Stratified bootstrap: resample within defined strata while preserving their structure.
• Paired bootstrap: resample paired observations together.
• Block bootstrap: resample contiguous time blocks for serially correlated data.
• Residual bootstrap: resample model residuals when the fitted regression structure is retained.

The principle is simple: resample the unit that was independently sampled. Resampling individual measurements from a batch study can artificially destroy batch variation and produce overly narrow intervals.

4. Randomization and permutation tests

A permutation test constructs a null distribution by rearranging observed labels or values in ways that are valid under the null hypothesis. Consider two independent groups with outcome means X̄₁ and X̄₂. The observed difference is

Under the null hypothesis of no group effect, group labels are exchangeable. Pool all observations, randomly assign the original group sizes, calculate the difference, and repeat B times. The two-sided randomization p-value can be estimated as

The +1 correction avoids a reported p-value of zero and gives a conservative finite-randomization estimate. The test is valid when the exchangeability assumption is appropriate. Under unequal variances, dependence, or covariate imbalance, unrestricted permutation may not be valid; restricted permutations, studentized statistics, paired sign permutations, or model-based methods may be needed.

For paired data, differences are calculated within pairs and signs may be randomly flipped under a symmetric null. For regression, permutation can involve residuals or treatment labels, but the design must preserve nuisance structure.

A p-value is not the probability that the null hypothesis is true. It is the probability, under the null and the specified randomization scheme, of obtaining a result at least as extreme as the observed result.

5. Cross-validation as a resampling method

Cross-validation estimates how well a fitted model predicts new observations. In K-fold cross-validation, the data are partitioned into K folds. Each fold is used once as a validation set while the remaining K−1 folds form the training set.

For loss function L and prediction f̂₋ₖ trained without fold k, the cross-validation error is

For squared-error loss, L(y, ŷ) = (y − ŷ)²; for classification, common choices include misclassification loss or log loss. The folds must be created before model fitting, and all preprocessing—including scaling, feature selection, and imputation—must be learned inside each training fold. Otherwise information leaks from validation data into the model.

Common designs include:

• Leave-one-out cross-validation: K = n; low bias but potentially high variance and high computation.
• 5- or 10-fold cross-validation: often practical compromises.
• Repeated K-fold: repeats fold assignment to assess sensitivity.
• Stratified K-fold: preserves class proportions.
• Group K-fold: keeps observations from the same subject, plant, or batch together.
• Time-series validation: trains on earlier observations and validates on later observations; random folds are usually inappropriate.

Cross-validation estimates predictive performance, not causal effect and not necessarily future performance under process drift. A final untouched test set is valuable when model selection is extensive.

6. Reproducible implementation

A reproducible computational analysis should record the data version, software and package versions, random-number generator, seed, preprocessing rules, resampling design, number of replicates, and analysis code. A seed makes a computation repeatable, but it does not make a flawed analysis valid.

A language-independent pseudocode pattern is:

set_seed(2025)
read_and_validate_data()
define_statistic(data)
for b in 1,...,B:
    resample according to the study design
    store statistic(resampled_data)
summarize stored statistics
plot distribution and diagnostics
save code, settings, results, and session information

Best practice includes pre-specifying the estimand and stopping rules, separating exploratory from confirmatory analyses, checking missingness, retaining failed iterations for diagnosis, and using independent code or analytical benchmarks for validation. Parallel computation is acceptable, but random-number streams must be managed so that results do not depend unpredictably on processor count.

Choosing B

There is no universal number of replicates. For a percentile confidence interval, the tails require more replicates than a central standard error. If B = 1,000, the nominal 2.5th percentile is represented by only about 25 observations. B = 10,000 or more is often reasonable for ordinary intervals, while rare-event probabilities may require substantially more. Report B and, where relevant, the Monte Carlo error of the reported quantile or p-value.

7. Illustrative case: supplier coating thickness

This is a fictional illustrative case.

Background and problem

A manufacturer receives coated metal panels from two suppliers. Coating thickness is measured on five panels per shipment, and several shipments are observed. Management wants to know whether Supplier B has a higher median thickness and whether a predictive model can identify shipments at risk of falling below the lower specification limit.

Analysis

Measurements within a shipment are correlated because they share coating conditions. Treating all individual measurements as independent would understate uncertainty. The shipment, not the individual panel, is the primary resampling unit.

The analyst first calculates the observed difference in supplier medians. A cluster bootstrap resamples complete shipments within each supplier, retaining all five panel measurements within selected shipments. For each replicate, the difference in medians is calculated. A percentile or BCa interval describes uncertainty in the supplier contrast.

For a hypothesis test, supplier labels are permuted at the shipment level, preserving the within-shipment structure. If operational covariates such as coating line are imbalanced, unrestricted permutation may be misleading; restricted permutation within line or a regression-based analysis is preferable.

To predict out-of-specification shipments, the analyst uses group-wise cross-validation, ensuring that panels from one shipment never appear in both training and validation sets. If future deployment is chronological, time-ordered validation is added. Preprocessing and variable selection occur within each training fold.

Outcome and interpretation

Suppose the bootstrap interval for the median difference excludes zero, while the permutation p-value is small under the prespecified design. The evidence supports a supplier difference under the target population and assumptions, but it does not establish that supplier identity causes the difference if coating line and material grade are confounded. Predictive performance is reported using a validation loss and interval across folds, not merely training accuracy.

Learning points

1. The resampling unit must match the independent sampling unit.
2. Bootstrap estimation and permutation testing answer different questions.
3. Predictive validation must reproduce the intended deployment environment.
4. Statistical significance, practical significance, and causal interpretation are distinct.

8. Limitations, ethics, and governance

Resampling cannot repair biased sampling, systematic measurement error, missing-not-at-random data, or an inappropriate estimand. A bootstrap distribution may be misleading when the sample is tiny, contains influential observations, or does not represent the target population. A permutation test may fail when exchangeability is violated. Cross-validation may be optimistic under leakage or pessimistic under a deployment environment unlike the study sample.

In quality management, results can influence supplier qualification, employee evaluation, warranty decisions, and customer safety. Analysts should document uncertainty, avoid selective reporting, protect confidential supplier and worker data, and maintain audit trails. Automated or AI-assisted code should be independently reviewed: generated code can silently use the wrong resampling unit, leak information, or report an inappropriate confidence interval. Governance should require human accountability, version control, access control, and reproducible execution.

Contemporary developments include probabilistic programming, approximate Bayesian computation, distributed bootstrap, conformal prediction, and resampling for complex causal and machine-learning estimands. These tools increase flexibility but also increase the need to state assumptions and validate computational procedures.

Chapter Summary

Simulation and resampling methods approximate sampling distributions when analytical solutions are difficult or unreliable. Monte Carlo simulation samples from a specified model; parametric bootstrap samples from a fitted model; nonparametric bootstrap resamples observed data; permutation tests construct null distributions through valid rearrangements; and cross-validation estimates predictive performance through repeated training and validation. Each method requires assumptions about independence, exchangeability, clustering, time, and model structure. Correct resampling units, leakage-free preprocessing, adequate replication, Monte Carlo diagnostics, and transparent reporting are essential. Computational sophistication does not compensate for biased data or weak design.

Key Terms

Sampling distribution: Distribution of a statistic over repeated samples.

Monte Carlo simulation: Repeated random generation from a specified probability model.

Bootstrap: Resampling method used to approximate repeated sampling from observed data.

Parametric bootstrap: Bootstrap based on a fitted distributional model.

Nonparametric bootstrap: Bootstrap based on the empirical distribution of observations.

Bootstrap standard error: Standard deviation of bootstrap statistic values.

Bootstrap bias: Difference between the bootstrap mean and the original estimate.

Percentile interval: Confidence interval formed from bootstrap quantiles.

BCa interval: Bias-corrected and accelerated bootstrap interval.

Permutation test: Null-distribution test based on rearranging labels or values.

Exchangeability: Invariance of the relevant joint distribution under permitted rearrangements.

Randomization p-value: Tail probability calculated from a randomization distribution.

Cluster bootstrap: Resampling of intact groups or clusters.

Block bootstrap: Resampling of blocks to preserve serial dependence.

Cross-validation: Repeated partitioning for estimating predictive performance.

Data leakage: Use of validation information during model training.

Monte Carlo standard error: Simulation error due to a finite number of replicates.

Estimand: Precisely defined quantity targeted by an analysis.

Reproducibility: Ability to recreate results using documented data, code, and settings.

Review Questions

Short Answer Questions

1. Distinguish Monte Carlo simulation from nonparametric bootstrap.
2. What does the bootstrap standard error estimate?
3. Why is exchangeability important in a permutation test?
4. Define data leakage in cross-validation.
5. Why should the resampling unit match the independent sampling unit?

Long Answer Questions

1. Explain the assumptions, algorithm, interpretation, and limitations of a nonparametric bootstrap confidence interval.
2. Compare parametric bootstrap, nonparametric bootstrap, and permutation testing.
3. Describe how serial correlation and clustering affect resampling design.
4. Explain how cross-validation should be designed for a time-dependent quality-monitoring problem.
5. Develop a reproducibility protocol for a simulation-based quality study.

Analytical / Application Questions

1. A study samples ten products from each of six production lots. Explain why individual-level bootstrap resampling may be inappropriate and propose a better design.
2. A classification model has excellent training accuracy but poor performance in deployment. Design a leakage-free cross-validation investigation.
3. A supplier comparison has unequal variances and unequal shipment counts. Propose a defensible testing and uncertainty-estimation strategy, including assumptions and diagnostics.

Multiple Choice Questions with Answers

Question 1

Which method generates samples from a specified probability model to study a statistic?

A. Monte Carlo simulation

B. Nonparametric bootstrap

C. Cross-validation

D. Permutation testing

Correct Answer: A. Monte Carlo simulation

Question 2

In an ordinary nonparametric bootstrap sample of size n, observations are drawn:

A. Without replacement from a fitted normal distribution

B. With replacement from the observed data

C. Only from the largest cluster

D. By permuting the response labels once

Correct Answer: B. With replacement from the observed data

Question 3

The bootstrap standard error is usually estimated by the:

A. Mean of the original observations

B. Variance of the raw population

C. Standard deviation of bootstrap statistics

D. Number of bootstrap replicates

Correct Answer: C. Standard deviation of bootstrap statistics

Question 4

A permutation test is most directly justified when the relevant observations or labels are:

A. Exchangeable under the null hypothesis

B. Normally distributed in every group

C. All measured without error

D. Sorted in ascending order

Correct Answer: A. Exchangeable under the null hypothesis

Question 5

Which resampling design is generally appropriate for correlated observations within production lots?

A. Resample individual observations independently

B. Resample complete lots or use a cluster bootstrap

C. Delete all observations from small lots

D. Use only the largest observation per lot

Correct Answer: B. Resample complete lots or use a cluster bootstrap

Question 6

In K-fold cross-validation, each observation is normally used for validation:

A. In every fold

B. Exactly once

C. Only if it is an outlier

D. Never

Correct Answer: B. Exactly once

Question 7

Which action creates data leakage?

A. Fitting a scaler separately inside each training fold

B. Keeping a final test set untouched

C. Selecting features using the full dataset before cross-validation

D. Using group-wise folds for clustered data

Correct Answer: C. Selecting features using the full dataset before cross-validation

Question 8

If B independent simulation replicates produce statistics with standard deviation s, the approximate Monte Carlo standard error of their mean is:

A. s × B

B. s ÷ B

C. s × √B

D. s ÷ √B

Correct Answer: D. s ÷ √B

Question 9

What is the main risk of a parametric bootstrap?

A. It cannot calculate a statistic

B. It assumes the fitted distribution is adequate

C. It always destroys clustering

D. It requires no model assumptions

Correct Answer: B. It assumes the fitted distribution is adequate

Question 10

Why is a +1 correction commonly used in a permutation p-value numerator and denominator?

A. To guarantee a p-value of exactly 0.5

B. To avoid reporting a zero p-value from a finite randomization sample

C. To remove all Monte Carlo error

D. To make every test one-sided

Correct Answer: B. To avoid reporting a zero p-value from a finite randomization sample