Advanced Monte Carlo Methods
Learning Objectives
To be able to explain and apply the basic principles of bootstrap sampling
To be able to apply variance reduction techniques to simulation experiments
To be able to apply the main techniques recommended for generating multi-variate random variables
To be able to apply Markov Chain Monte Carlo (MCMC) techniques to the generation of random variables
What is Bootstrapping?
Bootstrapping is a statistical procedure that resamples from a sample to create many simulated samples.
It is generally a non-parametric Monte Carlo method that estimates characteristics of a population by resampling.
Recall a random sample \(x_1, x_2,\cdots,x_n\): the \(X_i\) are independent and every \(X_i\) has the same distribution, drawn from some population.
The fundamental principle: we assume the original random sample is essentially the population. Sampling from the sample is, in essence, sampling from the population.
If the original sample is representative of the true population, then statistics computed from resamples provide observations of the sampling distribution of the estimator.
The Bootstrap Principle and the Empirical Distribution
Let \(F\) be the CDF of \(X\), and let \(\theta = t(F)\) be a population parameter of interest (mean, median, correlation, \(\dots\)).
We usually do not know \(F\), but the sample \(\vec{x}=(x_1,\cdots,x_n)\) gives the empirical distribution \(F_n\), an estimate of \(F\).
Define \(X^{*}\) as a value drawn at random from \((x_1,\cdots,x_n)\), so \(P(X^{*}=x_i)=1/n\). A bootstrap sample \(\vec{x}^{*}=(x^{*}_1,\cdots,x^{*}_n)\) is drawn with replacement — values may repeat.
Let \(\hat{\theta}^{*}_b\) be the estimator evaluated on the \(b^{th}\) bootstrap sample, \(b=1,\cdots,B\). These are the bootstrap replicates.
The empirical distribution of \((\hat{\theta}^{*}_1,\cdots,\hat{\theta}^{*}_B)\) approximates the sampling distribution of \(\hat{\theta}\) — its quality depends entirely on the quality of the original sample.
Standard Error, Bias, and Confidence Intervals
Bootstrap standard error — the sample standard deviation of the replicates: \[ \widehat{se(\hat{\theta})} = \sqrt{\tfrac{1}{B-1}\textstyle\sum_{b=1}^{B}\big(\hat{\theta}^{*}_{b} - \bar{\hat{\theta^{*}}} \big)^2}, \qquad \widehat{bias}(\hat{\theta}) = \bar{\hat{\theta^{*}}} - \hat{\theta} \]
CIs derive from the distribution of \(\hat{\theta} - \theta\): \(\;\hat{\theta}-\gamma_{1-\alpha/2} \leq \theta \leq \hat{\theta} - \gamma_{\alpha/2}\). Five approximate intervals are used:
- standard normal (SNBCI), basic (BBCI), percentile (BPCI), bias-corrected & accelerated (BCa), and bootstrap-t (BTCI)
The percentile interval is popular, but BCa is recommended — it corrects the percentile interval for bias and skewness. Bootstrap-t is next, but needs an extra inner resampling step (except when \(\hat{\theta}=\bar{x}\)).
Guidance on the number of resamples \(B\): at least in the hundreds; \(B \approx 40n\) has been suggested.
Bootstrapping “by Hand” with DPopulation
Bootstrapping is just sampling again and again from a discrete population. The KSL
DPopulationclass is a natural starting point.Treat the original size-10 sample as the pseudo-population and resample size 10 with replacement.
// mainSample = doubleArrayOf(6.0, 7.0, 5.0, 1.0, 0.0, 4.0, 6.0, 0.0, 6.0, 1.0)
val samplePopulation = DPopulation(mainSample)
val bootStrapAverages = mutableListOf<Double>()
for (i in 1..10) {
val bootStrapSample = samplePopulation.sample(10) // resample size n, with replacement
bootStrapAverages.add(bootStrapSample.average())
}
// percentile-based 90% bootstrap CI on the original estimate (x̄ = 3.6)
val lcl = Statistic.percentile(bootStrapAverages.toDoubleArray(), 0.05)
val ucl = Statistic.percentile(bootStrapAverages.toDoubleArray(), 0.95)
val ci = Interval(lcl, ucl)- The resample size must equal the original sample size — the CI is a statement about the original estimate.
The KSL Bootstrap Class
- Requires the original
DoubleArrayand aBSEstimatorIfc(a functional interface; built-ins: average, variance, median, min, max). It computes the replicates, bias, standard error, and all five confidence intervals.
original estimate = 3.60 bias estimate = -0.0518 std. err. = 0.0401
percentile c.i. = [1.9025, 5.1975] BCa c.i. = [1.9000, 5.0629]Multiple Statistics and Case-Based Bootstrapping
BootstrapSampler(withMVBSEstimatorIfc) returns many estimates from one sample — e.g.,BasicStatistics()gives average, variance, min, max, skewness, kurtosis, lag-1 correlation/covariance.CaseBootstrapSampler(withCaseBootEstimatorIfc) samples case identifiers (row indices) — ideal for bootstrapping the rows of a regression data set.
statistic name = b0
original estimate = 10.9972 bias estimate = 0.0916
percentile c.i. = [9.5669, 12.9414] // covers the true intercept β0 = 10.0What is a Variance Reduction Technique (VRT)?
A VRT is a sampling strategy that reduces the variance of the sampling distribution of an estimator \(\hat{\theta}\) of some parameter \(\theta\).
For the crude estimator of the mean, \(\hat{\theta}=\bar{Y}\), with IID \(Y_i\): \(\;\text{Var}[\bar{Y}] = \sigma^2_Y / n\). A VRT tries to do better than this.
Variance reduction is about the estimator, not the process: \(\sigma^2_Y\) (the population variability) is unchanged — we change the variance of the sampling distribution.
Two families of techniques:
- Manipulating randomness: Common Random Numbers (CRN), Antithetic Variates (AV)
- Exploiting process knowledge: Indirect Estimation, Control Variates, Stratified Sampling, Conditional Expectation, Importance Sampling
Common Random Numbers (CRN)
- Used when comparing systems: \(\theta = \theta_1 - \theta_2\), with pairs \(D_i = X_i - Y_i\) and \(\hat{\theta}=\bar{D}=\bar{X}-\bar{Y}\). If each alternative sees the same randomness, differences reflect a true difference, not noise.
\[ \text{Var}[\bar{D}]=\frac{\sigma^{2}_{X}}{n} +\frac{\sigma^{2}_{Y}}{n} - 2\rho_{XY}\sigma_{X}\sigma_{Y} \;\;\Rightarrow\;\; V_{CRN} = V_{IND} - 2\rho_{XY}\sigma_{X}\sigma_{Y} \]
Reduction requires positive correlation within pairs (\(\rho_{XY}>0\)); it backfires under negative correlation. Induce positive correlation by feeding the same \(U_i\) through inverse transforms: \(X_i=F_X^{-1}(U_i),\,Y_i=F_Y^{-1}(U_i)\).
Synchronization best practices (the KSL does these automatically): inverse-transform generation, a dedicated stream per source of randomness, and resetting each stream’s sub-stream at the start of every replication.
Antithetic Variates (AV)
- Applies to a single model: induce negative correlation within pairs of runs \((Y_{2j-1}, Y_{2j})\).
\[ \text{Var}[\hat{Y}]=\frac{\sigma^2_Y}{n}(1+\rho) \;<\; \frac{\sigma^2_Y}{n} \quad\text{when } \rho<0 \]
Recipe: alternate the stream between runs — Run 1 uses \((u_1,\cdots,u_k)\), Run 2 uses \((1-u_1,\cdots,1-u_k)\), so \(X_{2j-1}=F^{-1}(u),\,X_{2j}=F^{-1}(1-u)\).
In Monte Carlo integration the antithetic sampler is built via
sampler.antitheticInstance(); each replication returns \(\hat{Y}_j=(Y_{2j-1}+Y_{2j})/2\):
- Within a DEDS model, simply set the
antitheticOptionproperty of theModelclass totrue.
Control Variates (CV): Theory
Observe a variable \(X\) with known mean \(\mu=E[X]\) that is correlated with \(Y\). Define \(Z_j = Y_j + C(X_j - \mu)\); then \(E[Z_j]=E[Y_j]\) — unbiased for any constant \(C\).
The optimal constant minimizes the variance: \(C^{*}=-\text{cov}(Y,X)/\text{Var}(X)\), giving \(\text{Var}[Z_j]=\text{Var}[Y_j](1-\rho^2_{XY})\).
\(\hat{C}^{*}\) has the form of a regression slope — so with \(E[Y|X{=}x]=\beta_0+\beta_1(x-\mu)\), the CV estimator of \(\theta=E[Y]\) is the regression intercept \(\hat{\beta}_0\). This extends to \(q\) controls via multiple regression.
Control Variates in the KSL
- The
ControlVariateDataCollectorcaptures the response and the control variates (each with its known mean) across replications, then builds the regression whose intercept is the CV estimate.
val cvCollector = ControlVariateDataCollector(model)
cvCollector.addResponse(palletWorkCenter.totalProcessingTime, "TotalTime")
cvCollector.addControlVariate(palletWorkCenter.processingTimeRV, (8.0 + 12.0 + 15.0) / 3.0, "PalletTime")
cvCollector.addControlVariate(palletWorkCenter.numPalletsRV, (100.0 * 0.8), "NumPallets")
model.simulate()
val regressionData = cvCollector.collectedData("TotalTime", 20) // 20 batches of size 5
val regressionResults = cvCollector.regressionResults(regressionData) Predictor parameter parameterSE TValue P-Values
0 Intercept 494.631368 2.792054 177.156788 0.000000- The CV estimate (intercept = 494.63) differs from the crude estimate (497.79) and has a smaller half-width. Check residual normality before trusting it.
Importance Sampling (IS)
Rewrite the integral with a proposal density \(w(x)\) and sample \(X_i \sim w(x)\): \[ \theta = \int_a^b g(x)\,dx = \int_a^b \frac{g(x)}{w(x)}w(x)\,dx, \qquad \hat{\theta}=\frac{1}{n}\sum_{i=1}^{n}\frac{g(X_i)}{w(X_i)} \]
Choose \(w(x)\) to put more samples where \(g(x)\) is large. For \(g(x)=x^2\) on \([0,1]\), using \(w(x)=2x\) instead of \(U(0,1)\) gives \(VR = \tfrac{1/72}{4/45}\ll 1\) — a dramatic reduction because \(w(x)\) is close in shape to \(g(x)\).
KSL support:
MC1DIntegrationandMCMultiVariateIntegrationlet you supply \(h(x)=g(x)/w(x)\), where \(w(x)\) is the sampler’s distribution.
Generating from a Bivariate Normal
\((X_1,X_2)\sim BVN(\mu_1,\sigma^2_1,\mu_2,\sigma^2_2,\rho)\) is generated from two independent standard normals \(Z_1,Z_2\): \[ X_1 = \sigma_1 Z_1 + \mu_1, \qquad X_2 = \sigma_2\big(\rho Z_1 + \sqrt{1-\rho^2}\,Z_2\big) + \mu_2 \]
The KSL
BivariateNormalRVimplements exactly this ingenerate():
- Key takeaway: we now have a procedure for generating correlated random variables — the basis for the general methods that follow.
The KSL Multi-Variate Sampling Framework
Multi-variate variables sub-class
MVRVariable(implementingMVRVariableIfcandMVSampleIfc) and overridegenerate(array: DoubleArray)— you are simply generating arrays of data.MVSampleIfcmirrors the 1-D interfaces:sample(): DoubleArray,sample(sampleSize): List<DoubleArray>,sampleByColumn(sampleSize). Thedimensionproperty gives the array size (2 for a BVN).The general \(d\)-dimensional multi-variate normal (
MVNormalRV) extends the BVN algorithm using the Cholesky decomposition \(\Sigma = CC^{T}\): \(\;X_i = \mu_i + \sum_{j=1}^{i} c_{ij}Z_j\).
Copulas and Sklar’s Theorem
Copula: a CDF \(C:[0,1]^d \to [0,1]\) with uniform marginals — the joint CDF of dependent uniforms \((U_1,\cdots,U_d)\).
Sklar’s Theorem: any joint CDF can be written as a copula of its marginals, \[ F(x_1,\cdots,x_d) = C\big(F_1(x_1),\cdots,F_d(x_d)\big) \] (unique when the marginals are continuous).
This separates the dependence structure (the copula) from the marginals. General recipe for correlated vectors:
- Specify a copula and its parameters
- Generate \(\vec{U}=(u_1,\cdots,u_d)\) from the copula
- Set \(Y_i = F_i^{-1}(u_i)\) — inverse transform each marginal
The Gaussian copula is extracted from the MVN: \(C^G_{\mathbf{P}}(\vec{u}) = \Phi_{\mathbf{P}}\big(\Phi^{-1}(u_1),\cdots,\Phi^{-1}(u_d)\big)\).
The Gaussian Copula in the KSL
MVGaussianCopulagenerates the dependent uniforms;MVGaussianCopulaRVcompletes the recipe by applying the user-supplied inverse-CDF marginals. The marginals need not be from the same family and may be discrete or continuous.
class MVGaussianCopulaRV(
val marginals: List<InverseCDFIfc>,
correlation: Array<DoubleArray>,
// ... stream parameters
) : MVRVariable(/* ... */) {
private val myCopula = MVGaussianCopula(correlation, stream)
override fun generate(array: DoubleArray) {
val u = myCopula.sample() // dependent uniforms from the copula
for (i in u.indices) {
array[i] = marginals[i].invCDF(u[i]) // inverse transform each marginal
}
}
}Autocorrelated Generation: NORTA and AR(1)
Ignoring correlation that is actually present in an input process can grossly distort performance estimates. The Normal-to-Anything (NORTA) transform generates correlated variates of any marginal: \[ X_i = F^{-1}\big(\Phi(Z_i)\big) \] Correlated \(Z_i\) produce correlated \(X_i\) (though \(\rho_x \neq \rho_z\) in general).
Induce correlation in the \(Z_i\) with a stationary AR(1) process: \(Z_i = \phi Z_{i-1} + \varepsilon_i\), with \(\varepsilon_i \sim N(0,1-\phi^2)\) and lag-1 correlation \(\phi\).
AR1CorrelatedRNStreampackages this as a stream you can hand to any random variable:
class AR1CorrelatedRNStream(
val lag1Corr: Double,
private val stream: RNStreamIfc = KSLRandom.nextRNStream(),
) : RNStreamIfc by stream {
// myX starts as N(0,1); errorVariance = 1 - lag1Corr^2
override fun randU01(): Double {
val e = KSLRandom.rNormal(0.0, errorVariance, stream)
myX = lag1Corr * myX + e
return Normal.stdNormalCDF(myX) // correlated U(0,1) via the normal CDF
}
}Why Markov Chain Monte Carlo (MCMC) Works
MCMC approximately generates variates from arbitrary distributions — even ones known only up to a multiplicative constant (crucial in Bayesian work) — and works for multi-variate targets.
The Monte Carlo estimator \(\hat{\theta}=\frac{1}{n}\sum_{i=1}^{n} h(X_i)\) converges to \(\theta\) even when the \(X_i\) are dependent, under mild conditions.
Idea: build a Markov chain whose limiting distribution is the target \(w(\vec{x})\) (often written \(\pi(\vec{x})\)). For an irreducible, aperiodic chain, \[ \lim_{n \to \infty}\frac{1}{n}\sum_{i=1}^{n}h(X_i)=\sum_{j=1}^{N}h(j)\pi_j \quad \text{w.p. } 1 \]
So instead of needing \(\pi_j\) in closed form, we run a chain that has \(\pi_j\) as its steady state.
Discrete Markov Chains
A chain has the Markov property — the next state depends only on the current state — with single-step transition probabilities \(p_{ij}=P\{X_{n+1}=j \mid X_n=i\}\).
DMarkovChain(inksl.utilities.random.markovchain) simulates a discrete-state chain from its transition matrix:
- The observed proportions \((0.278, 0.411, 0.311)\) converge to the true steady-state distribution \(\pi = (0.279, 0.410, 0.311)\).
The Metropolis-Hastings Algorithm
MCMC turns the usual problem around: given a desired \(\pi_j\), construct a transition scheme \(\mathbf{P}\) that produces it. Use an arbitrary proposal \(\mathbf{Q}\) and accept/reject moves.
Discrete acceptance probability — note it needs \(\pi\) only up to a constant (ratio form): \[ \alpha_{ij} = \min\!\bigg[\frac{\pi_j q_{ji}}{\pi_i q_{ij}},\,1\bigg] \]
Continuous / multi-variate: target \(f(\vec{x})\), proposal \(q(\vec{y}\mid\vec{x})\), acceptance ratio \[ \rho(\vec{x},\vec{y}) = \frac{q(\vec{x}\mid\vec{y})\,f(\vec{y})}{q(\vec{y}\mid\vec{x})\,f(\vec{x})}, \qquad \alpha(\vec{x},\vec{y})=\min[\rho,\,1] \]
Each step: propose \(\vec{Y}\sim q(\vec{y}\mid\vec{x})\); draw \(U\sim U(0,1)\); set \(X_{t+1}=\vec{Y}\) if \(U\le\alpha\), else stay at \(X_t\).
Samplers and the KSL Interfaces
MetropolisHastingsMVis supported by two interfaces:FunctionMVIfc— the target density \(f(\vec{x})\) viaf(x: DoubleArray)ProposalFunctionMVIfc—generateProposedGivenCurrent()andproposalRatio()for \(p_r(\vec{x},\vec{y})=q(\vec{x}\mid\vec{y})/q(\vec{y}\mid\vec{x})\)
Independence sampler: \(q(\vec{y}\mid\vec{x})=q(\vec{y})\), so \(\rho = \dfrac{q(\vec{x})\,f(\vec{y})}{q(\vec{y})\,f(\vec{x})}\). Good proposals resemble \(f\).
Random-walk sampler: symmetric proposal, so \(p_r = 1\) and \(\rho = f(\vec{y})/f(\vec{x})\), with \(\vec{Y}=\vec{x}+\vec{W}\), \(\vec{W}\) symmetric.
The MetropolisHastingsMV Class
MCMC Example: Independence Sampler
- Target: the ore-content joint density \(f(x_1,x_2)\) (zinc \(\in[0.5,1.5]\), iron \(\in[20,35]\)), implemented as a
FunctionMVIfc:
- Proposal
ExampleIndependencePFuses independentGeneralizedBetamarginals over the two ranges (proposalRatio= \(g(\vec{x})/g(\vec{y})\)). Run it:
Acceptance Statistics: avg = 0.3431
BatchStatistic{name='X_1', avg=0.9967, ci=[0.9452, 1.0482], lag-1 corr=0.497}
BatchStatistic{name='X_2', avg=26.081, ci=[25.635, 26.528], lag-1 corr=0.237}MCMC Example: Random-Walk Sampler
- A symmetric proposal makes
proposalRatioreturn1.0; the privategenYGivenX()resamples \(x_i+e_i\) until it lands back inside the variable’s interval.
object ExampleRandomWalkPF : ProposalFunctionMVIfc {
private val y1Interval = Interval(0.5, 1.5)
private val y2Interval = Interval(20.0, 35.0)
private val e1rv = UniformRV(-1.0, 1.0) // c1 = 1.0
private val e2rv = UniformRV(-5.0, 5.0) // c2 = 5.0
override val dimension: Int get() = 2
override fun proposalRatio(currentX: DoubleArray, proposedY: DoubleArray) = 1.0
override fun generateProposedGivenCurrent(currentX: DoubleArray): DoubleArray {
val y1 = genYGivenX(y1Interval, e1rv, currentX[0]) // y = x + e, kept in-bounds
val y2 = genYGivenX(y2Interval, e2rv, currentX[1])
return doubleArrayOf(y1, y2)
}
}Acceptance Statistics: avg = 0.7859
BatchStatistic{name='X_1', avg=0.9992, ci=[0.9940, 1.0044], lag-1 corr= 0.108}
BatchStatistic{name='X_2', avg=27.370, ci=[27.150, 27.590], lag-1 corr=-0.117}- Higher acceptance (0.79 vs 0.34) and lower batch correlation → better mixing than the independence sampler. Acceptance rates of 20–50% are common; the KSL reports batch statistics because MCMC output is correlated.
Summary: Key KSL Classes
Bootstrapping (
ksl.utilities.statistics):Bootstrap,BootstrapSampler,CaseBootstrapSampler— withBSEstimatorIfc,MVBSEstimatorIfc,CaseBootEstimatorIfc.Variance reduction: random number streams;
ReplicationDataCollectorandControlVariateDataCollector; Monte Carlo integration inksl.utilities.mcintegration(MC1DIntegration,antitheticOption).Multi-variate / correlated generation (
ksl.utilities.rvariable):MVRVariable/MVSampleIfc,BivariateNormalRV,MVNormalRV,MVGaussianCopula,MVGaussianCopulaRV,AR1CorrelatedRNStream.MCMC:
DMarkovChain(ksl.utilities.random.markovchain);MetropolisHastingsMVwithFunctionMVIfcandProposalFunctionMVIfc(ksl.utilities.random.mcmc).The KSL provides immediately useful implementations and a structure that makes it easy to add your own.