---
title: "Lecture 13: Wasserstein Geometry — Probability Distributions from Samples"
subtitle: "Empirical measures, convergence rates, and the curse of dimensionality"
format:
html:
code-fold: true
code-tools: true
code-copy: true
pdf:
documentclass: scrartcl
pdf-engine: xelatex
toc: true
number-sections: true
geometry:
- margin=1in
colorlinks: true
bibliography: source/ref.bib
---
## Learning Goals
By the end of this lecture, learners should be able to:
- Define the empirical measure $\hat\mu_n = \frac{1}{n}\sum \delta_{X_i}$ and explain why it is the most direct nonparametric estimator of a distribution.
- State the plug-in estimator of the Wasserstein distance and prove its consistency via Varadarajan's theorem and the SLLN.
- Describe the Fournier--Guillin upper bounds for the expected transport cost: $n^{-1/2}$ in low dimension and $n^{-p/d}$ in high dimension.
- Explain the Weed--Bach upper and lower Wasserstein dimensions, the upper rate $n^{-1/s}$ for $s > \dim_{W,p}^*(\mu)$, and when the matching lower dimension makes the exponent sharp.
- State the one-dimensional quantile representation and the Bobkov--Ledoux distributional limit $nW_2^2 \Rightarrow$ weighted Brownian bridge.
- Construct a one-sample Wasserstein goodness-of-fit test and apply it to probability-integral-transform diagnostics for financial return forecasts.
- Compare the empirical measure with smooth-density recovery, whose minimax rate under $W_p$ in dimensions $d\ge3$ is $n^{-(s+1)/(2s+d)}$ for a density of smoothness $s$.
## The Empirical Measure and Plug-in Estimation
Suppose an unknown probability distribution $\mu$ on a metric space $(\mathcal{X}, d)$ is observed only through i.i.d. samples
$$
X_1, \ldots, X_n \sim \mu.
$$
The most direct nonparametric estimator of $\mu$ is the **empirical measure**
$$
\hat\mu_n = \frac{1}{n}\sum_{i=1}^n \delta_{X_i},
$$
where $\delta_{X_i}$ denotes a unit point mass at $X_i$. If two distributions $\mu$ and $\nu$ are observed through independent samples
$$
X_1,\ldots,X_n \sim \mu, \qquad Y_1,\ldots,Y_m \sim \nu,
$$
then the **plug-in estimator** of the Wasserstein distance is
$$
W_p(\hat\mu_n, \hat\nu_m), \qquad \hat\nu_m = \frac{1}{m}\sum_{j=1}^m \delta_{Y_j}.
$$
This simple recipe is the starting point for essentially all statistical inference in Wasserstein space. However, Wasserstein recovery is much more sensitive to dimension, tails, support geometry, smoothness, and regularization than ordinary weak convergence. Understanding these sensitivities is the main subject of this lecture.
## Consistency of Empirical Measures
The topology on $\mathcal{P}_p(\mathcal{X})$ requires more than ordinary weak convergence. Villani's Definition 6.8 introduces the correct notion [@Villani2009].
::: {#def-weak-convergence-pp .definition title="Weak convergence in $\mathcal{P}_p(\mathcal{X})$"}
Let $(\mathcal{X}, d)$ be a Polish metric space, let $p \in [1, \infty)$, let $(\mu_n)_{n\ge 1}$ be a sequence in $\mathcal{P}_p(\mathcal{X})$, and let $\mu \in \mathcal{P}_p(\mathcal{X})$. We say that $\mu_n$ *converges weakly in $\mathcal{P}_p(\mathcal{X})$* to $\mu$ if any one of the following equivalent conditions holds for some (and then for every) base point $x_0 \in \mathcal{X}$:
(i) $\mu_n \Rightarrow \mu$ and
$$
\int_{\mathcal{X}} d(x_0, x)^p\, d\mu_n(x) \longrightarrow \int_{\mathcal{X}} d(x_0, x)^p\, d\mu(x).
$$
(ii) $\mu_n \Rightarrow \mu$ and
$$
\limsup_{n\to\infty} \int_{\mathcal{X}} d(x_0, x)^p\, d\mu_n(x) \le \int_{\mathcal{X}} d(x_0, x)^p\, d\mu(x).
$$
(iii) $\mu_n \Rightarrow \mu$ and
$$
\lim_{R\to\infty} \limsup_{n\to\infty} \int_{\{x: d(x_0,x) \ge R\}} d(x_0, x)^p\, d\mu_n(x) = 0.
$$
(iv) For every continuous $\varphi: \mathcal{X} \to \mathbb{R}$ satisfying $|\varphi(x)| \le C\{1 + d(x_0,x)^p\}$ for some constant $C < \infty$,
$$
\int_{\mathcal{X}} \varphi\, d\mu_n \longrightarrow \int_{\mathcal{X}} \varphi\, d\mu.
$$
:::
::: {.callout-note title="Polish spaces" collapse="true"}
**Definition.** A **Polish space** is a separable, completely metrizable topological space — equivalently, a topological space homeomorphic to a complete separable metric space.
**Why this matters in Wasserstein geometry.** Most results in optimal transport (existence of optimal couplings, metrization of weak convergence by $W_p$, Prokhorov compactness criteria) require the underlying space to be Polish. In statistical applications the ambient space is almost always Polish: $\mathbb{R}^d$, compact Riemannian manifolds, separable Hilbert spaces, and spaces of probability measures equipped with the Wasserstein distance itself all qualify.
**Key properties:**
- **Separability** ensures that $\mathcal{P}(\mathcal{X})$ (the space of Borel probability measures) is itself separable in the weak topology.
- **Completeness** guarantees that Cauchy sequences in the Wasserstein metric converge — that is, $(\mathcal{P}_p(\mathcal{X}), W_p)$ is a complete metric space.
- Together these properties imply that every Borel probability measure on a Polish space is Radon (inner regular), so no pathologies arise in the measure-theoretic foundations.
:::
Thus weak convergence in $\mathcal{P}_p(\mathcal{X})$ is ordinary weak convergence together with control of the $p$th moment, equivalently with uniform negligibility of the $p$th tails.
::: {#thm-wp-metrizes .theorem title="$W_p$ metrizes weak convergence in $\mathcal{P}_p(\mathcal{X})$"}
Let $(\mathcal{X}, d)$ be a Polish metric space and let $p \in [1, \infty)$. If $(\mu_n)_{n\ge 1} \subset \mathcal{P}_p(\mathcal{X})$ and $\mu \in \mathcal{P}_p(\mathcal{X})$, then
$$
W_p(\mu_n, \mu) \to 0 \quad \Longleftrightarrow \quad \mu_n \text{ converges weakly in } \mathcal{P}_p(\mathcal{X}) \text{ to } \mu.
$$
:::
Applying this theorem to empirical measures gives immediate consistency. If $\mu \in \mathcal{P}_p(\mathcal{X})$ and $X_1, X_2, \ldots$ are i.i.d. with law $\mu$, then the empirical measures $\hat\mu_n$ satisfy
$$
W_p(\hat\mu_n, \mu) \to 0 \qquad \text{a.s.}
$$
The proof combines two classical results:
- **Varadarajan's theorem:** $\hat\mu_n \Rightarrow \mu$ almost surely on Polish spaces.
- **Strong law of large numbers:** For any $x_0 \in \mathcal{X}$,
$$
\int_{\mathcal{X}} d(x_0,x)^p\, d\hat\mu_n(x) = \frac{1}{n}\sum_{i=1}^n d(x_0, X_i)^p \longrightarrow \int_{\mathcal{X}} d(x_0,x)^p\, d\mu(x) \quad \text{a.s.}
$$
Hence $\hat\mu_n$ converges weakly in $\mathcal{P}_p(\mathcal{X})$ to $\mu$ almost surely, and [Theorem 13.1](#thm-wp-metrizes) gives the Wasserstein consistency.
## Convergence Rates: The Fournier--Guillin Bound
For distributions on $\mathbb{R}^d$, @FournierGuillin2015 provide general moment bounds on the expected transport cost between an empirical measure and its population law.
::: {#def-fournier-guillin .definition title="Empirical transport cost and rate function"}
Let $d \ge 1$ and let $\mathcal{P}(\mathbb{R}^d)$ denote the Borel probability measures on $\mathbb{R}^d$. For $p > 0$ and $\mu, \nu \in \mathcal{P}(\mathbb{R}^d)$, define the order-$p$ transport cost
$$
T_p(\mu, \nu) = \inf_{\pi \in \Pi(\mu,\nu)} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|x - y\|_2^p\, d\pi(x,y).
$$
When $p \ge 1$, this is exactly the $p$th power of the Wasserstein distance: $T_p(\mu,\nu) = W_p^p(\mu,\nu)$. For $q > 0$, define the $q$th Euclidean moment
$$
M_q(\mu) = \int_{\mathbb{R}^d} \|x\|_2^q\, d\mu(x).
$$
For parameters $d \ge 1$, $p > 0$, and $q > p$, define the rate function
$$
\varepsilon_{p,q,d}(n) =
\begin{cases}
n^{-1/2} + n^{-(q-p)/q}, & p > d/2 \text{ and } q \ne 2p,\\[4pt]
n^{-1/2}\log(1+n) + n^{-(q-p)/q}, & p = d/2 \text{ and } q \ne 2p,\\[4pt]
n^{-p/d} + n^{-(q-p)/q}, & 0 < p < d/2 \text{ and } q \ne dp/(d-p).
\end{cases}
$$
The excluded equalities are critical cases in which the same proof produces additional logarithmic factors.
:::
::: {#thm-fournier-guillin .theorem title="Fournier--Guillin moment bound"}
Let $d \ge 1$, let $p > 0$, and let $\mu \in \mathcal{P}(\mathbb{R}^d)$. Suppose $M_q(\mu) < \infty$ for some $q > p$, and suppose the corresponding noncritical condition in [Definition 13.2](#def-fournier-guillin) holds. Then there exists a constant $C = C(p,d,q) < \infty$ such that, for every $n \ge 1$,
$$
\mathbb{E}\bigl\{T_p(\hat\mu_n, \mu)\bigr\} \le C\, M_q(\mu)^{p/q}\, \varepsilon_{p,q,d}(n).
$$
Equivalently, for $p \ge 1$,
$$
\mathbb{E}\bigl\{W_p^p(\hat\mu_n, \mu)\bigr\} \le C\, M_q(\mu)^{p/q}\, \varepsilon_{p,q,d}(n).
$$
:::
The factor $M_q(\mu)^{p/q}$ represents the dependence on the available tail moment. If $\mu$ has sufficiently many moments --- for example, if $\mu$ is compactly supported --- one may choose $q$ so large that the tail term $n^{-(q-p)/q}$ is of smaller order than the dimensional term. The leading behavior then becomes
$$
\mathbb{E}\, W_p^p(\hat\mu_n, \mu) \lesssim
\begin{cases}
n^{-1/2}, & d < 2p,\\[4pt]
n^{-1/2}\log(1+n), & d = 2p,\\[4pt]
n^{-p/d}, & d > 2p,
\end{cases}
$$
for $p \ge 1$.
These are worst-case upper bounds over broad moment classes, not exact rates for every distribution. They can be loose for regular laws in low dimension. For example, the one-dimensional result below gives $\mathbb{E}\{W_p^p(\hat\mu_n,\mu)\}=O(n^{-p/2})$ under some conditions, which is sharper than the $O(n^{-1/2})$ bound above when $p>1$.
::: {.callout-important title="The Wasserstein curse of dimensionality"}
In high dimension ($d > 2p$), the dominant rate is $n^{-p/d}$ for the $p$th power, equivalently about $n^{-1/d}$ for $W_p$ itself. This is exponentially worse than the parametric $n^{-1/2}$ rate we are accustomed to in Euclidean statistics. It reflects the sparsity of points in high-dimensional space: for a fixed sample size $n$, typical nearest-neighbor distances scale as $n^{-1/d}$.
:::
## The Weed--Bach Dimension: Intrinsic Complexity
Beyond the Euclidean case, @WeedBach2019 provide upper and lower bounds on the empirical Wasserstein rate for compact metric spaces, characterized by two intrinsic dimensions of the target measure.
::: {#def-weed-bach-dimensions .definition title="Wasserstein dimensions"}
Let $(\mathcal{X}, d)$ be a compact metric space. For $S \subset \mathcal{X}$ and $\varepsilon > 0$, let $N_\varepsilon(S)$ be the smallest integer $m$ for which $S$ can be covered by $m$ closed balls of diameter $\varepsilon$. Define the $\varepsilon$-dimension of $S$ by
$$
d_\varepsilon(S) = \frac{\log N_\varepsilon(S)}{-\log \varepsilon}.
$$
For a Borel probability measure $\mu$ on $\mathcal{X}$ and $0 \le \tau < 1$, define the $(\varepsilon,\tau)$-covering number
$$
N_\varepsilon(\mu,\tau) = \inf\{N_\varepsilon(S): S \subset \mathcal{X} \text{ Borel and } \mu(S) \ge 1-\tau\},
$$
and the $(\varepsilon,\tau)$-dimension $d_\varepsilon(\mu,\tau) = \log N_\varepsilon(\mu,\tau) / (-\log\varepsilon)$.
For $p \in [1,\infty)$, the **upper Wasserstein dimension of order $p$** is
$$
\dim_{W,p}^*(\mu) = \inf\left\{ s \in (2p,\infty): \limsup_{\varepsilon \downarrow 0} \, d_\varepsilon\!\left(\mu, \varepsilon^{sp/(s-2p)}\right) \le s \right\},
$$
with the convention $\inf\emptyset = \infty$. The **lower Wasserstein dimension** is
$$
\dim_W^*(\mu) = \lim_{\tau \downarrow 0} \liminf_{\varepsilon \downarrow 0} d_\varepsilon(\mu,\tau).
$$
:::
::: {#thm-weed-bach .theorem title="Weed--Bach asymptotic bound"}
Let $(\mathcal{X}, d)$ be a compact Polish metric space with $\operatorname{diam}(\mathcal{X}) \le 1$, let $p \in [1,\infty)$, and let $\mu$ be a Borel probability measure on $\mathcal{X}$. Let
$$
\hat\mu_n = \frac{1}{n}\sum_{i=1}^n \delta_{X_i}, \qquad X_1,\ldots,X_n \stackrel{\text{i.i.d.}}{\sim} \mu.
$$
Then the following two bounds hold.
(i) If $s > \dim_{W,p}^*(\mu)$, then $\mathbb{E}\{W_p(\mu,\hat\mu_n)\} \lesssim n^{-1/s}$.
(ii) If $t < \dim_W^*(\mu)$, then $W_p(\mu,\hat\mu_n) \gtrsim n^{-1/t}$.
:::
To see how this reproduces the Euclidean curse of dimensionality, suppose $\mathcal{X} = B(x_0, r) \subset \mathbb{R}^d$ is a compact ball and $\mu$ has a density $f$ satisfying $0 < m \le f(x) \le M < \infty$. Then
$$
\dim_W^*(\mu) = d,
$$
so the lower bound reads $W_p(\mu,\hat\mu_n) \gtrsim n^{-1/t}$ for every $t < d$, i.e., effectively $W_p(\mu,\hat\mu_n) \gtrsim n^{-1/d}$.
::: {.callout-note title="Key insight"}
The Weed--Bach result separates the *upper* and *lower* Wasserstein dimensions. When these coincide, the rate $n^{-1/s}$ is sharp: it tells us exactly how the intrinsic complexity of the target measure, not just the ambient dimension, controls empirical Wasserstein convergence.
:::
## One-Dimensional Distributions
The one-dimensional case is special because the Wasserstein distance is an $L^p$ distance between quantile functions [@BobkovLedoux2019].
### Quantile Representation
If $F$ is the CDF of $\mu$ and $F_n$ is the empirical CDF, then
$$
W_p^p(\hat\mu_n, \mu) = \int_0^1 \bigl|F_n^{-1}(u) - F^{-1}(u)\bigr|^p\, du.
$$
For two empirical distributions on the line,
$$
W_p^p(\hat\mu_n, \hat\nu_m) = \int_0^1 \bigl|F_n^{-1}(u) - G_m^{-1}(u)\bigr|^p\, du.
$$
This identity is the basis for both computation and asymptotic theory.
### Asymptotic Distribution
For a distribution $\mu$ with positive density $f$, define
$$
J_p(\mu) = \int_{-\infty}^\infty \frac{\{F(x)(1-F(x))\}^{p/2}}{f(x)^{p-1}}\, dx.
$$
For $p \ge 1$, @BobkovLedoux2019 [Theorem 5.3] show that $J_p(\mu) < \infty$ implies the standard one-dimensional bound
$$
\mathbb{E}\{W_p^p(\hat\mu_n, \mu)\} = O(n^{-p/2}).
$$
For $p = 2$, under additional regularity conditions, Theorem 5.2 of @BobkovLedoux2019 gives the distributional limit
$$
n W_2^2(\hat\mu_n, \mu) \Rightarrow \int_0^1 \frac{B(u)^2}{f(F^{-1}(u))^2}\, du,
$$
where $B$ is a standard Brownian bridge. This result is the Wasserstein analogue of the classical Kolmogorov--Smirnov and Cramér--von Mises theory: it provides a distributional limit that can be used for confidence intervals and hypothesis tests about the distribution.
::: {.callout-tip title="Why one dimension is special"}
The quantile representation shows that in one dimension, $(\mathcal{P}_2(\mathbb{R}), W_2)$ is isometric to a convex subset of $L^2(0,1)$ via $\mu \mapsto F^{-1}_\mu$. This linear structure is what makes the one-dimensional theory so much richer: empirical quantile processes are well-studied objects, and classical results about the Brownian bridge can be brought to bear directly.
:::
## Application: Backtesting Distributional Return Forecasts
### From a hypothesized forecast model to a uniformity test
Suppose a risk model issues a complete one-step-ahead predictive distribution for each out-of-sample return or profit-and-loss $R_t$. Let $\mathcal F_{t-1}$ contain all information available when the forecast for time $t$ is issued, and write $G_t$ for the model's predictive CDF. For each $r$, the forecast $G_t(r)$ is $\mathcal F_{t-1}$-measurable. The backtesting question is whether this hypothesized model supplies the correct conditional distribution:
$$
H_0^{\mathrm{model}}:\qquad
\Pr(R_t\le r\mid\mathcal F_{t-1})=G_t(r)
$$
for every out-of-sample date $t$ and every $r\in\mathbb R$, almost surely. Equivalently,
$$
R_t\mid\mathcal F_{t-1}\sim G_t
$$
at each forecast date. Thus the null does not require the predictive distribution to remain fixed over time. It requires the possibly time-varying $G_t$ to equal the true conditional distribution given the information available when the forecast was made.
This time variation prevents us from pooling the raw returns and comparing their empirical distribution with one fixed return distribution: even under a correct model, the conditional location, scale, and shape may change from one forecast date to the next. The **probability integral transform (PIT)** removes this variation by placing every realization on its own forecast's probability scale. Define
$$
U_t=G_t(R_t).
$$
To see what the transformation accomplishes, condition on $\mathcal F_{t-1}$. Conditional on that information, $G_t$ is a fixed CDF and, under $H_0^{\mathrm{model}}$, $R_t$ is a draw from it. For $0<u<1$, define the upper conditional quantile
$$
q_t(u)=\sup\{r:G_t(r)\le u\}.
$$
Continuity gives $G_t(q_t(u))=u$, and the events $\{G_t(R_t)\le u\}$ and $\{R_t\le q_t(u)\}$ agree up to a conditional probability-zero set. Therefore,
$$
\Pr(U_t\le u\mid\mathcal F_{t-1})
=\Pr\{G_t(R_t)\le u\mid\mathcal F_{t-1}\}
=\Pr\{R_t\le q_t(u)\mid\mathcal F_{t-1}\}
=G_t(q_t(u))
=u,\qquad 0<u<1.
$$
The endpoint cases follow by limits. Hence $U_t$ is uniform even after conditioning on everything known before time $t$.
This conditional statement also explains why the PIT values are **independent**, not merely marginally uniform. Let $A_1,\ldots,A_k$ be Borel subsets of $[0,1]$, and let $\lambda$ denote Lebesgue measure. Because $U_1,\ldots,U_{k-1}$ are determined by past observations and forecasts, the event $\{U_1\in A_1,\ldots,U_{k-1}\in A_{k-1}\}$ belongs to $\mathcal F_{k-1}$. The law of iterated expectations gives
$$
\begin{aligned}
&\Pr(U_1\in A_1,\ldots,U_k\in A_k)\\
&\quad=
\mathbb E\!\left[
\mathbf 1_{\{U_1\in A_1,\ldots,U_{k-1}\in A_{k-1}\}}
\Pr(U_k\in A_k\mid\mathcal F_{k-1})
\right]\\
&\quad=
\lambda(A_k)\Pr(U_1\in A_1,\ldots,U_{k-1}\in A_{k-1}).
\end{aligned}
$$
Iterating this argument yields
$$
\Pr(U_1\in A_1,\ldots,U_k\in A_k)
=\prod_{j=1}^k\lambda(A_j).
$$
Thus no separate independence assumption on the returns is needed: correct specification of each full conditional distribution implies
$$
H_0^{\mathrm{PIT}}:\qquad
U_t\stackrel{\mathrm{i.i.d.}}{\sim}\operatorname{Unif}(0,1).
$$
The PIT therefore converts a sequence of time-varying conditional return distributions into a fixed, distribution-free null target and is a standard tool for evaluating density forecasts of asset returns [@DieboldGuntherTay1998]. Continuity matters here. If a predictive distribution has atoms, the ordinary PIT need not be uniform; a randomized PIT is then required to recover an exact uniform null.
### Finite-sample formula and decision rule
We can now use the empirical-distribution asymptotics developed above. Let $\hat\lambda_n=n^{-1}\sum_{t=1}^n\delta_{U_t}$ be the empirical PIT distribution. The statistic derived above is
$$
T_n^{\mathrm{PIT}}
=nW_2^2\!\left(\hat\lambda_n,\operatorname{Unif}(0,1)\right),
$$
with asymptotic null distribution $Q:=\int_0^1 B(u)^2\,du$.
If $U_{(1)}\le\cdots\le U_{(n)}$ are the ordered PIT values, the statistic has the exact finite-sample formula
$$
T_n^{\mathrm{PIT}}
=\frac{1}{12n}
+\sum_{i=1}^n
\left(U_{(i)}-\frac{2i-1}{2n}\right)^2.
$$
Thus the one-dimensional Wasserstein test is exactly the classical Cramér--von Mises uniformity statistic, now interpreted as squared transport distance between the empirical calibration distribution and its uniform target.
For a level-$\alpha$ test, let $c_{1-\alpha}$ be the $(1-\alpha)$-quantile of $Q$. Reject the i.i.d.-uniform PIT null when
$$
T_n^{\mathrm{PIT}}>c_{1-\alpha}.
$$
The critical value can be obtained by simulating Brownian bridges. For moderate $n$, a finite-sample Monte Carlo calibration is equally simple: repeatedly generate $n$ independent uniforms, recompute $T_n^{\mathrm{PIT}}$, and compare the observed statistic with this simulated null distribution. A rejection is evidence against conditional distributional calibration. A non-rejection, however, does not establish serial independence; that component requires a separate diagnostic because $T_n^{\mathrm{PIT}}$ ignores time ordering.
### Interactive demo: a distributional return backtest
The demo below generates one reproducible evaluation sample from the standardized return model
$$
R_t\stackrel{\mathrm{i.i.d.}}{\sim}\mathcal N(0,1),
$$
but evaluates those returns using the predictive distribution $G_t=\mathcal N(m_f,s_f^2)$. It then computes $U_t=G_t(R_t)$ and applies the Wasserstein--Cramér--von Mises test. The forecast is correctly calibrated when $m_f=0$ and $s_f=1$. A nonzero $m_f$ creates directional bias; $s_f<1$ makes the predictive distribution too narrow; and $s_f>1$ makes it too wide.
The PIT histogram, uniform QQ plot, and observed statistic all use the **same single simulated return sample**. For each choice of $n$ and $\alpha$, a separate batch of 999 samples from the uniform null is used only to estimate the critical value and Monte Carlo $p$-value. Move the seed control to generate a different evaluation sample.
```{ojs}
//| echo: false
//| output: false
pit13_n_control = Inputs.range([25, 1000], {step: 25, value: 250, label: "Evaluation sample size n"})
pit13_mean_control = Inputs.range([-0.8, 0.8], {step: 0.05, value: 0, label: "Forecast mean m_f"})
pit13_scale_control = Inputs.range([0.5, 1.8], {step: 0.05, value: 1, label: "Forecast scale s_f"})
pit13_alpha_control = Inputs.select([0.01, 0.05, 0.1], {value: 0.05, format: x => `α = ${x}`, label: "Test level"})
pit13_seed_control = Inputs.range([1, 100], {step: 1, value: 13, label: "Sample seed"})
pit13_n = Generators.input(pit13_n_control)
pit13_mean = Generators.input(pit13_mean_control)
pit13_scale = Generators.input(pit13_scale_control)
pit13_alpha = Generators.input(pit13_alpha_control)
pit13_seed = Generators.input(pit13_seed_control)
pit13_controls = html`
<style>
.pit13-control-grid { display:grid; grid-template-columns:repeat(3,minmax(0,1fr)); gap:6px 20px; margin:0 auto 14px; font:0.85em system-ui,sans-serif; }
.pit13-control-grid > * { min-width:0; margin:0; }
.pit13-control-grid input[type="number"] { width:7.5rem !important; }
@media (max-width:700px) { .pit13-control-grid { grid-template-columns:1fr; } }
</style>
<div class="pit13-control-grid">
<div>${pit13_n_control}</div>
<div>${pit13_mean_control}</div>
<div>${pit13_scale_control}</div>
<div>${pit13_alpha_control}</div>
<div>${pit13_seed_control}</div>
</div>`
function pit13Rng(seed) {
let state = seed | 0;
return function() {
state = state + 0x6D2B79F5 | 0;
let t = Math.imul(state ^ state >>> 15, 1 | state);
t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
return ((t ^ t >>> 14) >>> 0) / 4294967296;
};
}
function pit13Normal(rng) {
const u1 = Math.max(rng(), Number.EPSILON);
return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * rng());
}
function pit13NormalCdf(x) {
const z = Math.abs(x);
const t = 1 / (1 + 0.2316419 * z);
const poly = t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429))));
const upper = Math.exp(-0.5 * z * z) * poly / Math.sqrt(2 * Math.PI);
return x >= 0 ? 1 - upper : upper;
}
function pit13Statistic(sortedPits) {
const n = sortedPits.length;
let statistic = 1 / (12 * n);
for (let i = 0; i < n; i++) {
const target = (2 * (i + 1) - 1) / (2 * n);
statistic += (sortedPits[i] - target) ** 2;
}
return statistic;
}
function pit13Simulate(n, forecastMean, forecastScale, seed) {
const rng = pit13Rng(seed);
const returns = Array.from({length: n}, () => pit13Normal(rng));
const pits = returns.map(r => pit13NormalCdf((r - forecastMean) / forecastScale));
const sortedPits = pits.slice().sort((a, b) => a - b);
return {returns, pits, sortedPits, statistic: pit13Statistic(sortedPits)};
}
function pit13NullCalibration(n, alpha, repetitions, seed) {
const rng = pit13Rng(seed);
const statistics = new Array(repetitions);
for (let b = 0; b < repetitions; b++) {
const uniforms = Array.from({length: n}, () => rng()).sort((a, b) => a - b);
statistics[b] = pit13Statistic(uniforms);
}
statistics.sort((a, b) => a - b);
const criticalIndex = Math.min(repetitions - 1, Math.ceil((1 - alpha) * repetitions) - 1);
return {statistics, criticalValue: statistics[criticalIndex]};
}
pit13_result = pit13Simulate(pit13_n, pit13_mean, pit13_scale, 1009 * pit13_seed + 17)
pit13_null = pit13NullCalibration(pit13_n, pit13_alpha, 999, 20260313 + 37 * pit13_n)
pit13_exceedances = pit13_null.statistics.filter(x => x >= pit13_result.statistic).length
pit13_pvalue = (pit13_exceedances + 1) / (pit13_null.statistics.length + 1)
pit13_reject = pit13_result.statistic > pit13_null.criticalValue
```
```{ojs}
//| label: fig-pit-return-backtest
//| fig-cap: "Interactive Wasserstein backtest for a Gaussian return forecast"
//| echo: false
html`
<div style="font-family:system-ui,sans-serif; max-width:900px; margin:0 auto;">
${pit13_controls}
<div style="display:grid; grid-template-columns:repeat(auto-fit,minmax(150px,1fr)); gap:10px; margin:0 0 14px;">
<div style="padding:10px 12px; background:#f8f9fa; border-radius:6px;">
<div style="font-size:0.78em; color:#6c757d;">Observed statistic</div>
<div style="font-size:1.25em; font-weight:600;">${pit13_result.statistic.toFixed(3)}</div>
</div>
<div style="padding:10px 12px; background:#f8f9fa; border-radius:6px;">
<div style="font-size:0.78em; color:#6c757d;">MC critical value</div>
<div style="font-size:1.25em; font-weight:600;">${pit13_null.criticalValue.toFixed(3)}</div>
</div>
<div style="padding:10px 12px; background:#f8f9fa; border-radius:6px;">
<div style="font-size:0.78em; color:#6c757d;">MC p-value</div>
<div style="font-size:1.25em; font-weight:600;">${pit13_pvalue.toFixed(3)}</div>
</div>
<div style="padding:10px 12px; color:${pit13_reject ? "#9c1c1c" : "#176b2c"}; background:${pit13_reject ? "#fff0f0" : "#ebfbee"}; border-radius:6px;">
<div style="font-size:0.78em;">Decision at α = ${pit13_alpha}</div>
<div style="font-size:1.1em; font-weight:700;">${pit13_reject ? "Reject calibration" : "Do not reject"}</div>
</div>
</div>
<div style="display:flex; gap:16px; flex-wrap:wrap;">
<div style="flex:1 1 390px; min-width:280px;">
<h4 style="margin:0 0 6px;">PIT histogram</h4>
<svg viewBox="0 0 420 255" role="img" aria-label="Histogram of probability integral transform values" style="width:100%; height:auto; border:1px solid #dee2e6; border-radius:4px;">
${(() => {
const M = {top:18, right:15, bottom:38, left:45};
const W = 420 - M.left - M.right, H = 255 - M.top - M.bottom;
const bins = 10, counts = Array(bins).fill(0);
pit13_result.pits.forEach(u => counts[Math.min(bins - 1, Math.floor(u * bins))]++);
const expected = pit13_n / bins;
const yMax = Math.max(expected * 1.35, ...counts) * 1.08;
const xS = u => M.left + u * W;
const yS = y => M.top + H - y / yMax * H;
let svg = "";
counts.forEach((count, j) => {
const x = xS(j / bins) + 1;
const width = W / bins - 2;
svg += `<rect x="${x}" y="${yS(count)}" width="${width}" height="${M.top + H - yS(count)}" fill="#4c78a8" opacity="0.82"/>`;
});
svg += `<line x1="${M.left}" y1="${yS(expected)}" x2="${M.left + W}" y2="${yS(expected)}" stroke="#e8590c" stroke-width="2" stroke-dasharray="5 4"/>`;
svg += `<text x="${M.left + W - 3}" y="${yS(expected) - 5}" text-anchor="end" font-size="9" fill="#e8590c">uniform expectation</text>`;
[0, 0.25, 0.5, 0.75, 1].forEach(u => {
svg += `<text x="${xS(u)}" y="${M.top + H + 17}" text-anchor="middle" font-size="9">${u.toFixed(2).replace(/0+$/, "").replace(/\.$/, "")}</text>`;
});
[0, 0.5, 1].forEach(frac => {
const value = frac * yMax;
svg += `<text x="${M.left - 7}" y="${yS(value) + 3}" text-anchor="end" font-size="9">${Math.round(value)}</text>`;
});
svg += `<text x="${M.left + W / 2}" y="${M.top + H + 32}" text-anchor="middle" font-size="10">PIT value</text>`;
svg += `<text x="12" y="${M.top + H / 2}" text-anchor="middle" font-size="10" transform="rotate(-90 12 ${M.top + H / 2})">Count</text>`;
svg += `<rect x="${M.left}" y="${M.top}" width="${W}" height="${H}" fill="none" stroke="#adb5bd"/>`;
return svg;
})()}
</svg>
</div>
<div style="flex:1 1 390px; min-width:280px;">
<h4 style="margin:0 0 6px;">Uniform QQ plot</h4>
<svg viewBox="0 0 420 255" role="img" aria-label="Quantile-quantile plot of PIT values against the uniform distribution" style="width:100%; height:auto; border:1px solid #dee2e6; border-radius:4px;">
${(() => {
const M = {top:18, right:15, bottom:38, left:45};
const W = 420 - M.left - M.right, H = 255 - M.top - M.bottom;
const xS = u => M.left + u * W;
const yS = u => M.top + H - u * H;
const stride = Math.max(1, Math.floor(pit13_n / 250));
let svg = `<line x1="${xS(0)}" y1="${yS(0)}" x2="${xS(1)}" y2="${yS(1)}" stroke="#e8590c" stroke-width="2"/>`;
for (let i = 0; i < pit13_n; i += stride) {
const expected = (i + 0.5) / pit13_n;
svg += `<circle cx="${xS(expected)}" cy="${yS(pit13_result.sortedPits[i])}" r="2.1" fill="#4c78a8" opacity="0.72"/>`;
}
[0, 0.25, 0.5, 0.75, 1].forEach(u => {
const label = u.toFixed(2).replace(/0+$/, "").replace(/\.$/, "");
svg += `<text x="${xS(u)}" y="${M.top + H + 17}" text-anchor="middle" font-size="9">${label}</text>`;
svg += `<text x="${M.left - 7}" y="${yS(u) + 3}" text-anchor="end" font-size="9">${label}</text>`;
});
svg += `<text x="${M.left + W / 2}" y="${M.top + H + 32}" text-anchor="middle" font-size="10">Uniform quantile</text>`;
svg += `<text x="12" y="${M.top + H / 2}" text-anchor="middle" font-size="10" transform="rotate(-90 12 ${M.top + H / 2})">Empirical PIT quantile</text>`;
svg += `<rect x="${M.left}" y="${M.top}" width="${W}" height="${H}" fill="none" stroke="#adb5bd"/>`;
return svg;
})()}
</svg>
</div>
</div>
<p style="margin:10px 2px 0; font-size:0.84em; color:#6c757d;">
Returns are simulated once from N(0,1). The orange histogram line is the expected bin count under uniform calibration; the orange QQ line is the 45° reference.
</p>
</div>`
```
::: {.callout-tip title="Try these calibration failures"}
- Set $s_f=0.65$: the forecast is too narrow, so PIT values accumulate near 0 and 1.
- Set $s_f=1.5$: the forecast is too wide, so PIT values accumulate near $1/2$.
- Set $m_f=0.4$: the forecast is centered above the true return distribution, shifting PIT values toward 0.
- Return to $m_f=0$, $s_f=1$, then change the seed: at level $\alpha$, a correctly calibrated forecast can still be rejected in roughly an $\alpha$ fraction of repeated samples.
:::
## Smooth-Density Recovery
The empirical measure is not always statistically optimal. If the target distribution is known to have a smooth density, this smoothness can be exploited to achieve faster convergence rates.
For densities on $[0,1]^d$ belonging to smooth Besov-type classes and bounded away from zero, @NilesWeedBerthet2022 show that the minimax rate for estimating the distribution under $W_p$ is
$$
n^{-(s+1)/(2s+d)}
$$
in dimensions $d \ge 3$, where $s$ represents the degree of smoothness of the density. In dimension $d=1$ the minimax rate is $n^{-1/2}$. In dimension $d=2$, their upper bound is $n^{-1/2}\log n$ and their lower bound is $n^{-1/2}$, leaving a logarithmic gap.
::: {.callout-important title="Smoothness improves rates dramatically"}
For $s = 2$ (twice-differentiable densities) in dimension $d = 3$, the rate is $n^{-3/7} \approx n^{-0.43}$, compared to the high-dimensional empirical-measure benchmark $n^{-1/3} \approx n^{-0.33}$. As $d$ grows, the smooth estimator's exponent is asymptotically $(s+1)$ times the empirical exponent $1/d$, although both rates deteriorate with dimension.
The gain requires genuine structural assumptions. The estimator in Theorem 1 of @NilesWeedBerthet2022 uses the density's smoothness level. Their Theorem 2 gives an estimator adaptive to smoothness, with an additional logarithmic factor when $d\ge2$ (and no rate penalty when $d=1$).
:::
## Estimating Wasserstein Distance Between Two Distributions
For two unknown distributions $\mu, \nu$, the plug-in estimator $W_p(\hat\mu_n, \hat\nu_m)$ is consistent whenever both empirical measures are consistent. The triangle inequality gives the basic deterministic error bound
$$
\bigl|W_p(\hat\mu_n, \hat\nu_m) - W_p(\mu, \nu)\bigr| \le W_p(\hat\mu_n, \mu) + W_p(\hat\nu_m, \nu).
$$
Therefore all one-sample recovery rates immediately imply two-sample rates for estimating $W_p(\mu, \nu)$.
## Interactive Exploration: Empirical Wasserstein Convergence
The following one-dimensional demo simulates empirical measures from one of four population distributions:
$$
\mathcal N(0,1), \qquad \operatorname{Exp}(1), \qquad
\tfrac12\mathcal N(-2,1)+\tfrac12\mathcal N(2,1), \qquad
\operatorname{Unif}(0,1).
$$
For the selected sample size $n$, the demo draws one reproducible sample using a fixed pseudorandom seed and forms its empirical measure $\hat\mu_n$. The **Empirical CDF vs. true CDF** panel compares the step CDF of this single sample (red) with the known population CDF evaluated on a fine grid (blue). The **QQ plot** also uses this same single sample: at probability levels $u=0.01,\ldots,0.99$, it plots the empirical quantile $\hat F_n^{-1}(u)$ against the corresponding population quantile $F^{-1}(u)$. Thus these two panels are single-realization diagnostics; they are not based on a Monte Carlo average or a Wasserstein barycenter.
The **Monte Carlo mean of $W_2$** panel answers a different question. At each displayed sample size, the program draws the requested number $R$ of independently seeded samples, computes $W_2(\hat\mu_n,\mu)$ for every sample, and plots their arithmetic mean. The vertical bars are approximate 95% Monte Carlo confidence intervals for this mean,
$$
\overline W_2 \pm 1.96\,\frac{S_{W_2}}{\sqrt R}.
$$
Increasing the Monte Carlo repetitions therefore reduces simulation noise in the estimated mean curve and its error bars; it does not average the empirical distributions or alter the single sample used by the CDF and QQ panels. The summary labeled **MC mean at this $n$** is computed from a separate Monte Carlo batch at the currently selected $n$. No Wasserstein barycenter is constructed anywhere in the demo.
The orange curve is an $n^{-1/2}$ reference slope, the standard one-dimensional rate under the regularity condition $J_2(\mu)<\infty$; it is not asserted to be the exact rate for every distribution in the menu.
```{ojs}
//| echo: false
//| output: false
n13_control = Inputs.range([10, 2000], {step: 10, value: 200, label: "Sample size n"})
dist_type_control = Inputs.select(["normal", "exponential", "bimodal", "uniform"], {value: "normal", label: "Distribution type"})
num_reps_control = Inputs.range([5, 100], {step: 5, value: 30, label: "Monte Carlo repetitions"})
n13 = Generators.input(n13_control)
dist_type = Generators.input(dist_type_control)
num_reps = Generators.input(num_reps_control)
empw_controls_view = html`
<style>
.empw-slider-grid { display:grid; grid-template-columns:repeat(3,minmax(0,1fr)); gap:6px 20px; max-width:850px; margin:0 auto 14px; font:0.85em system-ui,sans-serif; }
.empw-slider-grid > * { min-width:0; margin:0; }
.empw-slider-grid input[type="number"] { width:7.5rem !important; }
@media (max-width:700px) { .empw-slider-grid { grid-template-columns:1fr; } }
</style>
<div class="empw-slider-grid">
<div>${n13_control}</div>
<div>${dist_type_control}</div>
<div>${num_reps_control}</div>
</div>`
function normalCdf(x) {
const z = Math.abs(x);
const t = 1 / (1 + 0.2316419 * z);
const poly = t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429))));
const upperTail = Math.exp(-0.5 * z * z) * poly / Math.sqrt(2 * Math.PI);
return x >= 0 ? 1 - upperTail : upperTail;
}
function normalQuantile(u) {
const t = u < 0.5 ? u : 1 - u;
const s = Math.sqrt(-2 * Math.log(t));
const x = s - (2.515517 + 0.802853 * s + 0.010328 * s * s) / (1 + 1.432788 * s + 0.189269 * s * s + 0.001308 * s * s * s);
return u < 0.5 ? -x : x;
}
function trueCdf(x, dist) {
if (dist === "normal") return normalCdf(x);
if (dist === "exponential") return x < 0 ? 0 : 1 - Math.exp(-x);
if (dist === "bimodal") return 0.5 * normalCdf(x + 2) + 0.5 * normalCdf(x - 2);
return Math.max(0, Math.min(1, x));
}
function bimodalQuantile(u) {
let lo = -10, hi = 10;
for (let i = 0; i < 60; i++) {
const mid = (lo + hi) / 2;
if (trueCdf(mid, "bimodal") < u) lo = mid;
else hi = mid;
}
return (lo + hi) / 2;
}
function trueQuantile(u, dist) {
if (dist === "normal") return normalQuantile(u);
if (dist === "exponential") return -Math.log(1 - u);
if (dist === "bimodal") return bimodalQuantile(u);
return u;
}
function empiricalQuantile(samples, u) {
if (u <= 0) return samples[0];
if (u >= 1) return samples[samples.length - 1];
return samples[Math.ceil(u * samples.length) - 1];
}
w2IntegrationPoints = 2000
trueQuantileGrid = Object.fromEntries(
["normal", "exponential", "bimodal", "uniform"].map(function(dist) {
return [dist, Array.from({length: w2IntegrationPoints}, function(_, k) {
return trueQuantile((k + 0.5) / w2IntegrationPoints, dist);
})];
})
)
function empW2(n, dist, seed) {
let state = seed || 42;
function rng() {
state |= 0;
state = state + 0x6D2B79F5 | 0;
let t = Math.imul(state ^ state >>> 15, 1 | state);
t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
return ((t ^ t >>> 14) >>> 0) / 4294967296;
}
function standardNormal() {
const u1 = Math.max(rng(), Number.EPSILON);
return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * rng());
}
function oneSample() {
if (dist === "normal") return standardNormal();
if (dist === "exponential") return -Math.log(1 - rng());
if (dist === "bimodal") return standardNormal() + (rng() < 0.5 ? -2 : 2);
return rng();
}
const samples = Array.from({length: n}, oneSample).sort(function(a, b) { return a - b; });
let gridMin, gridMax;
if (dist === "uniform") {
gridMin = -0.1; gridMax = 1.1;
} else if (dist === "exponential") {
gridMin = -0.2; gridMax = Math.max(6, samples[n - 1] + 0.25);
} else {
const limit = dist === "bimodal" ? 5 : 4;
gridMin = Math.min(-limit, samples[0] - 0.25);
gridMax = Math.max(limit, samples[n - 1] + 0.25);
}
const targetQuantiles = trueQuantileGrid[dist];
let w2Sq = 0;
for (let k = 0; k < w2IntegrationPoints; k++) {
const u = (k + 0.5) / w2IntegrationPoints;
const delta = empiricalQuantile(samples, u) - targetQuantiles[k];
w2Sq += delta * delta;
}
return {w2: Math.sqrt(w2Sq / w2IntegrationPoints), gridMin, gridMax, samples, dist};
}
function mcW2(n, dist, reps, seed) {
let sum = 0, sumSq = 0;
for (let r = 0; r < reps; r++) {
const w = empW2(n, dist, (seed || 42) + r * 1000).w2;
sum += w;
sumSq += w * w;
}
const mean = sum / reps;
const variance = reps > 1 ? Math.max(0, (sumSq - reps * mean * mean) / (reps - 1)) : 0;
return {mean, std: Math.sqrt(variance)};
}
l13result = empW2(n13, dist_type, 42)
mcNs = [10, 25, 50, 100, 200, 400, 750, 1000, 1500, 2000]
mcResults = mcNs.map(function(sampleSize) { return mcW2(sampleSize, dist_type, num_reps, 1729 + sampleSize); })
l13CurrentMc = mcW2(n13, dist_type, num_reps, 104729)
```
```{ojs}
//| label: fig-empirical-wasserstein
//| fig-cap: "Interactive empirical Wasserstein convergence demo"
//| echo: false
html`
<div style="font-family:system-ui,sans-serif; max-width:900px; margin:0 auto;">
${empw_controls_view}
<div style="display:flex; gap:16px; flex-wrap:wrap;">
<div style="flex:1 1 380px; min-width:280px;">
<h4 style="margin:0 0 6px;">Empirical CDF vs. true CDF</h4>
<svg viewBox="0 0 420 240" role="img" aria-label="Empirical and true cumulative distribution functions" style="width:100%; height:auto; border:1px solid #dee2e6; border-radius:4px;">
${(() => {
const M = {top: 15, right: 15, bottom: 34, left: 45};
const W = 420 - M.left - M.right, H = 240 - M.top - M.bottom;
const samples = l13result.samples, sampleSize = samples.length;
const gMin = l13result.gridMin, gMax = l13result.gridMax;
function xS(v) { return M.left + ((v - gMin) / (gMax - gMin)) * W; }
function yS(v) { return M.top + H - v * H; }
let svg = '';
for (let i = 0; i < 5; i++) {
const v = gMin + i * (gMax - gMin) / 4;
svg += '<text x="' + xS(v) + '" y="' + (M.top + H + 17) + '" text-anchor="middle" font-size="9">' + v.toFixed(1) + '</text>';
}
[0, 0.25, 0.5, 0.75, 1].forEach(function(v) {
svg += '<text x="' + (M.left - 7) + '" y="' + (yS(v) + 3) + '" text-anchor="end" font-size="9">' + v.toFixed(2).replace(/0+$/, '').replace(/\.$/, '') + '</text>';
});
svg += '<text x="' + (M.left + W / 2) + '" y="' + (M.top + H + 31) + '" text-anchor="middle" font-size="10">Value</text>';
let cdfPath = '';
for (let i = 0; i <= 250; i++) {
const x = gMin + i * (gMax - gMin) / 250;
cdfPath += (i === 0 ? 'M' : 'L') + xS(x).toFixed(2) + ' ' + yS(trueCdf(x, l13result.dist)).toFixed(2);
}
svg += '<path d="' + cdfPath + '" fill="none" stroke="#1971c2" stroke-width="2.2"/>';
let ecdfPath = 'M' + xS(gMin).toFixed(2) + ' ' + yS(0).toFixed(2);
samples.forEach(function(x, i) {
ecdfPath += 'H' + xS(x).toFixed(2) + 'V' + yS((i + 1) / sampleSize).toFixed(2);
});
ecdfPath += 'H' + xS(gMax).toFixed(2);
svg += '<path d="' + ecdfPath + '" fill="none" stroke="#e03131" stroke-width="1.6" opacity="0.72"/>';
svg += '<text x="' + (M.left + 10) + '" y="' + (M.top + 14) + '" font-size="10" fill="#1971c2">True CDF</text>';
svg += '<text x="' + (M.left + 10) + '" y="' + (M.top + 29) + '" font-size="10" fill="#e03131">ECDF</text>';
svg += '<rect x="' + M.left + '" y="' + M.top + '" width="' + W + '" height="' + H + '" fill="none" stroke="#adb5bd"/>';
return svg;
})()}
</svg>
</div>
<div style="flex:1 1 380px; min-width:280px;">
<h4 style="margin:0 0 6px;">Monte Carlo mean of W₂ (log–log scale)</h4>
<svg viewBox="0 0 420 240" role="img" aria-label="Monte Carlo mean Wasserstein distance by sample size" style="width:100%; height:auto; border:1px solid #dee2e6; border-radius:4px;">
${(() => {
const M = {top: 18, right: 22, bottom: 34, left: 55};
const W = 420 - M.left - M.right, H = 240 - M.top - M.bottom;
const ns = mcNs;
const means = mcResults.map(function(result) { return result.mean; });
const ses = mcResults.map(function(result) { return result.std / Math.sqrt(num_reps); });
const referenceScale = means[0] * Math.sqrt(ns[0]);
const references = ns.map(function(sampleSize) { return referenceScale / Math.sqrt(sampleSize); });
const lowerValues = means.map(function(mean, i) { return Math.max(mean - 1.96 * ses[i], mean * 0.5); }).concat(references);
const upperValues = means.map(function(mean, i) { return mean + 1.96 * ses[i]; }).concat(references);
const xMin = ns[0], xMax = ns[ns.length - 1];
const yMin = Math.min.apply(null, lowerValues) * 0.8;
const yMax = Math.max.apply(null, upperValues) * 1.25;
function xS(v) { return M.left + (Math.log(v) - Math.log(xMin)) / (Math.log(xMax) - Math.log(xMin)) * W; }
function yS(v) { return M.top + H - (Math.log(v) - Math.log(yMin)) / (Math.log(yMax) - Math.log(yMin)) * H; }
let svg = '';
[10, 50, 200, 1000, 2000].forEach(function(v) {
svg += '<text x="' + xS(v) + '" y="' + (M.top + H + 17) + '" text-anchor="middle" font-size="9">' + v + '</text>';
});
for (let i = 0; i < 4; i++) {
const v = yMin * Math.pow(yMax / yMin, i / 3);
svg += '<text x="' + (M.left - 7) + '" y="' + (yS(v) + 3) + '" text-anchor="end" font-size="9">' + v.toPrecision(2) + '</text>';
}
svg += '<text x="' + (M.left + W / 2) + '" y="' + (M.top + H + 31) + '" text-anchor="middle" font-size="10">Sample size n</text>';
svg += '<text x="13" y="' + (M.top + H / 2) + '" text-anchor="middle" font-size="10" transform="rotate(-90,13,' + (M.top + H / 2) + ')">Mean W₂</text>';
let meanPath = '', referencePath = '';
ns.forEach(function(sampleSize, i) {
meanPath += (i === 0 ? 'M' : 'L') + xS(sampleSize).toFixed(2) + ' ' + yS(means[i]).toFixed(2);
referencePath += (i === 0 ? 'M' : 'L') + xS(sampleSize).toFixed(2) + ' ' + yS(references[i]).toFixed(2);
});
svg += '<path d="' + meanPath + '" fill="none" stroke="#1971c2" stroke-width="2.4"/>';
svg += '<path d="' + referencePath + '" fill="none" stroke="#e8590c" stroke-width="1.7" stroke-dasharray="6,4"/>';
ns.forEach(function(sampleSize, i) {
const low = Math.max(yMin, means[i] - 1.96 * ses[i]);
const high = Math.min(yMax, means[i] + 1.96 * ses[i]);
svg += '<line x1="' + xS(sampleSize) + '" y1="' + yS(high) + '" x2="' + xS(sampleSize) + '" y2="' + yS(low) + '" stroke="#868e96" stroke-width="1"/>';
svg += '<circle cx="' + xS(sampleSize) + '" cy="' + yS(means[i]) + '" r="3" fill="#1971c2"/>';
});
svg += '<text x="' + xS(150) + '" y="' + (yS(referenceScale / Math.sqrt(150)) - 7) + '" font-size="9" fill="#e8590c">n⁻¹ᐟ² reference</text>';
svg += '<rect x="' + M.left + '" y="' + M.top + '" width="' + W + '" height="' + H + '" fill="none" stroke="#adb5bd"/>';
return svg;
})()}
</svg>
</div>
</div>
<div style="margin-top:14px; max-width:700px;">
<h4 style="margin:0 0 6px;">QQ plot</h4>
<svg viewBox="0 0 700 230" role="img" aria-label="Empirical versus theoretical quantiles" style="width:100%; height:auto; border:1px solid #dee2e6; border-radius:4px;">
${(() => {
const M = {top: 14, right: 16, bottom: 38, left: 58};
const W = 700 - M.left - M.right, H = 230 - M.top - M.bottom;
const samples = l13result.samples;
const qq = [];
for (let i = 1; i < 100; i++) {
const u = i / 100;
qq.push({emp: empiricalQuantile(samples, u), theo: trueQuantile(u, l13result.dist)});
}
const allVals = qq.flatMap(function(point) { return [point.emp, point.theo]; });
const qMin = Math.min.apply(null, allVals), qMax = Math.max.apply(null, allVals);
const pad = Math.max(0.1, (qMax - qMin) * 0.08);
const x0 = qMin - pad, x1 = qMax + pad;
function xS(v) { return M.left + ((v - x0) / (x1 - x0)) * W; }
function yS(v) { return M.top + H - ((v - x0) / (x1 - x0)) * H; }
let svg = '<line x1="' + xS(x0) + '" y1="' + yS(x0) + '" x2="' + xS(x1) + '" y2="' + yS(x1) + '" stroke="#868e96" stroke-width="1.5" stroke-dasharray="6,3"/>';
for (let i = 0; i < 5; i++) {
const v = x0 + i * (x1 - x0) / 4;
svg += '<text x="' + xS(v) + '" y="' + (M.top + H + 18) + '" text-anchor="middle" font-size="9">' + v.toFixed(1) + '</text>';
svg += '<text x="' + (M.left - 7) + '" y="' + (yS(v) + 3) + '" text-anchor="end" font-size="9">' + v.toFixed(1) + '</text>';
}
svg += '<text x="' + (M.left + W / 2) + '" y="' + (M.top + H + 34) + '" text-anchor="middle" font-size="10">Theoretical quantiles</text>';
svg += '<text x="14" y="' + (M.top + H / 2) + '" text-anchor="middle" font-size="10" transform="rotate(-90,14,' + (M.top + H / 2) + ')">Empirical quantiles</text>';
qq.forEach(function(point) {
svg += '<circle cx="' + xS(point.theo) + '" cy="' + yS(point.emp) + '" r="2.6" fill="#1971c2" opacity="0.62"/>';
});
svg += '<rect x="' + M.left + '" y="' + M.top + '" width="' + W + '" height="' + H + '" fill="none" stroke="#adb5bd"/>';
return svg;
})()}
</svg>
</div>
<div style="margin-top:12px; padding:10px 14px; background:#f8f9fa; border-radius:6px; display:flex; gap:20px; flex-wrap:wrap;">
<div><span style="color:#495057;">n =</span> <strong>${l13result.samples.length}</strong></div>
<div><span style="color:#495057;">Distribution:</span> <strong>${dist_type}</strong></div>
<div><span style="color:#495057;">Current W₂:</span> <strong>${l13result.w2.toFixed(4)}</strong></div>
<div><span style="color:#495057;">MC mean at this n:</span> <strong>${l13CurrentMc.mean.toFixed(4)}</strong></div>
</div>
<p style="margin:10px 0 0; font-size:0.85em; color:#495057;">
Error bars show approximate 95% Monte Carlo confidence intervals for the simulated mean. The computed W₂ uses the empirical measure's step quantile function and numerical integration over 2,000 midpoint quantiles.
</p>
</div>
`
```
::: {.callout-tip title="Try these experiments"}
- **Increase $n$:** Watch the empirical CDF approach the true density, the $W_2$ distance decay, and the QQ plot align with the diagonal.
- **Change the distribution type:** Compare convergence across normal, exponential, bimodal, and uniform distributions.
- **Increase Monte Carlo repetitions:** Get smoother estimates of the expected $W_2$ curve.
- **Observe the bimodal case:** The QQ plot reveals the two-component mixture structure.
:::
## Key Takeaways
- The **empirical measure** $\hat\mu_n = \frac{1}{n}\sum_{i=1}^n \delta_{X_i}$ is the most direct nonparametric estimator of a distribution, and the **plug-in estimator** $W_p(\hat\mu_n, \hat\nu_m)$ provides a natural estimate of the Wasserstein distance.
- **Consistency** follows from Varadarajan's theorem (weak convergence) and the strong law of large numbers (moment convergence), together with Villani's theorem that $W_p$ metrizes weak convergence in $\mathcal{P}_p$.
- The **Fournier--Guillin bound** reveals the Wasserstein curse of dimensionality: $\mathbb{E}[W_p^p] \lesssim n^{-1/2}$ when $d < 2p$, but $\lesssim n^{-p/d}$ when $d > 2p$ — a catastrophic slowdown in high dimension.
- The **Weed--Bach dimension** provides an intrinsic characterization via the upper Wasserstein dimension $\dim_{W,p}^*(\mu)$, giving almost-tight upper and lower bounds $n^{-1/s}$.
- In **one dimension**, the quantile representation gives a complete theory: $W_p$ equals an $L^p$ distance between quantile functions, and $nW_2^2$ converges to a weighted Brownian bridge functional.
- For **smooth densities**, the minimax rate $n^{-(s+1)/(2s+d)}$ can be much faster than the empirical measure rate, especially when $d$ is moderate.
- **Two-sample estimation** reduces to the one-sample case via the triangle inequality.
## Exercises
1. **Consistency from Varadarajan + SLLN.** Complete the proof that $W_p(\hat\mu_n, \mu) \to 0$ almost surely when $\mu \in \mathcal{P}_p(\mathcal{X})$ on a Polish space $\mathcal{X}$. Specify exactly where Varadarajan's theorem and the strong law are used, and explain why convergence of the $p$th moment is needed. <a href="javascript:void(0)" onclick="showSolution('l13-sol-1')" class="solution-link">Show Solution</a>
2. **The curse of dimensionality in the Fournier--Guillin bound.** Let $\mu$ be the uniform distribution on $[0,1]^d$. For $p=2$, compute the leading Fournier--Guillin upper bound on $\mathbb{E}[W_2^2(\hat\mu_n,\mu)]$ for $d=1,2,5,10$. Using only the powers of $n$ and ignoring unknown constants, how much larger is the $d=10$ rate proxy than the low-dimensional $n^{-1/2}$ proxy at $n=10^6$? Why should these bounds not be read as exact rates for every $\mu$? <a href="javascript:void(0)" onclick="showSolution('l13-sol-2')" class="solution-link">Show Solution</a>
3. **Weed--Bach dimension of a uniform distribution on $[0,1]^d$.** Show that for $\mu = \operatorname{Unif}([0,1]^d)$ with the Euclidean metric, $\dim_W^*(\mu) = d$. (Hint: bound the covering number $N_\varepsilon([0,1]^d)$ from above and below using volume arguments, then compute $d_\varepsilon(\mu,\tau)$.) <a href="javascript:void(0)" onclick="showSolution('l13-sol-3')" class="solution-link">Show Solution</a>
4. **Quantile representation for $W_1$ on the line.** Let $\mu,\nu\in\mathcal P_1(\mathbb R)$, with CDFs $F,G$ and quantile functions $F^{-1},G^{-1}$. Prove that $W_1(\mu,\nu) = \int_0^1 |F^{-1}(u) - G^{-1}(u)|\, du = \int_{-\infty}^\infty |F(x) - G(x)|\, dx$. (The second equality shows that $W_1$ equals the $L^1$ distance between CDFs, a useful fact for two-sample testing.) <a href="javascript:void(0)" onclick="showSolution('l13-sol-4')" class="solution-link">Show Solution</a>
5. **Minimax rate for smooth densities.** Compare the empirical measure rate $n^{-1/d}$ with the Niles-Weed--Berthet minimax rate $n^{-(s+1)/(2s+d)}$ for $s = 2$ (twice-differentiable densities). Compute the exponent for $d = 3, 5, 10$. For which dimensions is the improvement from using smooth-density estimation most dramatic? <a href="javascript:void(0)" onclick="showSolution('l13-sol-5')" class="solution-link">Show Solution</a>
<style>
.solution-link { font-size: 0.9em; text-decoration: none; white-space: nowrap; margin-left: 0.3em; }
.solution-link:hover { text-decoration: underline; }
.solution-dialog { padding: 0; max-width: 720px; }
.solution-dialog-header { display: flex; justify-content: space-between; align-items: flex-start; border-bottom: 1px solid #dee2e6; padding: 1.25rem 1.5rem 1rem; background: #f8f9fa; border-radius: 8px 8px 0 0; }
.solution-dialog-header h4 { margin: 0; font-size: 1.15rem; }
.solution-dialog-close { background: none; border: 1px solid #adb5bd; border-radius: 4px; padding: 0.2rem 0.75rem; cursor: pointer; font-size: 0.9rem; color: #495057; white-space: nowrap; flex-shrink: 0; }
.solution-dialog-close:hover { background: #e9ecef; }
.solution-original { padding: 1rem 1.5rem; background: #f1f3f5; border-left: 4px solid #868e96; margin: 1rem 1.5rem; border-radius: 4px; font-size: 0.95rem; }
.solution-original strong { color: #495057; }
.solution-answer { padding: 0.5rem 1.5rem 1.5rem; }
.solution-answer strong { color: #2b8a3e; }
dialog { border: none; border-radius: 8px; box-shadow: 0 8px 32px rgba(0,0,0,0.22); padding: 0; max-width: 750px; width: 90vw; }
dialog::backdrop { background: rgba(0,0,0,0.45); }
</style>
<dialog id="l13-sol-1">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 1: Consistency Proof</h4>
<button onclick="closeSolution('l13-sol-1')" class="solution-dialog-close">Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Complete the proof that $W_p(\hat\mu_n, \mu) \to 0$ almost surely, specifying where Varadarajan's theorem and the SLLN are used.</div>
<div class="solution-answer">
<strong>Solution:</strong>
We need to show that $\hat\mu_n$ converges weakly in $\mathcal{P}_p(\mathcal{X})$ to $\mu$, then invoke Villani's Theorem 6.9.
**Step 1: Weak convergence.** Varadarajan's theorem states that for an i.i.d. sample from a Borel probability measure $\mu$ on a Polish space $\mathcal{X}$, the empirical measures $\hat\mu_n = \frac{1}{n}\sum_{i=1}^n \delta_{X_i}$ satisfy $\hat\mu_n \Rightarrow \mu$ almost surely. This gives the weak-convergence part of condition (i) in [Definition 13.1](#def-weak-convergence-pp).
**Step 2: Moment convergence.** Fix any $x_0 \in \mathcal{X}$. By the strong law of large numbers,
$$
\frac{1}{n} \sum_{i=1}^n d(x_0, X_i)^p \longrightarrow \mathbb{E}[d(x_0, X)^p] = \int_{\mathcal{X}} d(x_0,x)^p\, d\mu(x) \quad \text{a.s.},
$$
since $\mu \in \mathcal{P}_p(\mathcal{X})$ guarantees the expectation is finite. But the left-hand side is exactly $\int d(x_0,x)^p\, d\hat\mu_n(x)$. Hence condition (i) of [Definition 13.1](#def-weak-convergence-pp) holds almost surely.
**Step 3: Apply Villani's theorem.** Since $\hat\mu_n$ converges weakly in $\mathcal{P}_p(\mathcal{X})$ to $\mu$ almost surely, [Theorem 13.1](#thm-wp-metrizes) (Villani's Theorem 6.9) gives $W_p(\hat\mu_n, \mu) \to 0$ almost surely.
**Why moment convergence is necessary:** Without control of the $p$th moment, weak convergence alone does not guarantee $W_p$ convergence. For example, a sequence of measures that puts small mass at increasingly distant points can converge weakly to a standard normal while the Wasserstein distance remains bounded away from zero. The moment condition prevents mass from escaping to infinity, ensuring that the transport cost is controlled.
</div>
</div>
</dialog>
<dialog id="l13-sol-2">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 2: Curse of Dimensionality</h4>
<button onclick="closeSolution('l13-sol-2')" class="solution-dialog-close">Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Compute the leading Fournier--Guillin bound on $\mathbb{E}[W_2^2(\hat\mu_n,\mu)]$ for $d=1,2,5,10$.</div>
<div class="solution-answer">
<strong>Solution:</strong>
Because the uniform law is compactly supported, its tail term can be made negligible. With $p=2$, the threshold is $2p=4$, so
$$
\mathbb{E}[W_2^2(\hat\mu_n,\mu)] \lesssim
\begin{cases}
n^{-1/2}, & d=1,2,\\
n^{-2/5}, & d=5,\\
n^{-1/5}, & d=10.
\end{cases}
$$
At $n=10^6$, the ratio of the $d=10$ and low-dimensional rate proxies is
$$
\frac{n^{-1/5}}{n^{-1/2}}=n^{3/10}=(10^6)^{3/10}=10^{1.8}\approx 63.1.
$$
This compares only powers of $n$; the theorem's unknown constants need not agree across dimensions. Moreover, these are general upper bounds, not equalities. For instance, the uniform law on $[0,1]$ satisfies the one-dimensional regularity condition and in fact has $\mathbb{E}[W_2^2]=O(n^{-1})$, sharper than the general $O(n^{-1/2})$ Fournier--Guillin bound.
</div>
</div>
</dialog>
<dialog id="l13-sol-3">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 3: Weed--Bach Dimension</h4>
<button onclick="closeSolution('l13-sol-3')" class="solution-dialog-close">Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Show that for $\mu = \operatorname{Unif}([0,1]^d)$, $\dim_W^*(\mu) = d$.</div>
<div class="solution-answer">
<strong>Solution:</strong>
**Upper bound on covering numbers.** The unit cube $[0,1]^d$ can be covered by $N_\varepsilon([0,1]^d) \le C_d \varepsilon^{-d}$ sets of diameter $\varepsilon$. For an explicit construction, partition each coordinate axis into $\lceil\sqrt d/\varepsilon\rceil$ intervals. The resulting subcubes have side length at most $\varepsilon/\sqrt d$ and hence diameter at most $\varepsilon$. Thus $N_\varepsilon([0,1]^d)\le \lceil\sqrt d/\varepsilon\rceil^d\lesssim_d\varepsilon^{-d}$, so $d_\varepsilon(\mu,\tau)\le d+o(1)$ for fixed $\tau<1$.
**Lower bound on covering numbers.** For any set $S \subset [0,1]^d$ with $\mu(S) \ge 1-\tau$, we have $\operatorname{Leb}(S) \ge 1-\tau$. Cover $S$ by balls of diameter $\varepsilon$. By a volume argument, at least $(1-\tau)/(C'_d \varepsilon^d)$ balls are needed. Hence $N_\varepsilon(\mu,\tau) \gtrsim_d (1-\tau)\varepsilon^{-d}$, and $d_\varepsilon(\mu,\tau) \ge d + o(1)$ for each fixed $\tau<1$.
**Conclusion.** For any $\tau < 1$, $d_\varepsilon(\mu,\tau) \to d$ as $\varepsilon \to 0$. Hence
$$
\lim_{\tau \downarrow 0} \liminf_{\varepsilon \downarrow 0} d_\varepsilon(\mu,\tau) = d,
$$
so $\dim_W^*(\mu) = d$. By [Theorem 13.3](#thm-weed-bach), this gives $W_p(\mu,\hat\mu_n) \gtrsim n^{-1/t}$ for every $t < d$, i.e., essentially $n^{-1/d}$.
This confirms that the lower Wasserstein dimension of a uniform distribution on $[0,1]^d$ equals the Euclidean dimension, matching our intuition from the Fournier--Guillin bound.
</div>
</div>
</dialog>
<dialog id="l13-sol-4">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 4: Quantile Representation for $W_1$</h4>
<button onclick="closeSolution('l13-sol-4')" class="solution-dialog-close">Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Prove that $W_1(\mu,\nu) = \int_0^1 |F^{-1}(u)-G^{-1}(u)|\, du = \int_{-\infty}^\infty |F(x)-G(x)|\, dx$.</div>
<div class="solution-answer">
<strong>Solution:</strong>
**First equality ($W_1$ as $L^1$ of quantiles).** In one dimension, the optimal transport map for $W_1$ is given by the quantile matching: the optimal coupling is the comonotone coupling, which pairs quantiles in order. Under this coupling,
$$
W_1(\mu,\nu) = \int_0^1 |F^{-1}(u) - G^{-1}(u)|\, du.
$$
This follows from the general identity $W_p(\mu,\nu) = (\int_0^1 |F^{-1}(u) - G^{-1}(u)|^p\, du)^{1/p}$ for $p \ge 1$ (set $p=1$).
**Second equality ($L^1$ of CDFs).** Let $U\sim\operatorname{Unif}(0,1)$ and set $X=F^{-1}(U)$ and $Y=G^{-1}(U)$. For any real $a,b$,
$$
|a-b|=\int_{\mathbb R}\left|\mathbf 1\{a\le x\}-\mathbf 1\{b\le x\}\right|\,dx.
$$
Apply this identity to $X,Y$ and use Tonelli's theorem. For each $x$, the generalized-inverse property gives $\mathbf 1\{F^{-1}(U)\le x\}=\mathbf 1\{U\le F(x)\}$ almost surely, and similarly for $G$. Because the same $U$ appears in both indicators,
$$
\mathbb E\left|\mathbf 1\{U\le F(x)\}-\mathbf 1\{U\le G(x)\}\right|=|F(x)-G(x)|.
$$
Therefore
$$
\int_{-\infty}^\infty |F(x) - G(x)|\, dx = \int_0^1 |F^{-1}(u) - G^{-1}(u)|\, du,
$$
so both expressions equal $W_1(\mu,\nu)$.
**Interpretation:** This result is practically useful because it shows that $W_1$ can be computed as the area between the two CDFs, which is a simple numerical integral — no optimization over couplings needed.
</div>
</div>
</dialog>
<dialog id="l13-sol-5">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 5: Smooth-Density Improvement</h4>
<button onclick="closeSolution('l13-sol-5')" class="solution-dialog-close">Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Compare $n^{-1/d}$ with $n^{-(s+1)/(2s+d)}$ for $s=2$, $d=3,5,10$.</div>
<div class="solution-answer">
<strong>Solution:</strong>
For $s=2$ (twice-differentiable densities), the Niles-Weed--Berthet rate is $n^{-(s+1)/(2s+d)} = n^{-3/(4+d)}$.
**Comparison of exponents:**
| $d$ | Empirical: $-1/d$ | Smooth: $-3/(4+d)$ | Smooth/empirical proxy at $n=10^6$ |
|---|---|---|---|
| 3 | $-1/3\approx-0.333$ | $-3/7\approx-0.429$ | smooth $\times$ 0.27 |
| 5 | $-1/5=-0.200$ | $-1/3\approx-0.333$ | smooth $\times$ 0.16 |
| 10 | $-1/10=-0.100$ | $-3/14\approx-0.214$ | smooth $\times$ 0.21 |
**Key observations:**
1. **Dimension $d=3$:** The exponent improves from $1/3$ to $3/7$. At $n=10^6$, comparing only the powers of $n$, the smooth-rate proxy is about $10^{-4/7}\approx0.27$ times the empirical-rate proxy.
2. **Dimension $d=5$:** The empirical rate $n^{-0.2}$ is quite slow, while the smooth rate is $n^{-1/3}$. At $n=10^6$, the smooth-rate proxy is about $10^{-0.8}\approx0.16$ times the empirical-rate proxy.
3. **Dimension $d=10$:** The empirical rate $n^{-0.1}$ needs a factor of $10^{10}$ in sample size to reduce its rate proxy by ten. The smooth rate $n^{-3/14}$ needs a factor of about $10^{14/3}\approx4.6\times10^4$.
4. **The real benefit:** Smoothness changes the exponent from $-1/d$ to $-(s+1)/(2s+d)$. In high dimensions, $1/d \approx 0$ while $(s+1)/(2s+d) \approx (s+1)/d$, a $(s+1)$-fold improvement in the exponent. For $s=2$, this is a 3-fold improvement, which translates to enormous savings in sample size.
Among the three fixed-$n$ proxies in the table, the largest absolute exponent gain occurs at $d=5$, whose ratio is smallest. In relative terms, however, the exponent gain grows with dimension: at $d=10$ the smooth exponent is more than twice the empirical exponent, and the ratio tends to $s+1=3$ as $d\to\infty$.
These comparisons use only rate powers and ignore different unknown constants. Smooth-density estimation also requires the density to belong to the assumed smoothness class and to be bounded away from zero. The estimator achieving the displayed rate uses the smoothness level; Theorem 2 of @NilesWeedBerthet2022 provides adaptation, with an extra logarithmic factor when $d\ge2$.
</div>
</div>
</dialog>
<script>
function showSolution(id) {
const dialog = document.getElementById(id);
if (dialog) {
dialog.showModal();
if (window.MathJax && MathJax.typesetPromise) {
MathJax.typesetPromise([dialog]).catch(function(err) { console.log('MathJax typeset error:', err); });
}
}
}
function closeSolution(id) {
const dialog = document.getElementById(id);
if (dialog) { dialog.close(); }
}
document.addEventListener('click', function(e) {
if (e.target.tagName === 'DIALOG') { e.target.close(); }
});
</script>
## Further Reading
- @FournierGuillin2015 — The foundational paper on convergence rates of empirical measures in Wasserstein distance, giving general moment bounds in $\mathbb{R}^d$.
- @WeedBach2019 — Upper and lower bounds for empirical Wasserstein convergence in compact metric spaces via intrinsic Wasserstein dimensions.
- @BobkovLedoux2019 — Comprehensive treatment of one-dimensional Wasserstein distance, including quantile representations, distributional limits, and the Brownian bridge approximation.
- @DieboldGuntherTay1998 — Probability-integral-transform diagnostics for density forecasts, with an application to financial risk management.
- @Villani2009 — The standard reference for optimal transport theory; Theorem 6.9 provides the key result that $W_p$ metrizes weak convergence in $\mathcal{P}_p$.
- @NilesWeedBerthet2022 — Minimax rates for smooth density estimation under Wasserstein loss, showing how smoothness can break the curse of dimensionality.
## Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
👉 **[Lecture 13 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-13-quiz.qmd)**