Monte Carlo, Latin Hypercube, and Stratified Sampling

Why the Sampling Design Matters

In PBPK uncertainty analysis, uncertain inputs are represented by distributions, and a sampling design determines which input combinations are evaluated. Even when two designs use the same number of draws from the same distributions, they can cover the input space differently. This can affect how efficiently the resulting sample represents the uncertainty relevant to a model prediction.

This tutorial compares three designs: random Monte Carlo sampling, Latin hypercube sampling, and stratified sampling. The focus is on how they select input values and what their designs guarantee. The consequences of passing those samples through a PBPK model depend on the model response; selecting a well-balanced input sample does not by itself guarantee accurate estimates of every output feature.

Definition: Let \(\boldsymbol{\Theta}=(\Theta_1,\ldots,\Theta_d)\) denote the uncertain PBPK inputs, with a specified joint distribution. A sampling design generates \(N\) input vectors from, or in a controlled way representative of, that target distribution. Random Monte Carlo draws are independent draws from the joint distribution. Latin hypercube sampling (LHS) divides each input's marginal probability scale into \(N\) equal-probability intervals and selects one value from each interval, then pairs the marginal samples into input vectors. Stratified sampling divides the target distribution into strata, samples within each stratum, and combines the stratum results using their probabilities as weights.

Random Monte Carlo Sampling

For random Monte Carlo sampling, each vector \(\boldsymbol{\Theta}^{(i)}\), for \(i=1,\ldots,N\), is drawn independently from the specified joint distribution of the inputs. If the inputs are independent, each component can be drawn from its marginal distribution. If they are dependent, the sampling procedure must draw from a joint distribution that represents that dependence; independently drawing from each marginal would instead impose independence.

For a continuous marginal distribution with cumulative distribution function \(F_j\), one common construction is to draw \(U_j\sim\mathrm{Uniform}(0,1)\) and set \(\Theta_j=F_j^{-1}(U_j)\), where \(F_j^{-1}\) is the quantile function. This construction describes a marginal draw. It does not, by itself, define how to generate a dependent vector of inputs.

Monte Carlo sampling is flexible: it works with continuous or discrete distributions and can be paired with a joint sampling procedure appropriate to the inputs. Its basic design guarantee is that each vector has the target joint distribution and vectors are independent. It does not guarantee that a finite sample includes values in every part of the input range, or that each marginal sample is evenly spaced.

Worked Example: A Small Monte Carlo Sample

Suppose a hypothetical uncertain input \(\Theta\) is Uniform(2, 10), measured in arbitrary input units, and the output of interest is \(g(\Theta)=\Theta\). Four independent Monte Carlo draws are generated using uniforms \(0.12, 0.38, 0.67,\) and \(0.91\). The transformation is \(\Theta=2+8U\).

\begin{aligned} \theta^{(1)}&=2+8(0.12)=2.96,\\ \theta^{(2)}&=2+8(0.38)=5.04,\\ \theta^{(3)}&=2+8(0.67)=7.36,\\ \theta^{(4)}&=2+8(0.91)=9.28. \end{aligned}

The Monte Carlo estimate of the mean output is the average of the four sampled outputs:

\widehat{\mu}_g=\frac{2.96+5.04+7.36+9.28}{4} =\frac{24.64}{4}=6.16.

The exact mean of a Uniform(2, 10) input is \((2+10)/2=6\), so this particular estimate is 0.16 input units above the exact mean. That difference is a possible result of random sampling; four draws are not required to balance across the range.

Latin Hypercube Sampling

LHS imposes a coverage condition on each marginal distribution. For sample size \(N\), divide the probability scale of each input into \(N\) intervals, each with probability \(1/N\). Select one probability value at random within each interval, transform those values through the corresponding marginal quantile functions, and randomly permute the values across rows for each input. Each input then has exactly one draw in each marginal probability stratum.

For independent inputs, the marginal permutations are typically chosen independently. LHS therefore improves marginal coverage compared with an unconstrained random sample of the same size. It does not ensure that all combinations of marginal strata occur: in \(d\) dimensions, there are \(N^d\) combinations of \(N\) marginal strata, but LHS has only \(N\) rows. Nor does standard independent-permutation LHS represent a specified dependence structure. A dependence-aware construction must be used when the target inputs are dependent.

Worked Example: Constructing a Two-Input Latin Hypercube

Consider two independent hypothetical inputs, \(\Theta_A\sim\mathrm{Uniform}(2,10)\) and \(\Theta_B\sim\mathrm{Uniform}(10,30)\). Use \(N=4\). The probability strata are \([0,0.25), [0.25,0.50), [0.50,0.75),\) and \([0.75,1]\). Suppose the selected probability values, already arranged into four rows by randomized permutations, are as follows.

Row	\(U_A\)	\(U_B\)	\(\Theta_A=2+8U_A\)	\(\Theta_B=10+20U_B\)
1	0.10	0.60	2.80	22.0
2	0.35	0.10	4.80	12.0
3	0.60	0.95	6.80	29.0
4	0.95	0.35	9.60	17.0

For \(\Theta_A\), the four probability values occupy strata 1 through 4 once each. For \(\Theta_B\), the values occupy strata 3, 1, 4, and 2, again once each. The marginal sample means are \((2.80+4.80+6.80+9.60)/4=6.00\) and \((22+12+29+17)/4=20.0\), equal to the respective distribution means in this particular design. The particular row pairings are randomized; LHS guarantees marginal coverage, not a specific relationship between the inputs.

Stratified Sampling and Its Weighted Estimate

Stratified sampling is a general approach rather than one particular marginal design. Partition the target distribution into mutually exclusive strata \(S_1,\ldots,S_H\) that cover the target space. Let \(W_h=P(\boldsymbol{\Theta}\in S_h)\), so that \(\sum_{h=1}^H W_h=1\). Draw \(n_h\) samples from the target distribution conditional on belonging to stratum \(S_h\). If the quantity of interest is a model output \(g(\boldsymbol{\Theta})\), the natural estimator of its mean is a weighted combination of the stratum sample means.

The stratum weights matter. If strata have unequal probabilities, taking an unweighted average of all sampled outputs generally gives too much influence to strata that were sampled more heavily than their probability warrants. LHS can be viewed as a special marginal stratification with equal-probability strata and one draw per stratum for each input; general stratified designs can instead use different partitions, multiple draws per stratum, and unequal stratum probabilities.

Result: If \(g(\boldsymbol{\Theta})\) has a finite mean, and samples are drawn from the correct conditional distribution within each stratum, then \(\widehat{\mu}_g=\sum_{h=1}^H W_h\bar{g}_h\) is an unbiased estimator of \(\mu_g=E[g(\boldsymbol{\Theta})]\). Here \(\bar{g}_h\) is the average of the \(n_h\) sampled outputs in stratum \(h\). If draws are independent within and across strata and the conditional variances are finite, then \(\mathrm{Var}(\widehat{\mu}_g)=\sum_{h=1}^H W_h^2\sigma_h^2/n_h\), where \(\sigma_h^2=\mathrm{Var}(g(\boldsymbol{\Theta})\mid\boldsymbol{\Theta}\in S_h)\).

Why the Stratified Estimator Is Unbiased

Write \(\mu_h=E[g(\boldsymbol{\Theta})\mid\boldsymbol{\Theta}\in S_h]\). Since the strata partition the target distribution, the law of total expectation gives \(\mu_g=\sum_h W_h\mu_h\). Correct conditional sampling ensures that each sampled output in stratum \(h\) has mean \(\mu_h\), so \(E[\bar{g}_h]=\mu_h\). Therefore,

E[\widehat{\mu}_g] =E\left[\sum_{h=1}^H W_h\bar{g}_h\right] =\sum_{h=1}^H W_h E[\bar{g}_h] =\sum_{h=1}^H W_h\mu_h =\mu_g.

With independent draws within each stratum, \(\mathrm{Var}(\bar{g}_h)=\sigma_h^2/n_h\). If the stratum estimates are independent, the variance of their weighted sum is consequently \(\sum_h W_h^2\mathrm{Var}(\bar{g}_h)=\sum_h W_h^2\sigma_h^2/n_h\), as stated. This result explains why a useful partition often separates regions with substantially different output behavior: it can reduce within-stratum variation. Stratification does not guarantee a reduction for every response or every partition.

Worked Example: A Weighted Estimate From Four Strata

For illustration, let \(X\sim\mathrm{Uniform}(0,1)\), and consider \(g(X)=X^2\). Partition the range into four equal-probability strata. Draw one hypothetical value from each: \(0.10, 0.40, 0.60,\) and \(0.90\). Each stratum has weight \(W_h=0.25\), so the weighted estimate is:

\begin{aligned} \widehat{\mu}_g &=0.25\left[(0.10)^2+(0.40)^2+(0.60)^2+(0.90)^2\right]\\ &=0.25(0.01+0.16+0.36+0.81)\\ &=0.25(1.34)=0.335. \end{aligned}

The exact mean is \(E[X^2]=\int_0^1 x^2\,dx=1/3\approx0.3333\), so this estimate is about 0.00167 higher. The example demonstrates the weighting and calculation, not a guarantee that stratification will be closer than Monte Carlo for every finite sample.

Choosing a Design and Checking It

Design	What it guarantees	Useful when	Important limitation
Random Monte Carlo	Independent input vectors from the target joint distribution	A flexible baseline is needed, particularly with a complex joint sampler	Finite samples can leave gaps or cluster by chance
Latin hypercube	One draw in each equal-probability stratum of every marginal	Marginal coverage is valuable and a moderate number of model evaluations is available	Marginal balance does not guarantee joint-space coverage or preserve dependence
General stratified sampling	Specified allocation across defined strata, with weighted estimation when needed	Important regions can be identified and sampled deliberately	Requires valid stratum probabilities and correct conditional sampling

Before running a PBPK model, verify that the sample matches the design you intended. For LHS, check that every marginal probability stratum contains exactly one sample. For general stratification, check stratum membership, sample counts, and weights. For any design with dependent inputs, inspect whether the sampled joint distribution reflects the intended dependence rather than merely the correct marginals. Also check that transformations produce valid values—for example, positive inputs remain positive and bounded inputs stay within their bounds.

A broken stratum assignment, an incorrect quantile transformation, or an unintended loss of dependence is a sampling failure, not evidence about the PBPK model. Conversely, apparent marginal balance is not proof that the joint sample is correct. Keep the design and the target distribution conceptually separate: the target describes the uncertainty being represented; the design describes how a finite set of input vectors is selected.

Common Mistakes and Limitations

Assuming LHS guarantees every combination. LHS stratifies each marginal, not every multidimensional cell. The correction is to check joint coverage separately or use a design tailored to the relevant joint regions.
Applying independent LHS to dependent inputs without checking the result. Independent permutations generally produce a different dependence structure from the target. Use a dependence-aware sampling construction and verify the resulting joint sample.
Ignoring weights in unequal-probability strata. An unweighted sample average can estimate the wrong target mean when the allocation differs from stratum probability. Use \(\sum_h W_h\bar{g}_h\).
Expecting stratification to improve every output summary. Its benefit depends on how the chosen strata relate to the model response. A partition helpful for one output or endpoint may be unhelpful for another.
Treating a designed sample as a guarantee of accuracy. Balanced input coverage is not a universal guarantee of precise output estimates, especially for sharp thresholds, rare regions, or strong interactions.

Monte Carlo is often the straightforward choice when the joint distribution is complicated and a reliable joint sampler is available. LHS is attractive when improved marginal coverage is desired without specifying a detailed partition of the full input space. General stratification is useful when meaningful regions can be defined and their probabilities and conditional distributions are known. These approaches can also be combined, such as using randomized sampling within carefully chosen strata. The choice should be made for the target distribution and model task, not on the assumption that one design is always superior.

Check Your Understanding

Answer the questions, then compare your reasoning with the brief solutions below.

With \(N=5\), what marginal probability-stratum condition does LHS impose on each input?
Why can independent marginal LHS fail to represent a specified dependence between two PBPK inputs?
Two strata have probabilities \(0.8\) and \(0.2\), and their sampled output means are 4 and 10. What is the weighted estimate of the overall mean?
Does an LHS sample ensure that every combination of marginal strata is present? Explain.
Name one check that is particularly important when using general stratified sampling.

Solutions. (1) Each input has one sampled value in each of five equal-probability marginal strata. (2) Independent marginal permutations generally determine pairings that need not have the target joint dependence. (3) \(0.8(4)+0.2(10)=3.2+2.0=5.2\). (4) No; LHS guarantees marginal stratum coverage, not all multidimensional combinations. (5) Check that samples are in the intended strata, that the stratum probabilities and weights are correct, and that conditional draws follow the target distribution within each stratum.