Experimental Design

Motivation

Online platforms run randomized experiments on subjects (users, sessions, orders) that arrive one at a time and must be assigned to a control or a treatment group on arrival, before later subjects are seen. The precision of a treatment-effect estimate depends on how similar the two groups are in their observed covariates. When the groups are badly balanced, the estimator's variance grows, and a large literature on covariate balancing (blocking, rerandomization, matched pairs) addresses this for offline experiments, where all covariates are known in advance.

J. Zhao and Z. Zhou (arXiv:2201.12936, Management Science, 2025) cast the online problem as an online matching problem. The cost of an assignment is the discrepancy, the weight of a minimum-cost perfect matching between the two groups, and they compare three designs against an oblivious adversary that fixes the covariate sequence in advance:

• the clairvoyant matched-pair design;
• the completely randomized design;
• their online pigeonhole design.

For one continuous covariate the pigeonhole design attains $O(T^{1/4}(\log T)^{1/2})$ against $\Theta(T^{1/2})$ for complete randomization. This mission concerns the multi-dimensional result, Theorem 9.

Setting

$T$ subjects arrive in order; $T$ is even. Subject $t$ carries a covariate vector $x_t$ with $q$ discrete coordinates, the $i$-th in a finite set $S_i\subseteq[0,1]$ of size $m_i$, and $p\ge2$ continuous coordinates in $[0,1]$. Distances are Euclidean in $\mathbb R^{q+p}$. An assignment puts each subject in control ($A$) or treatment ($B$).

The discrepancy of a split with $|A|=|B|$ is

$$d_{A,B}=\min_{\sigma:A\to B\text{ bijective}}\ \sum_{t\in A}\lVert x_t-x_{\sigma(t)}\rVert_2 ,$$

the value of the assignment problem (1) of the paper.

The completely randomized design picks $A$ uniformly among the $\binom{T}{T/2}$ subsets of size $T/2$; $\mathrm{ALG}^C(x)$ is the mean of $d_{A,B}$.

A pigeonhole design fixes a partition of the covariate space into pigeonholes. When subject $t$ arrives:

1. if one group already has $T/2$ subjects, $t$ joins the other group;
2. otherwise $t$ joins the group that has fewer subjects in its own pigeonhole, and on a tie a fair coin decides.

Theorem 9 uses the pigeonholes

$$\mathcal P(\mathcal S)=S_1\times\cdots\times S_q\times\Big\{[0,c^{1/p}T^{-\varphi}),\ \dots,\ [1-c^{1/p}T^{-\varphi},1]\Big\}^p ,$$

uniform cells of edge $c^{1/p}T^{-\varphi}$ on each continuous axis, for parameters $\varphi\in(0,1)$ and $c>1$. The expected discrepancy is the mean of $d_{A,B}$ over the design's coins.

Formalization targets

Goal: Theorem 9 (p. 26)

There are $C>0$ and $T_0$, depending only on $p,q,S_1,\dots,S_q,\varphi,c$, such that for every even $T\ge T_0$ and every admissible arrival sequence,

$$\mathbb E[d_{A,B}]\le C\Big(T^{1-\varphi}+T^{(p-1)\varphi}\,\lambda_p(T)\Big),\qquad \lambda_2(T)=(\log T)^{3/2},\quad \lambda_p(T)=1\ (p\ge3).$$

Milestones

1. Lemma EC.1: Hoeffding for a binomial, $\Pr(|H(n)-pn|\le\sqrt{n\log n})\ge1-2/n^2$.
2. Lemma 5, $p=2$: if $T=m^2$ and the subjects occupy each cell of the $m\times m$ grid exactly once, then $\mathrm{ALG}^C(x)=O(T^{1/2}\log T)$.
3. Lemma 5, $p\ge3$: on the full grid of $T=m^p$ cells, $\mathrm{ALG}^C(x)=O(T^{(p-1)/p})$.
4. Balance: every run of a pigeonhole design has $|A|=|B|=T/2$.
5. Component 1: the subjects paired within their pigeonhole before the stopping time $\tau$ contribute at most $\tfrac12c^{1/p}\sqrt p\,T^{1-\varphi}$, on every run.
6. Unmatched subjects: at most one unmatched subject per pigeonhole, $H^0_\tau+H^1_\tau$ at most the number of pigeonholes, and $|H^0_\tau-H^1_\tau|=T-\tau$.
7. Component 2: the remaining subjects contribute $O(T^{(p-1)\varphi}\lambda_p(T))$ in expectation.

Corollary 2 (the case $\varphi=1/p$) is an extra theorem after the goal.

Significance

At $\varphi=1/p$ Theorem 9 gives $O(T^{(p-1)/p})$ for $p\ge3$. This is the order $\Theta(T^{(p-1)/p})$ of the matched-pair design (Theorem 7), which sees all covariates in advance. An online design that commits to each assignment on arrival therefore loses at most a constant factor in dimension $p\ge3$, and at most a $(\log T)^{3/2}$ factor at $p=2$. For the completely randomized design no comparable upper bound is known for general sequences (p. 21). Lemma 5 is the paper's bound for complete randomization on the special family of full grids, and is of independent interest as a random-labelling matching estimate.

The results are proved on paper; none of them has a machine-checked proof, and the platform holds no item on them. A formalization would close the gaps of the printed argument:

• Lemma 5's induction for $p=2$ ends at $\sqrt T(\log T)^{3/2}$ rather than the stated $\sqrt T\log T$.
• Component 2 applies Lemma 5 to unmatched subjects, which occupy distinct cells but do not fill the grid.
• The case $q\ge1$ is argued per discrete support, although the balancing rule acts across supports.

Every statement in the mission is believed true.

Difficulty

The matched part (Component 1) is deterministic: paired subjects share a pigeonhole of diameter $O(T^{-\varphi})$. The difficulty lies in the subjects left over: the unmatched ones and the subjects after $\tau$. Their labels come from coins, but which subjects end up unmatched, and when $\tau$ occurs, depend on the same coins. So the leftover set is not a fixed set with independent fair labels, and the conditional binomial law used on p. ec22 does not hold literally.

A second obstacle is the random-matching estimate itself. Bounding the expected optimal matching between randomly labelled points of a grid at the rates $T^{1/2}\log T$ ($p=2$) and $T^{(p-1)/p}$ ($p\ge3$) is a multi-scale problem of Ajtai–Komlós–Tusnády type. A single-scale count of imbalances per cell gives only $T^{1/2}$ times a power of $T$ that is too large.

Formalization scope

• Subjects are indexed $0,\dots,T-1$; an assignment is a map to Bool (false = control). The covariate space is EuclideanSpace ℝ (Fin q ⊕ Fin p). Admissible sequences have discrete coordinates in $S_i$ and continuous ones in $[0,1]$, with $S_i\subseteq[0,1]$ (the paper's normalisation).
• Discrepancy: an infimum over bijections between the groups. An unbalanced split has the junk value $0$. The designs, the matched set and its complement are balanced on every run, so the junk is never reached.
• Cells: the printed edge $c^{1/p}T^{-\varphi}$ gives a partition only when $c^{-1/p}T^{\varphi}$ is an integer. The formalization uses $m=\lceil c^{-1/p}T^{\varphi}\rceil$ cells per axis of edge $1/m\le c^{1/p}T^{-\varphi}$, with the last cell closed. A pigeonhole is keyed by the exact discrete values and the cell index on each axis. The partition is the printed one, not a quantified family.
• Randomness: the pigeonhole design uses one coin per subject, consulted only on a tie. Its expected discrepancy is the average over all $2^T$ coin vectors, which is exactly the paper's law. The completely randomized design averages over all $\binom{T}{T/2}$ control sets.
• Asymptotics: $O(f(T))$ means $\exists C>0,\exists T_0,\forall T\ge T_0$ even, $\forall x$ admissible. $C$ and $T_0$ are chosen before $T$ and $x$ and depend only on the fixed parameters. $\log$ is the natural logarithm.
• Lemma 5's family $\mathcal X$: $T=m^p$ subjects in $[0,1]^p$, one per cell, stated as injectivity of the cell map.
• Stopping objects: $\tau$ is the first $n$ at which a group holds $T/2$ subjects. Within a pigeonhole the subjects before $\tau$ are paired in arrival order, and the last one is unmatched when their number is odd.
• Corrections to the page: the second cell printed as $[T^{-\varphi},2c^{1/p}T^{-\varphi})$ on p. ec21 is read as $[c^{1/p}T^{-\varphi},2c^{1/p}T^{-\varphi})$. The family's $x'_t\notin\mathbb X$ is read as $x_{t'}\notin\mathbb X$. The bound $H\le\frac1cT^{p\varphi}$ becomes $H\le(\prod_i m_i)\,m^p$ under the ceiling and the discrete supports.

The following formalizations would trivialize the mission and are ruled out:

• a bound on the discrepancy of one run, or of the best run, instead of the average over all coin vectors;
• constants that depend on $T$ or $x$;
• an unbalanced design that exploits the junk value $0$;
• an arbitrary partition in place of the printed grid.

A complete development needs a binomial Hoeffding bound and an exact analysis of the design's recursion (balance, alternation within a pigeonhole before $\tau$). It also needs a multi-scale bound on random bichromatic matchings of grid points. That bound is reusable for Lemma 5 and for matching estimators generally. Contributions on any milestone are welcome, as are proofs of Lemma 5 for arbitrary subsets of grid points, which is what Component 2 needs.

Selected references

• J. Zhao and Z. Zhou, Pigeonhole Design: Balancing Sequential Experiments from an Online Matching Perspective, Management Science, 2025; arXiv:2201.12936v6, 2024. https://arxiv.org/abs/2201.12936
• M. Ajtai, J. Komlós and G. Tusnády, On optimal matchings, Combinatorica 4 (1984) 259–264. https://doi.org/10.1007/BF02579135
• N. Fournier and A. Guillin, On the rate of convergence in Wasserstein distance of the empirical measure, Probability Theory and Related Fields 162 (2015) 707–738. https://doi.org/10.1007/s00440-014-0583-7
• W. Hoeffding, Probability inequalities for sums of bounded random variables, Journal of the American Statistical Association 58 (1963) 13–30. https://doi.org/10.1080/01621459.1963.10500830

Motivation

Online platforms run experiments on subjects that arrive one at a time: a user opens an app, a customer reaches a checkout page, a patient enters a trial. Each subject must be assigned to control or treatment at once, before the next one arrives, and the quality of the final estimate depends on how well the two groups resemble each other in their observed covariates. Classical sequential randomization such as Efron's biased coin (Efron 1971) balances the group sizes; Zhao and Zhou (arXiv:2201.12936, Management Science 2025) instead measure imbalance by an online matching cost and ask how small it can be made against an adversarial sequence of covariates.

Their pigeonhole design partitions the covariate space into cells and balances the two groups inside each cell. With one continuous covariate the paper shows (Lemma 3 and Corollary 1) that $\sqrt T$ uniform cells achieve an expected cost of order $O(T^{1/4}(\log T)^{1/2})$, against $\Theta(T^{1/2})$ for complete randomization (Lemma 2). This mission concerns the matching lower bound, Lemma 4: no choice of pigeonholes does better than $T^{1/4}$.

Setting

A horizon of $T$ subjects, with $T$ even, arrives at times $t=1,\dots,T$. Subject $t$ carries a covariate $x_t\in[0,1]$, and the whole sequence $x=(x_1,\dots,x_T)$ is fixed in advance by an oblivious adversary. An assignment puts each subject in the control group $A$ or the treated group $B$.

The discrepancy of the split (display (1), p. 11) is the cost of a minimum-weight perfect matching between the two groups:

$$d_{A,B}=\min_{e:A\to B\ \text{bijective}}\ \sum_{a\in A}|x_a-x_{e(a)}| .$$

A pigeonhole design (Definition 3, p. 16) takes a partition of $[0,1]$ into pigeonholes. When subject $t$ arrives:

1. if one group already holds $T/2$ subjects, $t$ joins the other group;
2. otherwise, with $n^0_{k,t-1}$ and $n^1_{k,t-1}$ the control and treated subjects already in the pigeonhole $k$ of $x_t$, subject $t$ joins control if $n^0_{k,t-1}<n^1_{k,t-1}$, treatment if $n^0_{k,t-1}>n^1_{k,t-1}$, and either group with probability $1/2$ in case of a tie.

The design's expected discrepancy $\mathbb E[d_{A,B}]$ on $x$ is the average over these fair coins.

An interval partition is given by breakpoints $0=l_0<l_1<\dots<l_K=1$, $K\ge1$; its pigeonholes are $[l_0,l_1),[l_1,l_2),\dots,[l_{K-1},l_K]$, and the cell lengths are $L_k=l_k-l_{k-1}$.

Formalization targets

Goal: Lemma 4 (p. 18)

There are an absolute constant $c>0$ and a threshold $T_0$ such that for every even $T\ge T_0$ and every interval partition of $[0,1]$, some arrival sequence $x\in[0,1]^T$ gives

$$\mathbb E\bigl[d_{A,B}\bigr]\ \ge\ c\,T^{1/4}.$$

The partition is arbitrary and the hard sequence is chosen after it; $c$ and $T_0$ do not depend on $T$ or on the partition.

Milestones (proof of Lemma 4, pp. ec15–ec16)

1. Lemma EC.2 (p. ec3), the binomial mean absolute deviation: for $X\sim\mathrm{Bin}(n,p)$,

$$\mathbb E|X-np|=\sqrt{\tfrac{2p(1-p)}{\pi}\,n}+O(n^{-1/2}).$$

2. Per-pigeonhole bound. On the alternating instance, which places $n_k/2$ consecutive pairs $(l_{k-1},l_k-\epsilon)$ in cell $k$ (the pairs $(l_{K-1},l_K)$ in the last cell), the expected discrepancy is at least $c\sum_k(L_k-\epsilon)\sqrt{n_k}$.
3. Cauchy–Schwarz display. With $n_k=L_k^2T/\sum_iL_i^2$,

$$\sum_k L_k\sqrt{n_k}=\sqrt{\Bigl(\sum_k L_k^2\Bigr)\Bigl(\sum_k n_k\Bigr)}=\sqrt{T\sum_kL_k^2}.$$

4. Jensen display. If $\sum_kL_k=1$ then $\sqrt{T\sum_kL_k^2}\ge\sqrt T\sqrt{K(1/K)^2}=\sqrt{T/K}$, with equality iff every $L_k=1/K$.
5. Few pigeonholes, $K\le C_0T^{1/2}$ for a fixed $C_0>0$: some instance gives $\mathbb E[d_{A,B}]\ge cT^{1/4}$.
6. Many pigeonholes, $K\ge C_0T^{1/2}$ for a fixed $C_0>0$: some instance gives $\mathbb E[d_{A,B}]\ge cT^{1/4}$.

Milestones 5 and 6 (with $C_0=1$) together give the goal.

Significance

Lemma 4 makes the upper bound of Lemma 3 meaningful as an optimality statement: among all pigeonhole designs on one continuous covariate, $\sqrt T$ uniform cells are optimal up to the factor $(\log T)^{1/2}$. It also locates the source of the $T^{1/4}$ rate in a tension between two kinds of randomness. Long cells let fair coins inside a cell create imbalance across the cell's length. Many cells leave one unmatched subject per cell, and these must be matched across cells. The rate cannot be improved by tuning the partition alone; doing better requires designs outside the pigeonhole class.

The paper's proof is a short sketch, and no part of it has been machine-checked. A formal proof would supply an explicit hard instance for every interval partition, including the rounding of the counts $n_k$ to even integers and the instance for the many-cell case, which the paper leaves unspecified. Lemma EC.2 is a classical result (Blyth 1980) that also has no formal proof in Mathlib.

Difficulty

The algebra of milestones 3 and 4 is elementary. The difficulty lies in the two cases.

In the few-cell case, the per-cell bound of milestone 2 must be computed from the design's actual run. That requires showing that on the alternating instance the first member of each pair always meets a tie in its cell, and that the global $T/2$ cap never alters an assignment. The matching cost must then be bounded below. On the line the optimal matching is determined by the order of the points, so a cell's internal imbalance cannot be cancelled by subjects of other cells, and this needs proof. The ideal counts $L_k^2T/\sum_iL_i^2$ are not even integers, and the rounding must not destroy the $\sqrt{T/K}$ bound.

The many-cell case is the least specified step. The paper's sketch (an odd number of subjects per cell, then a coupling with complete randomization) does not say where the unmatched subjects lie. As printed, it does not handle a partition that mixes one long interval with many short ones; such a partition must be handled by a combination of the two constructions. Lemma EC.2 needs a sharp local estimate of the binomial distribution: the error term is $O(n^{-1/2})$, which is finer than what a plain central limit theorem gives.

Formalization scope

• Covariates and time. $x:\mathrm{Fin}\,T\to\mathbb R$ with every $x_t\in[0,1]$. Subjects and cells are indexed from $0$, so the paper's subject $t$ is index $t-1$ and its cell $k\in[K]$ is index $k-1$.
• Discrepancy. Display (1) is an assignment problem, and it is encoded as the infimum over bijections between the two groups. If the groups differ in size this infimum is Lean's junk value $0$. The pigeonhole design always ends with $T/2$ subjects in each group, so this value is never reached; it could only make a lower bound harder.
• Randomness. The design reads a coin vector $\{0,1\}^T$, and coin $t$ is used only when subject $t$ meets a tie. The expected discrepancy is the average over all $2^T$ coin vectors, which is exactly the paper's law. Step 1 (the $T/2$ cap) is checked before the pigeonhole rule.
• Partitions. Lemma 4 says "any pigeonhole design", and Definition 3 allows any partition. The proof treats partitions of $[0,1]$ into $K$ intervals, and the goal is stated for those, encoded by breakpoints. Partitions into arbitrary measurable sets are not covered.
• Asymptotics. $\Omega(T^{1/4})$ is read as: there exist $c>0$ and $T_0$, chosen before $T$ and the partition, such that for every even $T\ge T_0$ and every partition there exists a hard sequence. $T^{1/4}$ is the real power. In the two cases, $K=O(T^{1/2})$ and $K=\Omega(T^{1/2})$ are read as $K\le C_0\sqrt T$ and $K\ge C_0\sqrt T$ for an arbitrary fixed $C_0>0$, on which $c$ and $T_0$ may depend. In Lemma EC.2, $p$ is fixed before the constant.
• Printed slips corrected. The proof's $L_k=l_K-l_{K-1}$ is read as $l_k-l_{k-1}$. Its "due to Lemma EC.4" is read as Lemma EC.2, because the count is binomial, not hypergeometric. In milestone 2 the per-cell constant is absolute and the gap factor is the exact width $L_k-\epsilon$; the counts are arbitrary even integers, since $L_k^2T/\sum_iL_i^2$ is not an integer in general.
• Ruled out. None of the following counts as a proof of the goal: the junk discrepancy $0$ of an unbalanced split; a bound on one run of the design, or on its best run, instead of the average over all coin vectors; constants that depend on $T$, $K$ or the partition; a hard sequence chosen before the partition.
• Infrastructure. A complete development needs: the one-dimensional optimal matching (sorted matching, or the $L^1$ distance between counting functions); a lemma computing the design's run on structured instances; and the binomial mean absolute deviation. The first and third are reusable well beyond this mission, and contributions of them are welcome. An extension of the goal to arbitrary finite partitions of $[0,1]$ is also welcome, but it is not claimed here.

Selected references

• J. Zhao and Z. Zhou, Pigeonhole Design: Balancing Sequential Experiments from an Online Matching Perspective, arXiv:2201.12936v6, 2024; Management Science, 2025. https://arxiv.org/abs/2201.12936
• C. R. Blyth, Expected absolute error of the usual estimator of the binomial parameter, The American Statistician 34(3), 155–157, 1980. https://doi.org/10.1080/00031305.1980.10483022
• B. Efron, Forcing a sequential experiment to be balanced, Biometrika 58(3), 403–417, 1971. https://doi.org/10.1093/biomet/58.3.403

Motivation

Online platforms run controlled experiments in which subjects (users, sessions, sellers) arrive one at a time and must be assigned to control or treatment on arrival, before later subjects are seen. A good design makes the two groups look alike in their observed covariates, because covariate imbalance inflates the variance of treatment-effect estimates. Classical offline remedies such as matched-pair designs need all covariates in advance, and the simplest online design, complete randomization, ignores covariates altogether.

Zhao and Zhou (arXiv:2201.12936, Management Science 2025) measure imbalance by the cost of the best one-to-one matching between the treated and the control subjects, which links sequential experimental design to online matching. They propose the pigeonhole design: partition the covariate space into cells and, inside each cell, assign each arrival to the group that is currently short there. This mission formalizes their main upper bound for one continuous covariate.

Setting

A horizon $T$, assumed even, of subjects arrives sequentially. Subject $t$ carries a covariate $x_t\in[0,1]$, and the whole sequence $x=(x_1,\dots,x_T)$ is fixed in advance by an adversary who does not see the randomization (an oblivious adversary). Each subject is assigned, at arrival and irrevocably, to the control group $A$ or the treated group $B$, and in the end $|A|=|B|=T/2$.

The discrepancy of a split is the cost of a minimum-weight perfect matching between the two groups (display (1) of the paper):

$$d_{A,B}=\min_{e:A\to B\ \text{bijective}}\ \sum_{a\in A}|x_a-x_{e(a)}| .$$

The completely randomized design chooses $A$ uniformly among the $\binom{T}{T/2}$ subsets of size $T/2$; its expected discrepancy is written $\mathrm{ALG}^{\mathsf C}(x)$.

A pigeonhole design takes a partition of $[0,1]$ into cells. When subject $t$ arrives: (1) if one group already holds $T/2$ subjects, $t$ joins the other; (2) otherwise, if the cell of $x_t$ currently holds fewer control than treated subjects, $t$ joins control, if fewer treated, $t$ joins treatment, and on a tie a fair coin decides. Its expected discrepancy is the average of $d_{A,B}$ over the coins.

The uniform pigeonholes with parameter $\eta\in(0,1)$ are the $T^{\eta}$ cells $[0,T^{-\eta}),[T^{-\eta},2T^{-\eta}),\dots,[1-T^{-\eta},1]$.

Formalization targets

Goal: Lemma 3

For every fixed $\eta\in(0,1)$ there are $C>0$ and $T_0$ such that, for every even $T\ge T_0$ and every $x\in[0,1]^T$, the pigeonhole design with uniform pigeonholes satisfies

$$\mathbb E\bigl[d_{A,B}\bigr]\ \le\ C\Bigl(T^{\frac12-\frac\eta2}+T^{\frac\eta2}(\log T)^{\frac12}\Bigr).$$

The constant may depend on $\eta$ but on nothing else. Corollary 1, the case $\eta=1/2$ giving $O(T^{1/4}(\log T)^{1/2})$, is included as a separate statement.

Milestones

In attack order:

1. Hoeffding's inequality for a binomial with deviation $\sqrt{n\log n}$ (Lemma EC.1).
2. The value of $\mathrm{ALG}^{\mathsf C}$ on the two-point instance $(0,\dots,0,1,\dots,1)$ (Lemma EC.4, with corrected constant).
3. Permutation invariance of the completely randomized design (Lemma EC.5).
4. $\mathrm{ALG}^{\mathsf C}$ is $\Theta(T^{1/2})$ in the worst case (Lemma 2).
5. Every run of a pigeonhole design is balanced (remark after Definition 3).
6. A median-crossing swap does not decrease the single-cell design's expected discrepancy (Lemma EC.7, corrected).
7. The single-cell design has expected discrepancy $O(T^{1/2})$ (Proposition EC.1).
8. At the stopping time $\tau$, each cell has at most one unmatched subject and $|H^0_\tau-H^1_\tau|=T-\tau$.
9. The subjects matched inside cells contribute $O(T^{1/2-\eta/2})$ (Component 1).
10. The unmatched subjects and the last $T-\tau$ subjects contribute $O(T^{\eta/2}(\log T)^{1/2})$ (Component 2).

Significance

Lemma 3 shows that a simple, time-invariant, covariate-aware online rule beats complete randomization by a polynomial factor: with $\eta=1/2$ the expected discrepancy is $O(T^{1/4}(\log T)^{1/2})$, against $\Theta(T^{1/2})$ for complete randomization (Lemma 2). A companion result of the paper (Lemma 4) shows that no pigeonhole design does better than $\Omega(T^{1/4})$, so the bound is tight up to the logarithm. The same two-component argument underlies the paper's multi-dimensional result.

The paper proves these statements in an e-companion; none is machine-checked. A formalization adds the following. It fixes the exact model behind each expectation, namely an average over coin vectors with the balancing rule included. It repairs the places where the printed argument does not close as written. The printed constant of Lemma EC.4 is wrong, and the stated value here is the corrected one. Lemma EC.7 is false as printed in its strict form and trivial in its weak form, and the stated version is the one the proof actually establishes. The printed proof of Lemma 2 passes through a lemma (EC.6) that is false. Component 2 relies on a conditional binomial law that does not hold literally. A complete development would certify Lemma 3 independently of these slips.

Difficulty

The natural first idea is to bound each cell separately and add the results. This handles subjects that are paired inside their own cell, but it leaves out the balancing rule. Once one group is full, every later subject is forced into the other group regardless of its covariate. These forced subjects can be far from everything they are matched to. The number of forced subjects, $T-\tau$, equals the imbalance among the cells' leftover subjects, which is random and read at a random time. Bounding its expectation requires a concentration argument that is uniform over time.

A second difficulty is the per-cell bound itself. Inside one cell the design pairs consecutive arrivals with random orientation, but the number of pairs in a cell is random and depends on the stopping time. A bound of the form $c\,\ell\sqrt n$ for a cell of length $\ell$ holding $n$ subjects must therefore survive a coin-dependent truncation. Finally, the discrepancy is a minimum over matchings, so every upper bound must exhibit a matching, and every comparison between instances (Lemma EC.7) is an exchange argument over all coin outcomes.

Formalization scope

All statements live in the namespace PigeonholeDesign.Upper, and the conventions are fixed as follows:

• Time and assignments. Subjects are indexed by Fin T, from $0$, so the paper's subject $t$ is index $t-1$, and its odd index $i$ is an even 0-based position. An assignment is Fin T → Bool, with true meaning treated.
• Discrepancy. The minimum over bijections between the two groups of the summed distances, which equals the assignment problem (1). When the groups differ in size this infimum is Lean's junk value $0$. The balance milestone shows that the designs never produce such a split.
• Randomness. The pigeonhole design is run on a coin vector in $\{0,1\}^T$, and coin $t$ is read only when subject $t$ meets a tie. Every expectation is the average over all $2^T$ coin vectors (respectively over all $\binom{T}{T/2}$ control sets), never one run or the best run.
• Partition. The printed partition has $T^{\eta}$ cells, which is not an integer in general. It is pinned as $K=\lceil T^{\eta}\rceil$ equal cells of length $1/K\le T^{-\eta}$, with the last cell closed. This coincides with the printed one when $T^\eta$ is an integer.
• Asymptotics. "$O(f(T))$" means that there exist $C>0$ and $T_0$, chosen before $T$ and $x$, depending only on $\eta$ (an absolute $C$ where the page says "absolute constant"). "$\Theta$" means both bounds. The logarithm is natural.
• Corrected statements. Lemma EC.4 uses $1/\sqrt{2\pi}$ in place of the printed $1/\sqrt\pi$. Lemma EC.7 is the non-strict inequality for the specific swap of positions $i+1$ and $j+1$. Component 2 states only its conclusion. Lemma EC.6 is not posed.

Several formalizations would make the statements easy for the wrong reason, and the statements rule them out. The junk value $0$ of an unbalanced split never applies, because the designs are balanced. A bound on a single run, or on the best run, in place of the average over all coin vectors is not the claim. Constants may not depend on $T$ or $x$. The partition is fixed, and the statements do not range over arbitrary partitions.

A complete development needs:

• finite combinatorics of the design run;
• one-dimensional matching (on the line the sorted matching is optimal);
• binomial and hypergeometric tail and mean-absolute-deviation estimates;
• a central-binomial asymptotic for Lemma EC.4.

The design definitions and the balance lemma are reusable by the paper's lower-bound and multi-dimensional missions. Contributions of any milestone, of alternative proofs that bypass the faulty printed steps, and of sharper constants are welcome.

Selected references

• J. Zhao and Z. Zhou, Pigeonhole Design: Balancing Sequential Experiments from an Online Matching Perspective, arXiv:2201.12936v6 (2024); Management Science (2025). https://arxiv.org/abs/2201.12936
• W. Hoeffding, Probability inequalities for sums of bounded random variables, Journal of the American Statistical Association 58 (1963) 13–30. https://doi.org/10.1080/01621459.1963.10500830
• C. R. Blyth, Expected absolute error of the usual estimator of the binomial parameter, The American Statistician 34 (1980) 155–157. https://doi.org/10.1080/00031305.1980.10483022
• D. Bertsimas and J. N. Tsitsiklis, Introduction to Linear Optimization, Athena Scientific (1997), Chapter 7.8 (assignment problem).

Motivation

A simulation metamodel is a cheap surrogate for an expensive stochastic simulation: after running the simulation at a few design settings, one predicts the mean response at settings never simulated. Kriging, the interpolation method of geostatistics and of the design and analysis of computer experiments (DACE; Sacks, Welch, Mitchell and Wynn 1989; Santner, Williams and Notz 2003), treats the unknown response surface as a realization of a random field and predicts by the best linear predictor. Kriging was developed for deterministic computer codes, where an observation is the exact response. A stochastic simulation returns the response plus sampling noise, and the noise may be correlated across settings when the experimenter uses common random numbers (CRN).

Ankenman, Nelson and Staum, Stochastic Kriging for Simulation Metamodeling (Proc. 2008 Winter Simulation Conference, pp. 362–370), extend kriging to this setting. The mission formalizes §2 of that paper: the form of the optimal linear predictor, called stochastic kriging, and the formula for its mean squared error, which shows how much the simulation noise costs. The source used throughout is the WSC 2008 proceedings paper, not the later Operations Research 58(2) 2010 article, whose numbering differs.

Setting

Fix a probability space $(\Omega,\mathcal F,P)$ and a set of design settings $\mathbf x$. On replication $j=1,2,\dots$ at setting $\mathbf x$ the simulation outputs (display (3), with constant trend)

$$\mathcal Y_j(\mathbf x)=\beta_0+\mathsf M(\mathbf x)+\varepsilon_j(\mathbf x).$$

The constant $\beta_0$ is the overall mean. $\mathsf M$ is a mean-zero random field, the extrinsic uncertainty imposed by the modeller. $\varepsilon_j(\mathbf x)$ is mean-zero sampling noise, the intrinsic uncertainty of the simulation; its variance may depend on $\mathbf x$, and noises at different settings may be correlated (CRN). All of these have finite second moments, and the field is uncorrelated with the noise.

An experiment design runs $n_i\ge1$ replications at each of $k$ design points $\mathbf x_1,\dots,\mathbf x_k$. The data are the sample means (4)

$$\bar{\mathcal Y}(\mathbf x_i)=\frac1{n_i}\sum_{j=1}^{n_i}\mathcal Y_j(\mathbf x_i),\qquad\bar{\mathcal Y}=(\bar{\mathcal Y}(\mathbf x_1),\dots,\bar{\mathcal Y}(\mathbf x_k))^\top,$$

and the target at a point $\mathbf x_0$ is the noise-free response $\mathsf Y(\mathbf x_0)=\beta_0+\mathsf M(\mathbf x_0)$.

Write $\Sigma_{\mathsf M}(\mathbf x,\mathbf x')=\mathrm{Cov}[\mathsf M(\mathbf x),\mathsf M(\mathbf x')]$. Then $\Sigma_{\mathsf M}$ is the $k\times k$ matrix of these covariances at the design points, $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ is the vector $(\mathrm{Cov}[\mathsf M(\mathbf x_0),\mathsf M(\mathbf x_i)])_{i}$, and $\Sigma_\varepsilon$ is the $k\times k$ covariance matrix of the averaged noises $\sum_{j=1}^{n_i}\varepsilon_j(\mathbf x_i)/n_i$. A linear predictor (5) is $\lambda_0+\lambda^\top\bar{\mathcal Y}$ with arbitrary weights $(\lambda_0,\lambda)\in\mathbb R\times\mathbb R^k$; its mean squared error is $\mathrm E[(\lambda_0+\lambda^\top\bar{\mathcal Y}-\mathsf Y(\mathbf x_0))^2]$. The optimal MSE $\mathrm{MSE}^\star$ is the infimum of this quantity over all weights.

Formalization targets

Goal: display (7)

There is a rule $\Xi$, assigning a positive definite matrix $\Xi(\Sigma_{\mathsf M},\Sigma_\varepsilon)$ to every pair of positive definite matrices, such that in every model as above with $\Sigma_{\mathsf M}\succ0$ and $\Sigma_\varepsilon\succ0$, and at every $\mathbf x_0$,

$$\begin{aligned}\mathrm{MSE}^\star&=\Sigma_{\mathsf M}(\mathbf x_0,\mathbf x_0)-\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top[\Sigma_{\mathsf M}+\Sigma_\varepsilon]^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\\&=\Big[\Sigma_{\mathsf M}(\mathbf x_0,\mathbf x_0)-\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\Sigma_{\mathsf M}^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\Big]+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\,\Xi\,\Sigma_{\mathsf M}(\mathbf x_0,\cdot).\end{aligned}$$

The bracket is the usual kriging MSE. The goal also states that $\mathrm{MSE}^\star$ strictly exceeds it whenever $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\neq0$. The matrix $\Xi$ is existential and depends only on the two covariance matrices, not on $\mathbf x_0$.

Milestone: display (6)

If $\Sigma_{\mathsf M}+\Sigma_\varepsilon\succ0$, the stochastic kriging predictor

$$\widehat{\mathsf Y}(\mathbf x_0)=\beta_0+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top[\Sigma_{\mathsf M}+\Sigma_\varepsilon]^{-1}(\bar{\mathcal Y}-\beta_0\mathbf 1_k)$$

is of the form (5), attains the smallest MSE among all predictors of that form, and is the only one that does.

Further items (not milestones)

• (9)–(10): the two-point example, with closed forms for the predictor and its MSE in terms of $\tau^2$, $r_{12}$, $r_0$, $\mathsf V$, $\rho$ and $n$.
• (14): the equicorrelated example, $\mathrm{MSE}^\star=\tau^2\big(1-k r_0^2/(1+(k-1)r+\gamma/n)\big)$.
• §3, p. 365: under the Gaussian Assumption 1, $\widehat{\mathsf Y}(\mathbf x_0)=\mathrm E[\mathsf Y(\mathbf x_0)\mid\bar{\mathcal Y}]$ almost surely.

Significance

Display (7) quantifies the price of simulation noise. With exact observations the best achievable error is the kriging MSE; with noisy observations it is larger by a positive definite form in the cross-covariance vector. Every point $\mathbf x_0$ that is correlated with the design therefore loses accuracy. The two examples make this concrete. In (10) the MSE increases with the intrinsic correlation $\rho$, which is why the paper concludes that common random numbers, a standard variance-reduction device for comparing systems, do not help prediction. Display (14) is the baseline against which the paper measures the cost of estimating the noise variance. The conditional-expectation statement shows that under Gaussian assumptions stochastic kriging is optimal among all predictors, not only linear ones.

The results are classical in substance: they are the best-linear-prediction formulas of second-order random-field theory, applied to data whose covariance is $\Sigma_{\mathsf M}+\Sigma_\varepsilon$. The paper states them without proof ("we can show"). No machine-checked version exists. Formalizing them yields a reusable, measure-theoretic best-linear-predictor theorem for random vectors defined from a model rather than postulated, together with the Loewner-order fact behind the positive definiteness of $\Xi$.

Difficulty

The algebra is short; the difficulty is in deriving the second-moment structure from the model, not assuming it. The covariance of $\bar{\mathcal Y}$ must be computed from the sample-mean definition, using bilinearity of covariance over finite sums of square-integrable variables, and shown to equal $\Sigma_{\mathsf M}+\Sigma_\varepsilon$. The cross-covariance with $\mathsf Y(\mathbf x_0)$ must be shown to equal $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$. The MSE of an arbitrary affine predictor then has to be written as a quadratic in $(\lambda_0,\lambda)$, and its infimum identified.

A first idea is to take "a random vector with mean $\beta_0\mathbf 1$ and covariance $\Sigma$" as the hypothesis. That proves a different, more abstract statement and leaves the model's content unproved.

For the second line of (7), the positive definiteness of $\Xi$ does not follow from a scalar argument. It is a matrix inequality, $(A+B)^{-1}\prec A^{-1}$ for $A,B\succ0$, and it fails without $\Sigma_\varepsilon\succ0$: when $\Sigma_\varepsilon=0$ the two lines coincide.

Formalization scope

• The model. Settings form an arbitrary type X; design points are x : Fin k → X, so the paper's index $i$ is i.val + 1. Replication $j\ge1$ is index j - 1 of a sequence ε : ℕ → X → Ω → ℝ.
• Moments and covariances. Field and noise are in $L^2$ and have mean zero. Covariances are Mathlib's covariance. The optimal MSE is an infimum over all of $\mathbb R\times\mathbb R^k$, with no unbiasedness constraint and no sign constraint.
• Hypotheses not printed in §2. Three are added and disclosed:
  1. uncorrelatedness of $\mathsf M$ and $\varepsilon$, the second-order content of Assumption 1's "independent of $\mathsf M$", without which the covariance of $\bar{\mathcal Y}$ is not $\Sigma_{\mathsf M}+\Sigma_\varepsilon$;
  2. square integrability;
  3. $\Sigma_{\mathsf M}\succ0$ and $\Sigma_\varepsilon\succ0$ in the goal. Lean's matrix inverse returns $0$ for a singular matrix, and $\Xi\succ0$ is false when $\Sigma_\varepsilon$ is singular.
• Generality kept. No Gaussianity, no independence across replications and no "no CRN" restriction is imposed on (6) and (7).
• Examples. $|r_{12}|<1$ in (10) and $r\le1$ in (14) are added, since both are correlations.
• Non-trivializing. A positive semidefinite $\Xi$, or a $\Xi$ chosen after $\mathbf x_0$, would make the goal trivial or strictly weaker; the statement fixes $\Xi$ as a function of $(\Sigma_{\mathsf M},\Sigma_\varepsilon)$ before the model and the point.

Infrastructure needed and reusable beyond this mission:

• covariance of averages of $L^2$ random variables;
• the MSE of an affine predictor as a quadratic form;
• minimization of a positive definite quadratic;
• the Loewner antitonicity of the inverse;
• for the conditional-expectation item, Gaussian conditioning via uncorrelated-hence-independent residuals.

Contributions to any of these are welcome.

Selected references

• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Proceedings of the 2008 Winter Simulation Conference, IEEE, pp. 362–370, 2008. http://www.informs-sim.org/wsc08papers/042.pdf
• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Operations Research 58(2), 371–382, 2010 (later journal version, not used here). https://doi.org/10.1287/opre.1090.0754
• J. Sacks, W. J. Welch, T. J. Mitchell, H. P. Wynn, Design and Analysis of Computer Experiments, Statistical Science 4(4), 409–423, 1989. https://doi.org/10.1214/ss/1177012413
• T. J. Santner, B. J. Williams, W. I. Notz, The Design and Analysis of Computer Experiments, Springer, 2003. https://doi.org/10.1007/978-1-4757-3799-8

Motivation

Stochastic simulation models of queues, supply chains and manufacturing systems are often too slow to run at every input setting of interest. A metamodel fitted to a moderate number of simulation runs then stands in for the simulation, for instance inside an optimization or a sensitivity analysis. Kriging, developed in geostatistics and in the design and analysis of deterministic computer experiments (DACE), predicts an unknown response surface by treating it as a realization of a Gaussian random field. Deterministic kriging interpolates its data exactly, which is wrong for stochastic simulation, where each run returns the response plus sampling noise.

Ankenman, Nelson and Staum's stochastic kriging separates the two sources of uncertainty: the extrinsic uncertainty of the unknown surface and the intrinsic uncertainty of the simulation output. Its optimal predictor requires the intrinsic noise variances, which are never known in practice and are replaced by sample variances computed from the same replications. This mission formalizes the paper's first key result: this plug-in step introduces no prediction bias.

The source is the Winter Simulation Conference 2008 proceedings version of the paper (pp. 362–370). The later Operations Research 58(2) (2010) article numbers its results differently and is not the reference here.

Setting

Fix a probability space $(\Omega,\mathcal F,P)$ and design variables $\mathbf x\in\mathbb R^d$. On replication $j=1,2,\dots$ at setting $\mathbf x$ the simulation returns

$$\mathcal Y_j(\mathbf x)=\beta_0+\mathsf M(\mathbf x)+\varepsilon_j(\mathbf x),$$

where $\beta_0\in\mathbb R$ is a constant trend, $\mathsf M$ is a mean-zero random field on $\mathbb R^d$ (the unknown surface) and $\varepsilon_j(\mathbf x)$ is the noise of replication $j$. The quantity to predict at any point $\mathbf x_0$, simulated or not, is the noise-free response $\mathsf Y(\mathbf x_0)=\beta_0+\mathsf M(\mathbf x_0)$.

An experiment design is a list of pairs $(\mathbf x_i,n_i)$, $i=1,\dots,k$, with distinct design points $\mathbf x_i$ and $n_i$ replications at $\mathbf x_i$. The data are summarized by the sample means $\bar{\mathcal Y}(\mathbf x_i)=\frac1{n_i}\sum_{j=1}^{n_i}\mathcal Y_j(\mathbf x_i)$, collected in the vector $\bar{\mathcal Y}\in\mathbb R^k$, and the sample variances

$$\mathcal S^2(\mathbf x_i)=\frac1{n_i-1}\sum_{j=1}^{n_i}\big(\mathcal Y_j(\mathbf x_i)-\bar{\mathcal Y}(\mathbf x_i)\big)^2 .$$

The extrinsic covariances are the $k\times k$ matrix $\Sigma_{\mathsf M}=(\mathrm{Cov}[\mathsf M(\mathbf x_h),\mathsf M(\mathbf x_i)])_{h,i}$ and the vector $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)=(\mathrm{Cov}[\mathsf M(\mathbf x_0),\mathsf M(\mathbf x_i)])_{i}$.

Assumption 1. $\mathsf M$ is a stationary Gaussian random field: all its finite-dimensional laws are multivariate normal with mean $0$, the covariance is $\tau^2R(\mathbf x-\mathbf x')$ with $\tau^2>0$ and $R(\mathbf 0)=1$, and the covariance matrix at finitely many distinct points is positive definite. At each design point the noises $\varepsilon_1(\mathbf x_i),\varepsilon_2(\mathbf x_i),\dots$ are i.i.d. $N(0,\mathsf V(\mathbf x_i))$ with $\mathsf V(\mathbf x_i)>0$; noises at different design points are independent (no common random numbers); and the whole noise family is independent of $\mathsf M$.

The intrinsic variance at a design point is estimated by $\widehat{\mathsf V}(\mathbf x_i)=\mathcal S^2(\mathbf x_i)$, giving $\widehat\Sigma_\varepsilon=\mathrm{Diag}\{\widehat{\mathsf V}(\mathbf x_1)/n_1,\dots,\widehat{\mathsf V}(\mathbf x_k)/n_k\}$ and the plug-in stochastic kriging predictor

$$\widehat{\widehat{\mathsf Y}}(\mathbf x_0)=\beta_0+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\big[\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon\big]^{-1}\big(\bar{\mathcal Y}-\beta_0\mathbf 1_k\big).\tag{13}$$

Formalization targets

Goal: Theorem 1 (p. 366)

Under Assumption 1, with $n_i\ge2$ replications at every design point, the prediction error is integrable and

$$\mathrm E\Big[\widehat{\widehat{\mathsf Y}}(\mathbf x_0)-\mathsf Y(\mathbf x_0)\Big]=0 .$$

Milestones (p. 365)

1. Under Assumption 1, $(\mathsf Y(\mathbf x_0),\bar{\mathcal Y}(\mathbf x_1),\dots,\bar{\mathcal Y}(\mathbf x_k))$ is multivariate normal.
2. Under Assumption 1, $\mathcal S^2(\mathbf x_i)$ has a scaled chi-squared distribution: $(n_i-1)\mathcal S^2(\mathbf x_i)/\mathsf V(\mathbf x_i)\sim\chi^2_{n_i-1}$.
3. Under Assumption 1, $\mathcal S^2(\mathbf x_i)$ is strongly consistent for $\mathsf V(\mathbf x_i)$: the sample variance of the first $m$ replications converges to $\mathsf V(\mathbf x_i)$ almost surely as $m\to\infty$.

Significance

Theorem 1 says that the practitioner's shortcut, estimating the noise variances from the replications and plugging them into the optimal predictor, keeps the predictor unbiased; the cost of not knowing the intrinsic variance is paid entirely in mean squared error. This is what makes the paper's subsequent comparison of the plug-in MSE with the known-variance MSE the right measure of the penalty for estimation.

None of these statements has a machine-checked proof, and Mathlib has neither the chi-squared law of the normal sample variance nor the independence of the normal sample mean and sample variance. A complete development produces these classical facts of normal sampling theory for the first time in Lean, together with a reusable model of noisy observations of a Gaussian random field.

Difficulty

The predictor is a nonlinear function of the data: its weight vector $[\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon]^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ is random and depends on the same replications as $\bar{\mathcal Y}$. Linearity of expectation therefore does not apply directly: the obvious argument, which treats the weights as constants, fails, and the dependence between the estimated variances, the sample means and the field has to be controlled exactly. The joint independence structure of Assumption 1 — across replications, across design points, and between noise and field — is load-bearing; a weaker pairwise or pointwise independence does not suffice. Integrability of the error must also be established, since the weights depend on random variances.

Formalization scope

Design points are indexed by Fin k (the paper's $i$ is i.val + 1) and lie in EuclideanSpace ℝ (Fin d); replication $j\ge1$ of the paper is index j - 1 of a sequence, so sums run over Finset.range. The field is M : EuclideanSpace ℝ (Fin d) → Ω → ℝ with Mathlib's IsGaussianProcess; covariances are ProbabilityTheory.covariance; stationarity is $\mathrm{Cov}[\mathsf M(\mathbf y),\mathsf M(\mathbf y')]=\tau^2R(\mathbf y-\mathbf y')$ with $R(\mathbf 0)=1$. The noise law is HasLaw … (gaussianReal 0 (V (x i))); independence is one iIndepFun over the family indexed by Fin k × ℕ plus one IndepFun between that whole family and the whole field. The variance function is ℝ≥0-valued with $\mathsf V(\mathbf x_i)>0$ assumed; replication counts satisfy $n_i\ge2$ wherever $\mathcal S^2$ appears and $n_i\ge1$ where only means appear. The chi-squared law $\chi^2_\nu$ is gammaMeasure (ν/2) (1/2).

The covariance quantities $\Sigma_{\mathsf M}$ and $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ in (13) are the true ones; only $\Sigma_\varepsilon$ is estimated, and $\widehat\Sigma_\varepsilon$ is built from the same replications as $\bar{\mathcal Y}$. A formalization in which $\widehat\Sigma_\varepsilon$ is deterministic or independent of the data, in which the matrix $\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon$ may be singular (Lean's Matrix.inv returns $0$ there, reducing the claim to $\mathrm E[\mathsf M(\mathbf x_0)]=0$), or in which the expectation is asserted without integrability, would trivialize the goal and is ruled out: positive definiteness of $\Sigma_{\mathsf M}$ at distinct design points is part of Assumption 1, and integrability is part of the conclusion. $\mathbf x_0$ may coincide with a design point.

A complete development needs: the law of linear images of independent Gaussian vectors; independence of the sample mean and the residual vector for i.i.d. normal samples (Cochran's theorem); the chi-squared law of the normal sample variance; the strong law of large numbers for the sample variance; and a measurability and boundedness argument for the random weight vector. The normal-sampling results are reusable well beyond this mission and are welcome as separate contributions.

Selected references

• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Proceedings of the 2008 Winter Simulation Conference, IEEE, pp. 362–370, 2008. http://www.informs-sim.org/wsc08papers/042.pdf
• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Operations Research 58(2):371–382, 2010. https://doi.org/10.1287/opre.1090.0754
• T. J. Santner, B. J. Williams, W. I. Notz, The Design and Analysis of Computer Experiments, Springer, 2003. https://doi.org/10.1007/978-1-4757-3799-8