---
title: "Lecture 14: Wasserstein Geometry — Distribution-on-Distribution Regression"
subtitle: "Regression when both predictor and response are probability distributions"
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:
- Survey the regression literature for distributional data, distinguishing transformation-based, Fréchet, and optimal-transport-based approaches.
- Formulate the distribution-on-distribution regression problem in terms of a monotone optimal-transport map.
- Derive the Fréchet least-squares criterion for transport-map regression and explain its strict convexity.
- Describe the weighted isotonic regression (PAVA) estimator for the regression map and state its convergence rate.
- Apply distribution-on-distribution regression to economic and financial contexts where both predictor and response are distributions.
## Overview: Two Approaches
### From Scalar Regression to Distributional Regression
Classical regression models a scalar response $Y \in \mathbb{R}$ as a function of a predictor $X \in \mathbb{R}^p$ plus noise: $Y = f(X) + \varepsilon$. When the observations are **probability distributions** rather than scalars or vectors, this paradigm must be rethought. The regression literature for distributional data has grown rapidly over the past decade, driven by applications in which the natural unit of observation is an entire histogram, density, or empirical measure.
Broadly, the literature splits into two strategies:
| Approach | Key idea | Representative work |
|---|---|---|
| **Transformation-based / extrinsic** | Map one-dimensional densities to a linear function space, apply FDA or time-series methods, and map valid results back | @PetersenMueller2016; @KokoszkaEtAl2019 |
| **Wasserstein-based** | Use Wasserstein Fréchet means, tangent coordinates, or optimal-transport maps | @PetersenMueller2019; @PetersenLiuDivani2021; @ChenLinMueller2023; @GhodratiPanaretos2022 |
### Transformation-Based Approaches
The transformation strategy addresses a fundamental challenge: probability densities do not form a vector space (they are constrained to be nonnegative and integrate to one), so linear operations on densities — addition, scalar multiplication, inner products — are not naturally defined.
@PetersenMueller2016 introduced a framework that circumvents this obstacle by mapping suitably regular one-dimensional densities to $L^2$ through transformations equipped with inverse maps. Two transformations studied in detail are:
1. **Log quantile density transform (LQDT).** For a positive density $f$ on $[0,1]$, define the quantile density $q(u) = Q'(u)=f\{Q(u)\}^{-1}$ and set $\psi_Q(f)(u) = \log q(u)$. The paper defines a normalized inverse that maps transformed functions back to densities on $[0,1]$. This supports tools such as FPCA and functional regression without treating the density itself as an unconstrained $L^2$ function.
2. **Log hazard transform.** For densities on $[0,1]$, the paper uses $\psi_H(f)(t)=\log[f(t)/\{1-F(t)\}]$ on a truncated interval $[0,1-\delta]$ because the hazard diverges at the right endpoint. The inverse is defined modulo the unobserved tail, for which a representative must be chosen.
Once transformed, an appropriate FDA model can be fitted in linear coordinates and valid fitted densities recovered with the inverse transformation. @KokoszkaEtAl2019 use LQD and compositional transformations for **forecasting density-valued time series**, rather than for i.i.d. distribution-on-distribution regression. Earlier papers by @Delicado2011 and @KneipUtikal2001 concern dimension reduction and FPCA for samples of densities; they are useful background, but are not regression methods.
The main advantage is access to linear FDA tools while respecting the density constraints after back-transformation. The tradeoff is geometric: the $L^2$ norm of transformed functions is generally not the Wasserstein distance, so the transformation determines which discrepancies the fit treats as small.
### Intrinsic / Wasserstein-Based Approaches
Rather than transforming distributions away from their natural geometry, intrinsic approaches build the regression directly on the **Wasserstein space** $\mathcal{W}_2 = (\mathcal{P}_2(\mathbb{R}), W_2)$. This preserves the metric structure of distributions, ensuring that distances between fitted and observed responses are measured in the same units (earth mover's distance) that the application domain cares about.
Below are three intrinsic formulations that recently emerge:
- **Fréchet regression** [@PetersenMueller2019] treats a metric-space-valued response $Y$ and a Euclidean predictor $X\in\mathbb R^p$. It targets $m(x) = \arg\min_{\omega \in \Omega} \mathbb{E}\{d^2(\omega, Y) \mid X = x\}$ and develops both global least-squares-type weights and local-linear kernel weights. With a distribution-valued response and $d=W_2$, this becomes Wasserstein-Fréchet response regression; it does not make the predictor distribution-valued.
- **Tangent-space Wasserstein regression** [@ChenLinMueller2023] is developed for univariate distributions. Predictor and response distributions are separately linearized at their Wasserstein Fréchet means using $\log_\mu(\nu)=F_\nu^{-1}\circ F_\mu-\mathrm{id}$, and a functional linear operator is estimated between the two tangent spaces. A fitted response is mapped back with $\exp_\mu(g)=(g+\mathrm{id})\#\mu$, with a projection/boundary correction when the fitted tangent function is not admissible. The paper also treats distribution-to-scalar regression. It remarks that a multivariate extension is possible in principle, but the univariate formulas, computation, parallel transport, and asymptotic analysis do not carry over directly.
- **Optimal-transport-map regression** [@GhodratiPanaretos2022], which is the focus of the remainder of this lecture, models the regression effect **directly as a monotone transport map** $T_0$ pushing predictor distributions to response distributions. Rather than passing through a tangent space, the map $T_0 : \Omega \to \Omega$ is estimated by weighted isotonic regression. This approach offers sharp interpretability — each quantile of the predictor maps to the corresponding quantile of the fitted response — at the cost of being restricted to one-dimensional distributions.
Related response-regression work addresses different data structures. @PetersenLiuDivani2021 develop tests and confidence bands for Wasserstein-Fréchet regression with vector predictors and density responses. @ZhouMueller2024 regress a univariate distribution response on Euclidean covariates using empirical measures, borrowing information across responses when some are represented by very few observations. Neither paper is an extension of distribution-on-distribution tangent regression.
## Transport-Map Regression
### Motivation
Consider a setting where both the predictor and the response are one-dimensional probability distributions supported on a compact interval $\Omega = [a,b] \subset \mathbb{R}$. Examples include:
- **Predictor:** income distribution of a country at year $t$; **Response:** income distribution at year $t+1$.
- **Predictor:** distribution of asset returns during calm markets; **Response:** distribution of returns during stress periods.
In each case, the natural question is: *how does the distribution deform from predictor to response?* The transport-map regression framework of @GhodratiPanaretos2022 answers this by modeling the deformation directly.
### The Regression Map as a Fréchet Mean
Let $(X_i, Y_i)$, $i = 1,\ldots,n$, be independent pairs of probability distributions in $\mathcal{P}_2(\Omega)$. First define the conditional Fréchet mean **operator**
$$
f(x)=\arg\min_{b\in\mathcal P_2(\Omega)}
\mathbb E\!\left\{W_2^2(b,Y)\mid X=x\right\},
\qquad x\in\mathcal P_2(\Omega).
$$
The model of @GhodratiPanaretos2022 imposes the shape constraint that one common nondecreasing map $T_0:\Omega\to\Omega$ generates this operator:
$$
f(x) = T_0\#x,
$$
where $T_0\#x$ is the pushforward of the predictor distribution $x$. Thus $f$ maps distributions to distributions, whereas $T_0$ maps locations in their common support to locations. For a predictor quantile $Q_x(u)$, the corresponding conditional-mean response quantile is $T_0\{Q_x(u)\}$.
### The Observation Model
At the observation level, the response distribution $Y_i$ is modeled as a **random monotone deformation** of the conditional mean $T_0\#X_i$:
$$
Y_i = \epsilon_i\#(T_0\#X_i),
$$
where $\epsilon_i : \Omega \to \Omega$ are i.i.d. random monotone maps satisfying the centering condition $\mathbb{E}\{\epsilon_i(t)\} = t$ for every $t \in \Omega$. This condition ensures that the noise deforms but does not systematically shift the distribution — the Fréchet mean of $Y_i$ given $X_i$ is indeed $T_0\#X_i$.
The key structural assumptions are:
1. **Monotonicity of $T_0$:** $T_0$ is nondecreasing on $\Omega$. Any measurable map defines a valid pushforward, but monotonicity makes $T_0$ the one-dimensional optimal transport from $X_i$ to $T_0\#X_i$.
2. **Random deformation noise:** The $\epsilon_i$ are i.i.d. monotone optimal maps, independent of the predictors, and satisfy the identity-mean condition almost everywhere.
3. **Predictor regularity and compact support:** Predictor measures are absolutely continuous on the compact interval $\Omega=[a,b]$. The regression and error maps belong to the paper's class of almost-everywhere differentiable nondecreasing maps.
### Fréchet Least-Squares Estimation
The population criterion is the expected squared Wasserstein distance between the fitted and observed response distributions:
$$
M(T) = \frac{1}{2}\,\mathbb{E}\bigl\{W_2^2(T\#X, Y)\bigr\}.
$$
The true regression map is the minimizer: $T_0 = \arg\min_{T \in \mathcal{T}} M(T)$, where $\mathcal{T}$ is the class of monotone maps $\Omega \to \Omega$.
The empirical counterpart replaces the expectation with a sample average. Let $Q_n$ be the empirical distribution of the predictor distributions, $H_n = n^{-1}\sum_{i=1}^n \delta_{X_i}$. The estimator is:
$$
\hat{T}_n = \arg\min_{T \in \mathcal{T}_{H_n}} \frac{1}{2n}\sum_{i=1}^n W_2^2(T\#X_i, Y_i),
$$
where $\mathcal{T}_{H_n}$ denotes the admissible monotone-map class, with maps identified up to $H_n$-null sets.
### Reduction to Weighted Least Squares via Optimal Transport
In one dimension, the Wasserstein distance has a closed form in terms of quantile functions. For any two distributions $\mu, \nu$ on $\Omega$:
$$
W_2^2(\mu, \nu) = \int_0^1 |F_\mu^{-1}(u) - F_\nu^{-1}(u)|^2 \, du.
$$
Moreover, the optimal transport map from $X_i$ to $Y_i$ is
$$
T_{X_i \to Y_i} = F_{Y_i}^{-1} \circ F_{X_i}.
$$
When $X_i$ is absolutely continuous, this is the unique $X_i$-almost-everywhere nondecreasing optimal map pushing $X_i$ forward to $Y_i$.
The crucial insight is that $W_2^2(T\#X_i, Y_i)$ can be expressed as an $L^2$ distance between $T$ and the pair-specific optimal transport map $T_{X_i \to Y_i}$. Specifically, on the support of $X_i$:
$$
W_2^2(T\#X_i, Y_i) = \int_\Omega |T(s) - T_{X_i \to Y_i}(s)|^2 \, dX_i(s).
$$
Indeed, monotonicity gives $F_{T\#X_i}^{-1}=T\circ F_{X_i}^{-1}$ almost everywhere, so the quantile formula and the substitution $s=F_{X_i}^{-1}(u)$ give the displayed identity directly. This argument does not require $T$ to have an inverse.
Substituting into the empirical criterion:
$$
M_n(T) = \frac{1}{2n}\sum_{i=1}^n \int_\Omega |T(s) - T_{X_i \to Y_i}(s)|^2 \, dX_i(s).
$$
Define the **aggregate predictor measure** $Q_n=n^{-1}\sum_iX_i$ and the finite signed measure
$$
A_n(B)=\frac1n\sum_{i=1}^n\int_B T_{X_i\to Y_i}(s)\,dX_i(s).
$$
Since $A_n\ll Q_n$, let $\bar T_n=dA_n/dQ_n$. Expanding the square yields the equivalent minimization
$$
\hat{T}_n = \arg\min_{T \text{ monotone}} \int_\Omega |T(s) - \bar{T}_n(s)|^2 \, dQ_n(s),
$$
This is a **weighted isotonic regression** problem. If the $X_i$ have densities $f_i$, then $\bar T_n(s)=\sum_i f_i(s)T_{X_i\to Y_i}(s)/\sum_i f_i(s)$ wherever the denominator is positive.
### Estimation via the Pool-Adjacent-Violators Algorithm (PAVA)
The discretized problem is solved by the **pool-adjacent-violators algorithm** (PAVA), whose roots include @AyerEtAl1955 and @BarlowBrunk1972. Choose grid points $s_1<\cdots<s_K$ and cells $I_k$ around them. For pair $i$, set
$$
y_{ik}=T_{X_i\to Y_i}(s_k),\qquad w_{ik}=X_i(I_k).
$$
Equivalently, use the aggregate weight $w_k=\sum_iw_{ik}$ and target $v_k=\sum_iw_{ik}y_{ik}/w_k$. PAVA then proceeds as follows:
1. **Initialization:** Set $\hat{T}(s_k) = v_k$ and block weight $W_k = w_k$ for each $k = 1,\ldots,K$.
2. **Forward scan:** For $k = 1,\ldots,K-1$, if $\hat{T}(s_k) > \hat{T}(s_{k+1})$, the monotonicity constraint is violated. **Pool** the two adjacent blocks: replace both values by their weighted average and merge their weights.
3. **Repeat** until the estimated values are nondecreasing.
The fitted **grid values** are constant within pooled blocks. A function on all of $\Omega$ still requires an interpolation convention; @GhodratiPanaretos2022 use monotone interpolation of the fitted grid values, so the final displayed map need not be a step function.
PAVA is computationally efficient ($O(K)$ after sorting) and gives the unique fitted vector at grid points having positive aggregate weight. In the continuous formulation, uniqueness is in the $L^2(Q_n)$ sense, not necessarily pointwise outside the support of $Q_n$.
::: {.proof title="Strict Convexity of the Criterion"}
The criterion $M_n(T)$ is of the form $\sum_{i=1}^n \int (T - g_i)^2 d\mu_i$ with $g_i = T_{X_i \to Y_i}$ and $\mu_i = X_i/n$. For any two maps $T_1 \neq T_2$ (differing on a set of positive $Q_n$-measure) and $\lambda \in (0,1)$:
$$
\begin{aligned}
M_n(\lambda T_1 + (1-\lambda)T_2)
&= \sum_i \int \bigl(\lambda T_1 + (1-\lambda)T_2 - g_i\bigr)^2 d\mu_i \\
&= \sum_i \int \bigl(\lambda(T_1 - g_i) + (1-\lambda)(T_2 - g_i)\bigr)^2 d\mu_i \\
&< \lambda \sum_i \int (T_1 - g_i)^2 d\mu_i + (1-\lambda) \sum_i \int (T_2 - g_i)^2 d\mu_i \\
&= \lambda M_n(T_1) + (1-\lambda) M_n(T_2),
\end{aligned}
$$
where the strict inequality follows from the strict convexity of $x \mapsto x^2$ and the fact that $T_1 \neq T_2$ on a set of positive aggregate measure. Hence $M_n$ admits a minimizer that is unique as an element of $L^2(Q_n)$; values on $Q_n$-null sets are not identified. Otherwise, if $T_1$ and $T_2$ are minimizers (so that $M_n(T_1) = M_n(T_2)$) but $T_1\neq T_2$ on a set of positive $Q_n$-measure, then the strict convexity argument above would yield a strictly smaller value at $\lambda T_1 + (1-\lambda)T_2$, contradicting minimality.
:::
### Convergence Rates
The convergence results of @GhodratiPanaretos2022 rest on the structural assumptions stated [earlier](#the-observation-model) — compact support, absolute continuity of predictors, monotonicity of $T_0$, and the identity-mean noise condition — together with two additional regularity conditions on the map class:
::: {.callout-important title="Assumptions for Convergence Rates"}
**Assumption 3.1 (Absolute continuity and compact support).** Any predictor distribution $\mu$ in the support of the law of $X$ is absolutely continuous with respect to Lebesgue measure and is supported on the compact interval $\Omega = [a,b]$.
**Assumption 3.2 (Bounded derivative / smoothness).** The true regression map $T_0$ and the random noise maps $\epsilon_i$ belong to the class
$$
\mathcal{T} := \{T : \Omega\rightarrow\Omega: 0\leq T^\prime(x) \leq \infty \text{ for }Q-\text{a.e. }x \in \Omega\},
$$
where $Q(A)=\mathbb{E}\{X(A)\}$ defines a measure on $\Omega$.
:::
Under these assumptions, @GhodratiPanaretos2022 establish:
- **Fully observed distributions** (distributions observed as continuous objects):
$$
\|\hat{T}_n - T_0\|_{L^2(Q)} = O_{\mathbb{P}}(n^{-1/3}).
$$
- **Partially observed distributions** (each $X_i$, $Y_i$ observed through $m$ i.i.d. samples):
$$
\|\hat{T}_n - T_0\|_{L^2(Q)} = O_{\mathbb{P}}(n^{-1/3} + r_m^{-1/2}),
$$
where $r_m^{-1}$ is the $W_2$ convergence rate of the smoothed estimator of each latent predictor and response distribution. For example, if distribution estimation has $W_2$ error $O_{\mathbb P}(m^{-1/2})$, then $r_m=m^{1/2}$ and the added term is $m^{-1/4}$, not $m^{-1/2}$.
The cube-root rate comes from the empirical-process complexity of the bounded monotone-map class $\mathcal{T}$, which is a Donsker class whose entropy integral yields the $n^{-1/3}$ scaling — the same rate as in classical isotonic regression and monotone density estimation.
### Residual Diagnostics
A natural definition of residuals in transport-map regression comes from the observation model $Y_i = \epsilon_i\#(T_0\#X_i)$. The **residual transport map** from the fitted response $\hat{Y}_i = \hat{T}_n\#X_i$ to the observed response $Y_i$ is:
$$
\hat{R}_i = T_{\hat{Y}_i \to Y_i} = F_{Y_i}^{-1} \circ F_{\hat{Y}_i}.
$$
Under correct model specification, these maps should fluctuate around the identity map $\mathrm{id}(t)=t$, although their empirical average is not constrained to equal the identity exactly. Systematic pointwise deviations can reveal lack of fit or clusters.
## Interactive Demo: Transport-Map Regression
The following demo simulates distribution-on-distribution regression in one dimension. Given $n$ pairs of distributions $(X_i, Y_i)$, we estimate the regression map $\hat{T}_n$ via weighted isotonic regression and display the fitted responses.
::: {.callout-note title="Simulation Design"}
- **Predictor distributions** $X_i$: smooth, full-support laws on $\Omega=[-5,15]$, generated through random strictly increasing quantile functions.
- **True regression map** $T_0(s)=s+2\sin\{2\pi(s+5)/20\}$, which is smooth, strictly increasing, and maps both endpoints to themselves.
- **Response distributions** $Y_i$: generated exactly as $Y_i=\epsilon_i\#(T_0\#X_i)$, with monotone endpoint-preserving error maps satisfying $\mathbb E\{\epsilon_i(t)\}=t$.
- **Estimation:** weighted PAVA using cell probabilities $X_i(I_k)$ on a grid of 101 points.
:::
```{ojs}
//| label: transportmap-regression-demo
n14_control = Inputs.range([10, 100], {step: 5, value: 40, label: "Sample size n"})
noise14_control = Inputs.range([0.05, 1.0], {step: 0.05, value: 0.3, label: "Deformation noise σ"})
n14 = Generators.input(n14_control)
noise14 = Generators.input(noise14_control)
tmr_controls_view = html`
<style>
.tmr-slider-grid { display:grid; grid-template-columns:repeat(2,minmax(0,1fr)); gap:6px 20px; max-width:620px; margin:0 auto 12px; font:0.85em system-ui,sans-serif; }
.tmr-slider-grid > * { min-width:0; margin:0; }
.tmr-slider-grid input[type="number"] { width:7.5rem !important; }
@media (max-width:500px) { .tmr-slider-grid { grid-template-columns:1fr; } }
</style>
<div class="tmr-slider-grid">
<div>${n14_control}</div>
<div>${noise14_control}</div>
</div>`
function runTransportMapReg(n, sigma, seed) {
const rng = (function(a) {
return function() { a |= 0; a = a + 0x6D2B79F5 | 0; var t = Math.imul(a ^ a >>> 15, 1 | a); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; };
})(seed || 137);
const a = -5, b = 15, L = b - a;
const K = 100;
const grid = Array.from({length: K + 1}, (_, k) => a + L * k / K);
// Endpoint-preserving monotone regression map:
// T0'(s) = 1 + 0.2π cos(2π(s-a)/L) > 0.
function T0(s) {
return s + 2 * Math.sin(2 * Math.PI * (s - a) / L);
}
// Each predictor is defined by a strictly increasing quantile function
// Q_i(u). The coefficient bounds ensure Q_i'(u) > 0 and Q_i(0)=a,
// Q_i(1)=b, so every X_i has full support on Ω.
function makePredictor(c1, c2) {
const quantile = u => a + L * (
u
+ c1 * Math.sin(2 * Math.PI * u) / (2 * Math.PI)
+ c2 * Math.sin(4 * Math.PI * u) / (4 * Math.PI)
);
const cdf = s => {
if (s <= a) return 0;
if (s >= b) return 1;
let lo = 0, hi = 1;
for (let iter = 0; iter < 45; iter++) {
const mid = (lo + hi) / 2;
if (quantile(mid) < s) lo = mid; else hi = mid;
}
return (lo + hi) / 2;
};
return {quantile, cdf};
}
// ε_i(t)=t+c_i L/(2π) sin(2π(t-a)/L) is endpoint-preserving and
// monotone because |c_i|<1. Symmetry of c_i gives E{ε_i(t)}=t.
function errorMap(t, coefficient) {
return t + coefficient * L * Math.sin(2 * Math.PI * (t - a) / L) / (2 * Math.PI);
}
const predictors = [];
const errorCoefficients = [];
const pairMaps = [];
for (let i = 0; i < n; i++) {
const c1 = 0.9 * (rng() - 0.5);
const c2 = 0.4 * (rng() - 0.5);
const err = 1.5 * sigma * (rng() - 0.5);
const predictor = makePredictor(c1, c2);
predictors.push(predictor);
errorCoefficients.push(err);
pairMaps.push(grid.map(s => errorMap(T0(s), err)));
}
// Discretize X_i with cell masses, as in the paper's Section 3.5.
const weights = new Array(K + 1).fill(0);
const targets = new Array(K + 1).fill(0);
for (let k = 0; k <= K; k++) {
const left = k === 0 ? a : (grid[k - 1] + grid[k]) / 2;
const right = k === K ? b : (grid[k] + grid[k + 1]) / 2;
for (let i = 0; i < n; i++) {
const wik = predictors[i].cdf(right) - predictors[i].cdf(left);
weights[k] += wik;
targets[k] += wik * pairMaps[i][k];
}
targets[k] = weights[k] > 0 ? targets[k] / weights[k] : grid[k];
}
// Correct weighted PAVA using a stack of pooled blocks.
function pava(vals, wts) {
const blocks = [];
for (let k = 0; k < vals.length; k++) {
blocks.push({start: k, end: k, weight: wts[k], sum: wts[k] * vals[k]});
while (blocks.length > 1) {
const right = blocks[blocks.length - 1];
const left = blocks[blocks.length - 2];
if (left.sum / left.weight <= right.sum / right.weight + 1e-12) break;
blocks.splice(blocks.length - 2, 2, {
start: left.start,
end: right.end,
weight: left.weight + right.weight,
sum: left.sum + right.sum
});
}
}
const result = new Array(vals.length).fill(0);
for (const block of blocks) {
const mean = block.sum / block.weight;
for (let k = block.start; k <= block.end; k++) result[k] = mean;
}
return result;
}
const T0vals = grid.map(s => T0(s));
const TnHat = pava(targets, weights);
function interp(values, s) {
if (s <= a) return values[0];
if (s >= b) return values[K];
const pos = (s - a) / L * K;
const k = Math.floor(pos);
const frac = pos - k;
return values[k] * (1 - frac) + values[k + 1] * frac;
}
// Descriptive Wasserstein R² using the sample response barycenter.
const U = 200;
const xQuants = [], yQuants = [], fitQuants = [];
for (let i = 0; i < n; i++) {
const xq = [], yq = [], fq = [];
for (let j = 0; j < U; j++) {
const u = (j + 0.5) / U;
const x = predictors[i].quantile(u);
xq.push(x);
yq.push(errorMap(T0(x), errorCoefficients[i]));
fq.push(interp(TnHat, x));
}
xQuants.push(xq);
yQuants.push(yq);
fitQuants.push(fq);
}
const baryQuant = Array.from({length: U}, (_, j) =>
yQuants.reduce((sum, q) => sum + q[j], 0) / n
);
let totalW2Sq = 0, resW2Sq = 0;
for (let i = 0; i < n; i++) {
for (let j = 0; j < U; j++) {
totalW2Sq += (yQuants[i][j] - baryQuant[j]) ** 2 / (n * U);
resW2Sq += (yQuants[i][j] - fitQuants[i][j]) ** 2 / (n * U);
}
}
const w2R2 = totalW2Sq > 1e-10 ? 1 - resW2Sq / totalW2Sq : 0;
// Empirical L²(Q_n) squared error, matching the paper's risk.
let ise = 0;
for (let k = 0; k <= K; k++) {
ise += (weights[k] / n) * (T0vals[k] - TnHat[k]) ** 2;
}
return {grid, T0vals, TnHat, targets, weights, n, sigma, w2R2, ise, pairMaps};
}
tmRegRes = runTransportMapReg(n14, noise14, 137);
html`
<div style="font-family: system-ui, sans-serif; max-width: 900px;">
<h4>Distribution-on-Distribution Regression via Transport Maps</h4>
<div style="display: flex; gap: 16px; flex-wrap: wrap;">
<div>
<svg width="420" height="300" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
const mg = {top: 20, right: 100, bottom: 35, left: 50};
const pW = 420 - mg.left - mg.right, pH = 300 - mg.top - mg.bottom;
const {grid, T0vals, TnHat} = tmRegRes;
const a = grid[0], b = grid[grid.length - 1];
const allV = [...T0vals, ...TnHat, ...grid];
const yMin = Math.min(...allV) - 0.5, yMax = Math.max(...allV) + 0.5;
const xS = s => mg.left + (s - a)/(b - a) * pW;
const yS = v => mg.top + pH - (v - yMin)/(yMax - yMin) * pH;
// Identity line
const idLine = grid.map((s,i) => `${i===0?'M':'L'} ${xS(s)} ${yS(s)}`).join(' ');
// True T0
const trueLine = grid.map((s,i) => `${i===0?'M':'L'} ${xS(s)} ${yS(T0vals[i])}`).join(' ');
// Estimated T_hat
const estLine = grid.map((s,i) => `${i===0?'M':'L'} ${xS(s)} ${yS(TnHat[i])}`).join(' ');
return `<line x1="${mg.left}" y1="${mg.top}" x2="${mg.left}" y2="${mg.top+pH}" stroke="#adb5bd"/>
<line x1="${mg.left}" y1="${mg.top+pH}" x2="${mg.left+pW}" y2="${mg.top+pH}" stroke="#adb5bd"/>
<text x="${mg.left+pW/2}" y="${mg.top+pH+25}" text-anchor="middle" font-size="11">s ∈ Ω (support)</text>
<text x="${mg.left-42}" y="${mg.top+pH/2}" text-anchor="middle" font-size="11" transform="rotate(-90,${mg.left-42},${mg.top+pH/2})">T(s)</text>
<path d="${idLine}" fill="none" stroke="#dee2e6" stroke-width="1" stroke-dasharray="4,4"/>
<path d="${trueLine}" fill="none" stroke="#2b8a3e" stroke-width="2.5" stroke-dasharray="6,3"/>
<path d="${estLine}" fill="none" stroke="#1971c2" stroke-width="2.5"/>
<rect x="${mg.left}" y="${mg.top}" width="${pW}" height="${pH}" fill="none" stroke="#dee2e6"/>
<g transform="translate(${mg.left+5},${mg.top+5})">
<line x1="0" y1="3" x2="15" y2="3" stroke="#dee2e6" stroke-width="1" stroke-dasharray="4,4"/><text x="18" y="7" font-size="10">Identity</text>
<line x1="0" y1="16" x2="15" y2="16" stroke="#2b8a3e" stroke-width="2.5" stroke-dasharray="6,3"/><text x="18" y="20" font-size="10">True T₀</text>
<line x1="0" y1="29" x2="15" y2="29" stroke="#1971c2" stroke-width="2.5"/><text x="18" y="33" font-size="10">Estimated T̂ₙ</text>
</g>`;
})()}
</svg>
</div>
<div>
<svg width="380" height="300" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
const mg = {top: 20, right: 20, bottom: 35, left: 45};
const pW = 380 - mg.left - mg.right, pH = 300 - mg.top - mg.bottom;
const {pairMaps, grid} = tmRegRes;
const a = grid[0], b = grid[grid.length - 1];
const xS = s => mg.left + (s - a)/(b - a) * pW;
const yS = v => mg.top + pH - (v - a)/(b - a) * pH;
// Show a few pair-specific OT maps as thin lines
const step = Math.max(1, Math.floor(pairMaps.length / 8));
const mapLines = [];
for (let i = 0; i < pairMaps.length; i += step) {
const alpha = 0.15 + 0.15 * (i / pairMaps.length);
const pathParts = [];
for (let k = 0; k < grid.length; k++) {
const cmd = k === 0 ? 'M' : 'L';
pathParts.push(cmd + ' ' + xS(grid[k]) + ' ' + yS(pairMaps[i][k]));
}
mapLines.push('<path d="' + pathParts.join(' ') + '" fill="none" stroke="#868e96" stroke-width="0.8" opacity="' + alpha + '"/>');
}
const idPathParts = [];
for (let i = 0; i < grid.length; i++) {
const cmd = i === 0 ? 'M' : 'L';
idPathParts.push(cmd + ' ' + xS(grid[i]) + ' ' + yS(grid[i]));
}
const idPathD = idPathParts.join(' ');
return `<line x1="${mg.left}" y1="${mg.top}" x2="${mg.left}" y2="${mg.top+pH}" stroke="#adb5bd"/>
<line x1="${mg.left}" y1="${mg.top+pH}" x2="${mg.left+pW}" y2="${mg.top+pH}" stroke="#adb5bd"/>
<text x="${mg.left+pW/2}" y="${mg.top+pH+25}" text-anchor="middle" font-size="11">s ∈ Ω</text>
<text x="${mg.left-38}" y="${mg.top+pH/2}" text-anchor="middle" font-size="11" transform="rotate(-90,${mg.left-38},${mg.top+pH/2})">T_{Xᵢ→Yᵢ}(s)</text>
<path d="${idPathD}" fill="none" stroke="#dee2e6" stroke-width="1" stroke-dasharray="4,4"/>
${mapLines.join('')}
<rect x="${mg.left}" y="${mg.top}" width="${pW}" height="${pH}" fill="none" stroke="#dee2e6"/>`;
})()}
</svg>
</div>
</div>
<div style="margin-top:12px;padding:12px;background:#f8f9fa;border-radius:6px;font-size:0.92em;">
<b>Descriptive Wasserstein R² = ${tmRegRes.w2R2.toFixed(3)}</b> |
<b>‖T̂ₙ−T₀‖²<sub>L²(Qₙ)</sub> = ${tmRegRes.ise.toFixed(4)}</b> |
n = ${tmRegRes.n} | σ<sub>noise</sub> = ${tmRegRes.sigma}
</div>
<p style="margin-top:8px;font-size:0.85em;color:#495057;">
<b>Left:</b> The true regression map T₀ (green dashed) and the PAVA estimate T̂ₙ (blue solid).
The gray dotted line is the identity — deviations capture the regression effect.<br>
<b>Right:</b> A sample of pair-specific optimal transport maps T<sub>Xᵢ→Yᵢ</sub> (thin gray curves).
These serve as pseudo-observations of the regression map; their vertical scatter at each s reflects the deformation noise.
</p>
</div>`
```
::: {.callout-tip title="Try these experiments"}
- **Increase $n$:** The estimate tends to approach $T_0$, although a realized error curve need not decrease at every slider step. A cube-root $L^2$ rate means that asymptotically an eightfold sample-size increase halves the order of the error.
- **Increase noise $\sigma$:** The pair-specific maps scatter more widely around $T_0$; the Wasserstein $R^2$ drops.
- **Inspect the map:** Where $T_0(s) > s$, the regression effect shifts mass to the right (distributions become "larger" at those quantiles). Where $T_0(s) < s$, mass shifts left.
:::
## Application: Income-Distribution Change — From One Year to the Next
### From Distribution Comparison to Distributional Regression
In [Lecture 11](lecture-11.qmd#sec-income-application), we used the Wasserstein distance to compare income distributions across countries and over time. That analysis treated each country-year as a point in Wasserstein space and asked: *how far apart are two income distributions?*
Distribution-on-distribution regression allows us to ask a dynamic aggregate question: **given a country's income distribution this year, what income distribution should we expect next year?** This concerns distributional change. It is not, by itself, a model of individual income mobility, because repeated cross-sectional distributions do not identify which people changed rank.
### Interpreting the Transport Map
The transport-map framework gives a direct quantile interpretation, with two important qualifications:
1. **$s$ is an income level, not a percentile.** For percentile $u$, the base-year income is $Q_{X_c}(u)$ and the fitted follow-up income is $\hat T_n\{Q_{X_c}(u)\}$. Thus an absolute quantile-change curve is $u\mapsto\hat T_n\{Q_{X_c}(u)\}-Q_{X_c}(u)$. A conventional growth-incidence curve would usually scale this change by $Q_{X_c}(u)$.
2. **Monotonicity preserves quantile order, not individual ranks.** If $u_1<u_2$, the fitted quantiles remain ordered. This is a statement about the two marginal distributions and their optimal quantile coupling, not evidence that the same individuals retain their ranks.
3. **Residuals identify anomalous distributional change.** The maps $\hat R_c=T_{\hat Y_c\to Y_c}$ show where the observed follow-up quantiles differ from the fitted distribution. They do not identify person-level transitions.
### Data Structure
Let $\{(X_c,Y_c)\}_{c=1}^C$ denote paired income distributions for $C$ countries, where $X_c$ is the base-year distribution and $Y_c$ is the follow-up distribution. The fitted common map answers: *if a base-year quantile has income level $s$, to which follow-up income level does the model send it?* This is a hypothetical application. @GhodratiPanaretos2022 instead apply the method to 1983 and 2013 age-at-death distributions for 37 countries and use residual maps to compare observed and fitted mortality change.
### Interactive Demo: Income-Distribution Change Across Countries
The demo below applies transport-map regression to a simulated panel of countries. Each country has a base-year income distribution $X_c$ (a perturbed lognormal mixture) and a follow-up distribution $Y_c$ generated by a common monotone map $T_0$ plus country-specific deformation noise. Adjust the inequality/growth shape $\gamma$, the number of countries, and the noise level; the fitted common map $\hat T_n$ is computed by weighted PAVA.
::: {.callout-note title="Simulation Design"}
- **Base-year distributions** $X_c$: random perturbations of three lognormal-mixture profiles (developed / emerging / developing), all supported on $[0.1,200]$ (in thousands of dollars).
- **True regression map:** $T_0(s) = s_{\min} + (s_{\max}-s_{\min})\bigl[x + \gamma x(1-x)(\tfrac12-x)\bigr]$, $x=(s-s_{\min})/(s_{\max}-s_{\min})$. Positive $\gamma$ raises middle incomes relative to the tails; negative $\gamma$ compresses the middle.
- **Follow-up distributions:** $Y_c = \epsilon_c\#(T_0\#X_c)$ with endpoint-preserving monotone noise satisfying $\mathbb E\{\epsilon_c(t)\}=t$.
- **Estimation:** weighted PAVA on a log-spaced grid using cell probabilities $X_c(I_k)$.
:::
```{ojs}
//| label: income-change-demo
incTmr_countries_control = Inputs.range([5, 50], {step: 5, value: 20, label: "Number of countries C"})
incTmr_noise_control = Inputs.range([0.0, 0.5], {step: 0.05, value: 0.15, label: "Deformation noise σ"})
incTmr_gamma_control = Inputs.range([-1.5, 1.5], {step: 0.1, value: 0.6, label: "Inequality/growth shape γ"})
incTmr_index_control = Inputs.range([0, 49], {step: 1, value: 0, label: "Highlighted country index"})
incTmr_countries = Generators.input(incTmr_countries_control)
incTmr_noise = Generators.input(incTmr_noise_control)
incTmr_gamma = Generators.input(incTmr_gamma_control)
incTmr_index = Generators.input(incTmr_index_control)
incTmr_controls_view = html`
<style>
.inctmr-slider-grid { display:grid; grid-template-columns:repeat(2,minmax(0,1fr)); gap:6px 20px; max-width:620px; margin:0 auto 12px; font:0.85em system-ui,sans-serif; }
.inctmr-slider-grid > * { min-width:0; margin:0; }
.inctmr-slider-grid input[type="number"] { width:7.5rem !important; }
@media (max-width:500px) { .inctmr-slider-grid { grid-template-columns:1fr; } }
</style>
<div class="inctmr-slider-grid">
<div>${incTmr_countries_control}</div>
<div>${incTmr_noise_control}</div>
<div>${incTmr_gamma_control}</div>
<div>${incTmr_index_control}</div>
</div>`
function incTmr_erf(x) {
const sign = x >= 0 ? 1 : -1;
x = Math.abs(x);
const a1 = 0.254829592, a2 = -0.284496736, a3 = 1.421413741;
const a4 = -1.453152027, a5 = 1.061405429, p = 0.3275911;
const t = 1 / (1 + p * x);
const y = 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-x * x);
return sign * y;
}
function incTmr_phi(z) { return 0.5 * (1 + incTmr_erf(z / Math.sqrt(2))); }
function incTmr_makeIncomeProfile(base, rng) {
const w1 = Math.max(0.5, Math.min(0.95, base.w1 + 0.06 * (rng() - 0.5)));
return {
label: base.label,
w1: w1, w2: 1 - w1,
mu1: base.mu1 + 0.12 * (rng() - 0.5),
sigma1: Math.max(0.3, base.sigma1 + 0.08 * (rng() - 0.5)),
mu2: base.mu2 + 0.12 * (rng() - 0.5),
sigma2: Math.max(0.4, base.sigma2 + 0.08 * (rng() - 0.5))
};
}
function incTmr_makeCDF(profile) {
const {w1, mu1, sigma1, w2, mu2, sigma2} = profile;
return function(x) {
if (x <= 0) return 0;
const z1 = (Math.log(x) - mu1) / sigma1;
let f = w1 * incTmr_phi(z1);
if (w2 > 0) {
const z2 = (Math.log(x) - mu2) / sigma2;
f += w2 * incTmr_phi(z2);
}
return f;
};
}
function incTmr_makeQuantile(profile) {
const cdf = incTmr_makeCDF(profile);
return function(u) {
const uu = Math.max(1e-10, Math.min(1 - 1e-10, u));
let lo = 0.1, hi = 500;
while (cdf(hi) < uu && hi < 10000) hi *= 2;
for (let k = 0; k < 60; k++) {
const mid = (lo + hi) / 2;
if (cdf(mid) < uu) lo = mid; else hi = mid;
}
return (lo + hi) / 2;
};
}
function incTmr_makeDensity(profile) {
const {w1, mu1, sigma1, w2, mu2, sigma2} = profile;
return function(x) {
if (x <= 0) return 0;
let d = 0;
const z1 = (Math.log(x) - mu1) / sigma1;
d += (w1 / (x * sigma1 * Math.sqrt(2 * Math.PI))) * Math.exp(-0.5 * z1 * z1);
if (w2 > 0) {
const z2 = (Math.log(x) - mu2) / sigma2;
d += (w2 / (x * sigma2 * Math.sqrt(2 * Math.PI))) * Math.exp(-0.5 * z2 * z2);
}
return d;
};
}
function incTmr_densityFromQuantile(Q, xGrid) {
const n = xGrid.length;
const cdf = xGrid.map(x => {
if (x <= Q(1e-10)) return 0;
if (x >= Q(1 - 1e-10)) return 1;
let lo = 0, hi = 1;
for (let k = 0; k < 55; k++) {
const mid = (lo + hi) / 2;
if (Q(mid) < x) lo = mid; else hi = mid;
}
return (lo + hi) / 2;
});
return cdf.map((u, i) => {
if (i === 0) return Math.max(0, cdf[1] / (xGrid[1] - xGrid[0]));
if (i === n - 1) return Math.max(0, (1 - cdf[n - 2]) / (xGrid[n - 1] - xGrid[n - 2]));
return Math.max(0, (cdf[i + 1] - cdf[i - 1]) / (xGrid[i + 1] - xGrid[i - 1]));
});
}
function incTmr_pava(vals, wts) {
const blocks = [];
for (let k = 0; k < vals.length; k++) {
blocks.push({start: k, end: k, weight: wts[k], sum: wts[k] * vals[k]});
while (blocks.length > 1) {
const right = blocks[blocks.length - 1];
const left = blocks[blocks.length - 2];
if (left.sum / left.weight <= right.sum / right.weight + 1e-12) break;
blocks.splice(blocks.length - 2, 2, {
start: left.start, end: right.end,
weight: left.weight + right.weight,
sum: left.sum + right.sum
});
}
}
const result = new Array(vals.length).fill(0);
for (const block of blocks) {
const mean = block.sum / block.weight;
for (let k = block.start; k <= block.end; k++) result[k] = mean;
}
return result;
}
function incTmr_runIncomeChangeDemo(n, sigma, gamma, highlightIdx, seed) {
const rng = (function(a) {
return function() { a |= 0; a = a + 0x6D2B79F5 | 0; var t = Math.imul(a ^ a >>> 15, 1 | a); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; };
})(seed || 314);
const sMin = 0.1, sMax = 200, L = sMax - sMin;
const K = 80;
const grid = Array.from({length: K + 1}, (_, k) => sMin * Math.pow(sMax / sMin, k / K));
function T0(s) {
const x = (s - sMin) / L;
const h = x * (1 - x) * (0.5 - x);
return sMin + L * (x + gamma * h);
}
function errorMap(s, coeff) {
const x = (s - sMin) / L;
return s + coeff * L * Math.sin(2 * Math.PI * x) / (2 * Math.PI);
}
const baseProfiles = [
{label: "Developed", w1: 0.88, mu1: Math.log(38), sigma1: 0.45, w2: 0.12, mu2: Math.log(95), sigma2: 0.55},
{label: "Emerging", w1: 0.82, mu1: Math.log(18), sigma1: 0.60, w2: 0.18, mu2: Math.log(70), sigma2: 0.70},
{label: "Developing", w1: 0.75, mu1: Math.log(8), sigma1: 0.70, w2: 0.25, mu2: Math.log(55), sigma2: 0.85}
];
const countries = [];
for (let i = 0; i < n; i++) {
const base = baseProfiles[i % baseProfiles.length];
const profile = incTmr_makeIncomeProfile(base, rng);
countries.push({
profile: profile,
quantile: incTmr_makeQuantile(profile),
cdf: incTmr_makeCDF(profile),
density: incTmr_makeDensity(profile),
errorCoeff: 1.8 * sigma * (rng() - 0.5)
});
}
const weights = new Array(K + 1).fill(0);
const targets = new Array(K + 1).fill(0);
for (let k = 0; k <= K; k++) {
const left = k === 0 ? sMin : (grid[k - 1] + grid[k]) / 2;
const right = k === K ? sMax : (grid[k] + grid[k + 1]) / 2;
for (let i = 0; i < n; i++) {
const wik = countries[i].cdf(right) - countries[i].cdf(left);
weights[k] += wik;
targets[k] += wik * errorMap(T0(grid[k]), countries[i].errorCoeff);
}
targets[k] = weights[k] > 0 ? targets[k] / weights[k] : grid[k];
}
const T0vals = grid.map(s => T0(s));
const TnHat = incTmr_pava(targets, weights);
function interp(values, s) {
if (s <= sMin) return values[0];
if (s >= sMax) return values[K];
const pos = Math.log(s / sMin) / Math.log(sMax / sMin) * K;
const k = Math.floor(pos);
const frac = pos - k;
return values[k] * (1 - frac) + values[k + 1] * frac;
}
const U = 200;
const yQuants = [], fitQuants = [];
for (let i = 0; i < n; i++) {
const yq = [], fq = [];
for (let j = 0; j < U; j++) {
const u = (j + 0.5) / U;
const x = Math.min(Math.max(sMin, countries[i].quantile(u)), sMax);
yq.push(errorMap(T0(x), countries[i].errorCoeff));
fq.push(interp(TnHat, x));
}
yQuants.push(yq);
fitQuants.push(fq);
}
const baryQuant = Array.from({length: U}, (_, j) => yQuants.reduce((sum, q) => sum + q[j], 0) / n);
let totalW2Sq = 0, resW2Sq = 0;
for (let i = 0; i < n; i++) {
for (let j = 0; j < U; j++) {
totalW2Sq += (yQuants[i][j] - baryQuant[j]) ** 2 / (n * U);
resW2Sq += (yQuants[i][j] - fitQuants[i][j]) ** 2 / (n * U);
}
}
const w2R2 = totalW2Sq > 1e-10 ? 1 - resW2Sq / totalW2Sq : 0;
let ise = 0;
for (let k = 0; k <= K; k++) {
ise += (weights[k] / n) * (T0vals[k] - TnHat[k]) ** 2;
}
const idx = Math.min(Math.max(0, highlightIdx), n - 1);
const sel = countries[idx];
const xGrid = Array.from({length: 300}, (_, i) => sMin * Math.pow(sMax / sMin, i / 299));
const basePdf = xGrid.map(x => sel.density(x));
const truePdf = incTmr_densityFromQuantile(u => { const x = Math.min(Math.max(sMin, sel.quantile(u)), sMax); return T0(x); }, xGrid);
const obsPdf = incTmr_densityFromQuantile(u => { const x = Math.min(Math.max(sMin, sel.quantile(u)), sMax); return errorMap(T0(x), sel.errorCoeff); }, xGrid);
const fitPdf = incTmr_densityFromQuantile(u => { const x = Math.min(Math.max(sMin, sel.quantile(u)), sMax); return interp(TnHat, x); }, xGrid);
const nUchange = 100;
const uGrid = Array.from({length: nUchange}, (_, j) => (j + 0.5) / nUchange);
const changeTrue = uGrid.map(u => { const x = Math.min(Math.max(sMin, sel.quantile(u)), sMax); return T0(x) - x; });
const changeHat = uGrid.map(u => { const x = Math.min(Math.max(sMin, sel.quantile(u)), sMax); return interp(TnHat, x) - x; });
return {
grid, T0vals, TnHat, targets, weights, n, sigma, gamma, w2R2, ise,
xGrid, basePdf, truePdf, obsPdf, fitPdf,
uGrid, changeTrue, changeHat,
selectedLabel: sel.profile.label, selectedIdx: idx
};
}
incTmrResult = incTmr_runIncomeChangeDemo(incTmr_countries, incTmr_noise, incTmr_gamma, Math.min(incTmr_index, incTmr_countries - 1), 314);
html`
<div style="font-family: system-ui, sans-serif; max-width: 900px;">
<h4>Income-Distribution Change: Transport-Map Regression</h4>
${incTmr_controls_view}
<div style="display: flex; gap: 16px; flex-wrap: wrap;">
<div style="flex: 1; min-width: 420px;">
<svg width="100%" height="260" viewBox="0 0 440 260" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
const margin = {top: 18, right: 16, bottom: 32, left: 50};
const plotW = 440 - margin.left - margin.right;
const plotH = 260 - margin.top - margin.bottom;
const {grid, T0vals, TnHat} = incTmrResult;
const a = grid[0], b = grid[grid.length - 1];
const allV = [...T0vals, ...TnHat, ...grid];
const yMin = Math.min(...allV) - 2, yMax = Math.max(...allV) + 2;
const xS = s => margin.left + (s - a) / (b - a) * plotW;
const yS = v => margin.top + plotH - (v - yMin) / (yMax - yMin) * plotH;
const idLine = grid.map((s, i) => `${i === 0 ? 'M' : 'L'} ${xS(s)} ${yS(s)}`).join(' ');
const trueLine = grid.map((s, i) => `${i === 0 ? 'M' : 'L'} ${xS(s)} ${yS(T0vals[i])}`).join(' ');
const estLine = grid.map((s, i) => `${i === 0 ? 'M' : 'L'} ${xS(s)} ${yS(TnHat[i])}`).join(' ');
return `
<text x="${margin.left + plotW / 2}" y="${margin.top - 4}" text-anchor="middle" font-size="10" font-weight="bold" fill="#37474f">Estimated Transport Map T̂ₙ</text>
<line x1="${margin.left}" y1="${margin.top}" x2="${margin.left}" y2="${margin.top + plotH}" stroke="#adb5bd"/>
<line x1="${margin.left}" y1="${margin.top + plotH}" x2="${margin.left + plotW}" y2="${margin.top + plotH}" stroke="#adb5bd"/>
<text x="${margin.left + plotW / 2}" y="${margin.top + plotH + 20}" text-anchor="middle" font-size="10" fill="#495057">Base-year income s (thousands)</text>
<text x="${margin.left - 40}" y="${margin.top + plotH / 2}" text-anchor="middle" font-size="10" fill="#495057" transform="rotate(-90,${margin.left - 40},${margin.top + plotH / 2})">Follow-up income T(s)</text>
<path d="${idLine}" fill="none" stroke="#dee2e6" stroke-width="1" stroke-dasharray="4,4"/>
<path d="${trueLine}" fill="none" stroke="#2b8a3e" stroke-width="2.5" stroke-dasharray="6,3"/>
<path d="${estLine}" fill="none" stroke="#1971c2" stroke-width="2.5"/>
<rect x="${margin.left}" y="${margin.top}" width="${plotW}" height="${plotH}" fill="none" stroke="#dee2e6"/>
<g transform="translate(${margin.left + 8},${margin.top + 10})">
<line x1="0" y1="3" x2="15" y2="3" stroke="#dee2e6" stroke-width="1" stroke-dasharray="4,4"/><text x="18" y="7" font-size="9" fill="#495057">Identity</text>
<line x1="0" y1="15" x2="15" y2="15" stroke="#2b8a3e" stroke-width="2.5" stroke-dasharray="6,3"/><text x="18" y="19" font-size="9" fill="#495057">True T₀</text>
<line x1="0" y1="27" x2="15" y2="27" stroke="#1971c2" stroke-width="2.5"/><text x="18" y="31" font-size="9" fill="#495057">Estimated T̂ₙ</text>
</g>`;
})()}
</svg>
</div>
<div style="flex: 1; min-width: 420px;">
<svg width="100%" height="260" viewBox="0 0 440 260" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
const margin = {top: 18, right: 16, bottom: 32, left: 56};
const plotW = 440 - margin.left - margin.right;
const plotH = 260 - margin.top - margin.bottom;
const {uGrid, changeTrue, changeHat} = incTmrResult;
const all = [...changeTrue, ...changeHat, 0];
const yMin = Math.min(...all) - 1, yMax = Math.max(...all) + 1;
const xS = u => margin.left + u * plotW;
const yS = d => margin.top + plotH - (d - yMin) / (yMax - yMin) * plotH;
const zeroLine = `M ${xS(0)} ${yS(0)} L ${xS(1)} ${yS(0)}`;
const trueLine = uGrid.map((u, i) => `${i === 0 ? 'M' : 'L'} ${xS(u)} ${yS(changeTrue[i])}`).join(' ');
const hatLine = uGrid.map((u, i) => `${i === 0 ? 'M' : 'L'} ${xS(u)} ${yS(changeHat[i])}`).join(' ');
const ticks = [0, 0.25, 0.5, 0.75, 1];
return `
<text x="${margin.left + plotW / 2}" y="${margin.top - 4}" text-anchor="middle" font-size="10" font-weight="bold" fill="#37474f">Quantile Change Curve</text>
<line x1="${margin.left}" y1="${margin.top}" x2="${margin.left}" y2="${margin.top + plotH}" stroke="#adb5bd"/>
<line x1="${margin.left}" y1="${margin.top + plotH}" x2="${margin.left + plotW}" y2="${margin.top + plotH}" stroke="#adb5bd"/>
<text x="${margin.left + plotW / 2}" y="${margin.top + plotH + 20}" text-anchor="middle" font-size="10" fill="#495057">Percentile u</text>
<text x="${margin.left - 46}" y="${margin.top + plotH / 2}" text-anchor="middle" font-size="10" fill="#495057" transform="rotate(-90,${margin.left - 46},${margin.top + plotH / 2})">Change in income ($k)</text>
${ticks.map(t => `<text x="${xS(t)}" y="${margin.top + plotH + 12}" text-anchor="middle" font-size="8" fill="#868e96">${(t * 100).toFixed(0)}</text>`).join('')}
<path d="${zeroLine}" fill="none" stroke="#dee2e6" stroke-width="1" stroke-dasharray="4,4"/>
<path d="${trueLine}" fill="none" stroke="#2b8a3e" stroke-width="2.5" stroke-dasharray="6,3"/>
<path d="${hatLine}" fill="none" stroke="#1971c2" stroke-width="2.5"/>
<rect x="${margin.left}" y="${margin.top}" width="${plotW}" height="${plotH}" fill="none" stroke="#dee2e6"/>
<g transform="translate(${margin.left + 8},${margin.top + 10})">
<line x1="0" y1="3" x2="15" y2="3" stroke="#2b8a3e" stroke-width="2.5" stroke-dasharray="6,3"/><text x="18" y="7" font-size="9" fill="#495057">True change</text>
<line x1="0" y1="15" x2="15" y2="15" stroke="#1971c2" stroke-width="2.5"/><text x="18" y="19" font-size="9" fill="#495057">Fitted change</text>
</g>`;
})()}
</svg>
</div>
</div>
<div style="margin-top: 14px;">
<svg width="100%" height="260" viewBox="0 0 900 260" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
const margin = {top: 18, right: 20, bottom: 35, left: 55};
const plotW = 900 - margin.left - margin.right;
const plotH = 260 - margin.top - margin.bottom;
const {xGrid, basePdf, truePdf, obsPdf, fitPdf} = incTmrResult;
const xMin = xGrid[0], xMax = xGrid[xGrid.length - 1];
const yMax = Math.max(...basePdf, ...truePdf, ...obsPdf, ...fitPdf) * 1.15;
const xS = x => margin.left + Math.log(x / xMin) / Math.log(xMax / xMin) * plotW;
const yS = y => margin.top + plotH - (y / yMax) * plotH;
function areaPath(pdf, color, opacity) {
const pts = xGrid.map((x, i) => `${xS(x)},${yS(pdf[i])}`).join(' ');
return `<polygon points="${xS(xMin)},${yS(0)} ${pts} ${xS(xMax)},${yS(0)}" fill="${color}" fill-opacity="${opacity}" stroke="none"/>`;
}
function linePath(pdf, color, width, dash) {
const pts = xGrid.map((x, i) => `${i === 0 ? 'M' : 'L'} ${xS(x)} ${yS(pdf[i])}`).join(' ');
return `<path d="${pts}" fill="none" stroke="${color}" stroke-width="${width}" stroke-dasharray="${dash || 'none'}"/>`;
}
const ticks = [0.5, 1, 2, 5, 10, 20, 50, 100, 200];
return `
<text x="${margin.left + plotW / 2}" y="${margin.top - 4}" text-anchor="middle" font-size="10" font-weight="bold" fill="#37474f">Highlighted Country: Base-Year and Follow-Up Densities</text>
<line x1="${margin.left}" y1="${margin.top}" x2="${margin.left}" y2="${margin.top + plotH}" stroke="#adb5bd"/>
<line x1="${margin.left}" y1="${margin.top + plotH}" x2="${margin.left + plotW}" y2="${margin.top + plotH}" stroke="#adb5bd"/>
<text x="${margin.left - 8}" y="${margin.top - 3}" text-anchor="end" font-size="9" fill="#495057">density</text>
<text x="${margin.left + plotW / 2}" y="${margin.top + plotH + 23}" text-anchor="middle" font-size="10" fill="#495057">Income (thousands, log scale)</text>
${ticks.filter(t => t >= xMin && t <= xMax).map(t => `<text x="${xS(t)}" y="${margin.top + plotH + 13}" text-anchor="middle" font-size="8" fill="#868e96">$${t}k</text>`).join('')}
${areaPath(basePdf, '#1971c2', 0.15)}
${areaPath(obsPdf, '#e03131', 0.12)}
${areaPath(fitPdf, '#2b8a3e', 0.10)}
${linePath(basePdf, '#1971c2', 2.2)}
${linePath(obsPdf, '#e03131', 2.2)}
${linePath(fitPdf, '#2b8a3e', 2.5, '6,3')}
${linePath(truePdf, '#7950f2', 2, '4,4')}
<rect x="${margin.left}" y="${margin.top}" width="${plotW}" height="${plotH}" fill="none" stroke="#dee2e6"/>
<g transform="translate(${margin.left + plotW - 160},${margin.top + 10})">
<line x1="0" y1="3" x2="15" y2="3" stroke="#1971c2" stroke-width="2.2"/><text x="18" y="7" font-size="9" fill="#495057">Base year X</text>
<line x1="0" y1="15" x2="15" y2="15" stroke="#7950f2" stroke-width="2" stroke-dasharray="4,4"/><text x="18" y="19" font-size="9" fill="#495057">True follow-up</text>
<line x1="0" y1="27" x2="15" y2="27" stroke="#e03131" stroke-width="2.2"/><text x="18" y="31" font-size="9" fill="#495057">Observed follow-up</text>
<line x1="0" y1="39" x2="15" y2="39" stroke="#2b8a3e" stroke-width="2.5" stroke-dasharray="6,3"/><text x="18" y="43" font-size="9" fill="#495057">Fitted follow-up</text>
</g>`;
})()}
</svg>
</div>
<div style="margin-top:12px;padding:12px;background:#f8f9fa;border-radius:6px;font-size:0.92em;">
<b>Descriptive Wasserstein R² = ${incTmrResult.w2R2.toFixed(3)}</b> |
<b>‖T̂ₙ−T₀‖²<sub>L²(Qₙ)</sub> = ${incTmrResult.ise.toFixed(3)}</b> |
γ = ${incTmrResult.gamma.toFixed(2)} |
σ = ${incTmrResult.sigma.toFixed(2)} |
n = ${incTmrResult.n} |
highlighted: ${incTmrResult.selectedLabel} (index ${incTmrResult.selectedIdx})
</div>
<p style="margin-top:8px;font-size:0.85em;color:#495057;">
<b>Top left:</b> The true regression map T₀ (green dashed) and the PAVA estimate T̂ₙ (blue solid). The gray dotted line is the identity — deviations show how base-year income levels are mapped to follow-up levels across the common support.<br>
<b>Top right:</b> Quantile change curve for the highlighted country: fitted vs true change in income at each percentile. This is the application analogue of a growth-incidence curve.<br>
<b>Bottom:</b> Base-year density (blue), true follow-up (purple dashed), observed follow-up with noise (red), and the fitted follow-up distribution from T̂ₙ (green dashed).
</p>
</div>`
```
::: {.callout-tip title="Try these experiments"}
- **Increase the number of countries:** The PAVA estimate $\hat T_n$ generally moves closer to $T_0$, although a single realization can fluctuate. The cube-root rate means that, asymptotically, an eightfold increase in $n$ halves the order of the $L^2$ error.
- **Increase the noise $\sigma$:** Pair-specific maps scatter more widely around $T_0$, and the descriptive Wasserstein $R^2$ drops.
- **Change $\gamma$:** Positive values lift the middle of the distribution relative to the tails (a "growth from the middle" pattern); negative values do the opposite. Watch how the quantile-change curve and the fitted follow-up density respond.
- **Switch highlighted country:** Different base-year profiles change the shape of the predicted follow-up distribution, but the estimated common map is shared across all countries.
:::
## Key Takeaways
- **Two strategies** for distributional regression: transformation-based (map to Hilbert space, use FDA) and intrinsic (work directly in Wasserstein geometry). They are complementary — choose based on whether the Wasserstein metric is meaningful for the application.
- **Transformation-based approaches** [@PetersenMueller2016; @KokoszkaEtAl2019] use transformations such as LQD or log hazard to work in $L^2$ and then map results back to densities. The induced Hilbert norm is generally not $W_2$.
- **Wasserstein approaches** use different predictor structures. @PetersenMueller2019 use Euclidean predictors and metric responses; @ChenLinMueller2023 use univariate distribution predictors and estimate a linear operator between tangent spaces; @GhodratiPanaretos2022 estimate one common monotone map acting directly on univariate distribution predictors.
- **Transport-map regression** models $Y_i = \epsilon_i\#(T_0\#X_i)$ with monotone $T_0$. Estimation proceeds by weighted isotonic regression (PAVA), minimizing $\sum_i \int |T(s) - T_{X_i \to Y_i}(s)|^2 dX_i(s)$. The criterion is strictly convex; the solution is unique.
- **Convergence rate** is $n^{-1/3}$ in $L^2(Q)$ for fully observed distributions under the paper's assumptions. The paper gives an additional distribution-estimation term under partial observation; it does not claim minimax optimality.
- **Interpretation requires care:** $T_0$ maps support values and therefore maps predictor quantiles after composition with $Q_X$. Quantile coupling describes marginal-distribution change, not individual mobility.
## Exercises
1. **Optimal transport map between Gaussians.** Let $X_i \sim N(\mu_i, \sigma_i^2)$ and $Y_i \sim N(\nu_i, \tau_i^2)$. Derive the pair-specific optimal transport map $T_{X_i \to Y_i}(s)$ in closed form. Show that it is linear with slope $\tau_i / \sigma_i$ and intercept $\nu_i - \mu_i \tau_i / \sigma_i$. <a href="javascript:void(0)" onclick="showSolution('l14-sol-1')" class="solution-link">📝 Show Solution</a>
2. **Monotonicity and model validity.** Why must $T_0$ be nondecreasing in the transport-map regression framework? What would go wrong if the estimated $\hat{T}_n$ were not monotone? <a href="javascript:void(0)" onclick="showSolution('l14-sol-2')" class="solution-link">📝 Show Solution</a>
3. **Choosing approaches.** Compare transformation-based (Petersen & Müller 2016), tangent-space (Chen, Lin, Müller 2023), and transport-map (Ghodrati & Panaretos 2022) regression for a problem where the response is a distribution on $\mathbb{R}^2$. Which approach would you choose and why? <a href="javascript:void(0)" onclick="showSolution('l14-sol-3')" 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-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="l14-sol-1"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 1: OT Map Between Gaussians</h4><button onclick="closeSolution('l14-sol-1')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Derive $T_{X_i \to Y_i}(s)$ for Gaussians.</div><div class="solution-answer"><strong>Solution:</strong> For $X_i \sim N(\mu_i, \sigma_i^2)$ and $Y_i \sim N(\nu_i, \tau_i^2)$, the CDFs are $F_{X_i}(s) = \Phi((s - \mu_i)/\sigma_i)$ and $F_{Y_i}^{-1}(u) = \nu_i + \tau_i \Phi^{-1}(u)$. Then
$$T_{X_i \to Y_i}(s) = F_{Y_i}^{-1}(F_{X_i}(s)) = \nu_i + \tau_i \Phi^{-1}\!\left(\Phi\!\left(\frac{s - \mu_i}{\sigma_i}\right)\right) = \nu_i + \frac{\tau_i}{\sigma_i}(s - \mu_i) = \frac{\tau_i}{\sigma_i}s + \left(\nu_i - \mu_i\frac{\tau_i}{\sigma_i}\right).$$
This is linear: the slope $\tau_i/\sigma_i$ captures the scale change, and the intercept adjusts for the mean shift. For equal variances ($\sigma_i = \tau_i$), the map is a pure translation: $T(s) = s + (\nu_i - \mu_i)$.</div></div></dialog>
<dialog id="l14-sol-2"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 2: Monotonicity and Model Validity</h4><button onclick="closeSolution('l14-sol-2')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Why must $T_0$ be nondecreasing?</div><div class="solution-answer"><strong>Solution:</strong> Any measurable $T$ produces a valid pushforward $T\#\mu$. The stronger requirement is that the same map represent the one-dimensional quadratic-cost <em>optimal</em> transport from every predictor $\mu$ to its fitted response. That map is the nondecreasing rearrangement, uniquely $\mu$-almost everywhere when $\mu$ is absolutely continuous. Without monotonicity, $T$ would still define a distribution but would no longer be the model's optimal deformation; the quantile identity $Q_{T\#\mu}=T\circ Q_\mu$ and the reduction of the loss to weighted isotonic regression would also fail. PAVA enforces the structural constraint.</div></div></dialog>
<dialog id="l14-sol-3"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 3: Choosing Approaches in Higher Dimensions</h4><button onclick="closeSolution('l14-sol-3')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Which approach for distributions on $\mathbb{R}^2$?</div><div class="solution-answer"><strong>Solution:</strong> None of the three papers, as developed, directly solves this $\mathbb R^2$ problem. Petersen and Müller (2016) study transformations of one-dimensional densities. Ghodrati and Panaretos (2022) rely on the one-dimensional nondecreasing rearrangement and PAVA. Chen, Lin, and Müller (2023) also develop and analyze their method for univariate distributions; they note that extension is possible in principle, but multivariate optimal maps lack the quantile closed form, their parallel-transport expression no longer applies, and computation and asymptotic analysis must change. One should therefore use a method specifically developed for multivariate distribution responses or predictors (for example, a multivariate or sliced-Wasserstein regression method), rather than claiming that one of these three implementations transfers unchanged.</div></div></dialog>
<dialog id="l14-sol-4"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 4: The $n^{-1/3}$ Rate</h4><button onclick="closeSolution('l14-sol-4')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Why $n^{-1/3}$ for isotonic regression?</div><div class="solution-answer"><strong>Solution:</strong> A useful heuristic approximates a sufficiently regular monotone curve by $K$ constant blocks. The integrated squared approximation error is of order $K^{-2}$, while estimating $K$ block levels contributes variance of order $K/n$. Balancing $K^{-2}$ and $K/n$ gives $K\asymp n^{1/3}$ and squared $L^2$ error $n^{-2/3}$, hence $L^2$ error $n^{-1/3}$. The paper's proof uses empirical-process entropy for its monotone-map class. This heuristic should not be read as a theorem that smoothness alone yields a root-$n$ nonparametric rate; the paper neither proves a minimax lower bound nor a faster smooth estimator.</div></div></dialog>
<dialog id="l14-sol-5"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 5: Income-Distribution Interpretation</h4><button onclick="closeSolution('l14-sol-5')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Interpret opposite fitted changes below and above the median.</div><div class="solution-answer"><strong>Solution:</strong> The fitted lower quantiles move upward and the fitted upper quantiles move downward, so the predicted distribution is compressed in this quantile-by-quantile sense. In particular, fitted gaps between a lower and an upper quantile shrink. Whether a named inequality index decreases should be checked from the entire fitted distribution and the index's scale convention; the sign pattern alone is not a universal theorem for every inequality measure. Most importantly, this comparison uses marginal quantiles. Without linked individual data, it says nothing about who moved, rank transitions, or individual income mobility.</div></div></dialog>
<script>
function showSolution(id) { const d = document.getElementById(id); if (d) { d.showModal(); if (window.MathJax && MathJax.typesetPromise) { MathJax.typesetPromise([d]).catch(function(e) { console.log('MathJax error:', e); }); } } }
function closeSolution(id) { const d = document.getElementById(id); if (d) d.close(); }
document.addEventListener('click', function(e) { if (e.target.tagName === 'DIALOG') e.target.close(); });
</script>
## Further Reading
- @PetersenMueller2016 — Foundational transformation-based approach: functional data analysis for density functions via log-quantile-density transform.
- @PetersenMueller2019 — Global and local Fréchet regression for metric-space responses with Euclidean predictors.
- @ChenLinMueller2023 — Univariate tangent-space framework for distribution-to-scalar and distribution-to-distribution regression.
- @GhodratiPanaretos2022 — Distribution-on-distribution regression via optimal transport maps with monotonicity constraints (the focus of this lecture).
- @ZhouMueller2024 — Distribution-response regression with Euclidean covariates and empirical measures, designed for unequal and sometimes sparse within-distribution sample sizes.
- @PetersenLiuDivani2021 — Wasserstein $F$-tests and simultaneous confidence bands for global Fréchet regression of density responses on vector predictors.
- @PanaretosZemel2020 — Comprehensive textbook on statistics and optimal transport, covering both theoretical foundations and methodological developments.
- @AyerEtAl1955 — An early source of the pooling construction that became the pool-adjacent-violators algorithm.
## Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
👉 **[Lecture 14 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-14-quiz.qmd)**