Distributional data analysis

Distributional data analysis is a branch of nonparametric statistics that is related to functional data analysis. It is concerned with random objects that are probability distributions, i.e., the statistical analysis of samples of random distributions where each atom of a sample is a distribution. One of the main challenges in distributional data analysis is that although the space of probability distributions is a convex space, it is not a vector space.

Notation

Let $\nu$ be a probability measure on $D$, where $D\subset \mathbb {R} ^{p}$ with $p\geq 1$. The probability measure $\nu$ can be equivalently characterized as cumulative distribution function $F$ or probability density function $f$ if it exists. For univariate distributions with $p=1$, quantile function $Q=F^{-1}$ can also be used.

Let ${\mathcal {F}}$ be a space of distributions $\nu$ and let $d$ be a metric on ${\mathcal {F}}$ so that $({\mathcal {F}},d)$ forms a metric space. There are various metrics available for $d$. For example, suppose $\nu _{1},\;\nu _{2}\in {\mathcal {F}}$, and let $f_{1}$ and $f_{2}$ be the density functions of $\nu _{1}$ and $\nu _{2}$, respectively. The Fisher-Rao metric is defined as $$d_{FR}(f_{1},f_{2})=\arccos \left(\int _{D}{\sqrt {f_{1}(x)f_{2}(x)}}dx\right).$$

For univariate distributions, let $Q_{1}$ and $Q_{2}$ be the quantile functions of $\nu _{1}$ and $\nu _{2}$. Denote the $L^{p}$-Wasserstein space as ${\mathcal {W}}_{p}$, which is the space of distributions with finite $p$-th moments. Then, for $\nu _{1},\;\nu _{2}\in {\mathcal {W}}_{p}$, the $L^{p}$-Wasserstein metric is defined as $$d_{W_{p}}(\nu _{1},\nu _{2})=\left(\int _{0}^{1}[Q_{1}(s)-Q_{2}(s)]^{p}ds\right)^{1/p}.$$

Mean and variance

For a probability measure $\nu \in {\mathcal {F}}$, consider a random process ${\mathfrak {F}}$ such that $\nu \sim {\mathfrak {F}}$. One way to define mean and variance of $\nu$ is to introduce the Fréchet mean and the Fréchet variance. With respect to the metric $d$ on ${\mathcal {F}}$, the Fréchet mean $\mu _{\oplus }$, also known as the barycenter, and the Fréchet variance $V_{\oplus }$ are defined as $${\begin{aligned}\mu _{\oplus }&=\operatorname {argmin} _{\mu \in {\mathcal {F}}}\mathbb {E} [d^{2}(\nu ,\mu )],\\V_{\oplus }&=\mathbb {E} [d^{2}(\nu ,\mu _{\oplus })].\end{aligned}}$$

A widely used example is the Wasserstein-Fréchet mean, or simply the Wasserstein mean, which is the Fréchet mean with the $L^{2}$-Wasserstein metric $d_{W_{2}}$. For $\nu ,\;\mu \in {\mathcal {W}}_{2}$, let $Q_{\nu },\;Q_{\mu }$ be the quantile functions of $\nu$ and $\mu$, respectively. The Wasserstein mean and Wasserstein variance is defined as $${\begin{aligned}\mu _{\oplus }^{*}&=\operatorname {argmin} _{\mu \in {\mathcal {W}}_{2}}\mathbb {E} \left[\int _{0}^{1}(Q_{\nu }(s)-Q_{\mu }(s))^{2}ds\right],\\V_{\oplus }^{*}&=\mathbb {E} \left[\int _{0}^{1}(Q_{\nu }(s)-Q_{\mu _{\oplus }^{*}}(s))^{2}ds\right].\end{aligned}}$$

Modes of variation

Modes of variation are useful concepts in depicting the variation of data around the mean function. Based on the Karhunen-Loève representation, modes of variation show the contribution of each eigenfunction to the mean.

Functional principal component analysis

Functional principal component analysis (FPCA) can be directly applied to the probability density functions. Consider a distribution process $\nu \sim {\mathfrak {F}}$ and let $f$ be the density function of $\nu$. Let the mean density function as $\mu (t)=\mathbb {E} \left[f(t)\right]$ and the covariance function as $G(s,t)=\operatorname {Cov} (f(s),f(t))$ with orthonormal eigenfunctions $\{\phi _{j}\}_{j=1}^{\infty }$ and eigenvalues $\{\lambda _{j}\}_{j=1}^{\infty }$.

By the Karhunen-Loève theorem, $f(t)=\mu (t)+\sum _{j=1}^{\infty }\xi _{j}\phi _{j}(t)$, where principal components $\xi _{j}=\int _{D}[f(t)-\mu (t)]\phi _{j}(t)dt$. The $j$th mode of variation is defined as $g_{j}(t,\alpha )=\mu (t)+\alpha {\sqrt {\lambda _{j}}}\phi _{j}(t),\quad t\in D,\;\alpha \in [-A,A]$ with some constant $A$, such as 2 or 3.

Transformation FPCA

Assume the probability density functions $f$ exist, and let ${\mathcal {F}}_{f}$ be the space of density functions. Transformation approaches introduce a continuous and invertible transformation $\Psi :{\mathcal {F}}_{f}\to \mathbb {H}$, where $\mathbb {H}$ is a Hilbert space of functions. For instance, the log quantile density transformation or the centered log ratio transformation are popular choices.

For $f\in {\mathcal {F}}_{f}$, let $Y=\Psi (f)$, the transformed functional variable. The mean function $\mu _{Y}(t)=\mathbb {E} \left[Y(t)\right]$ and the covariance function $G_{Y}(s,t)=\operatorname {Cov} (Y(s),Y(t))$ are defined accordingly, and let $\{\lambda _{j},\phi _{j}\}_{j=1}^{\infty }$ be the eigenpairs of $G_{Y}(s,t)$. The Karhunen-Loève decomposition gives $Y(t)=\mu _{Y}(t)+\sum _{j=1}^{\infty }\xi _{j}\phi _{j}(t)$, where $\xi _{j}=\int _{D}[Y(t)-\mu _{Y}(t)]\phi _{j}(t)dt$. Then, the $j$th transformation mode of variation is defined as $g_{j}^{TF}(t,\alpha )=\Psi ^{-1}\left(\mu _{Y}+\alpha {\sqrt {\lambda _{j}}}\phi _{j}\right)(t),\quad t\in D,\;\alpha \in [-A,A].$

Log FPCA and Wasserstein Geodesic PCA

Endowed with metrics such as the Wasserstein metric $d_{W_{2}}$ or the Fisher-Rao metric $d_{FR}$, we can employ the (pseudo) Riemannian structure of ${\mathcal {F}}$. Denote the tangent space at the Fréchet mean $\mu _{\oplus }$ as $T_{\mu _{\oplus }}$, and define the logarithm and exponential maps $\log _{\mu _{\oplus }}:{\mathcal {F}}\to T_{\mu _{\oplus }}$ and $\exp _{\mu _{\oplus }}:T_{\mu _{\oplus }}\to {\mathcal {F}}$. Let $Y$ be the projected density onto the tangent space, $Y=\log _{\mu _{\oplus }}(f)$.

In Log FPCA, FPCA is performed to $Y$ and then projected back to ${\mathcal {F}}$ using the exponential map. Therefore, with $Y(t)=\mu _{Y}(t)+\sum _{j=1}^{\infty }\xi _{j}\phi _{j}(t)$, the $j$th Log FPCA mode of variation is defined as $g_{j}^{Log}(t,\alpha )=\exp _{f_{\oplus }}\left(\mu _{f_{\oplus }}+\alpha {\sqrt {\lambda _{j}}}\phi _{j}\right)(t),\quad t\in D,\;\alpha \in [-A,A].$

As a special case, consider $L^{2}$-Wasserstein space ${\mathcal {W}}_{2}$, a random distribution $\nu \in {\mathcal {W}}_{2}$, and a subset $G\subset {\mathcal {W}}_{2}$. Let $d_{W_{2}}(\nu ,G)=\inf _{\mu \in G}d_{W_{2}}(\nu ,\mu )$ and $K_{W_{2}}(G)=\mathbb {E} \left[d_{W_{2}}^{2}(\nu ,G)\right]$. Let ${\text{CL}}({\mathcal {W}}_{2})$ be the metric space of nonempty, closed subsets of ${\mathcal {W}}_{2}$, endowed with Hausdorff distance, and define $\operatorname {CG} _{\nu _{0},k}({\mathcal {W}}_{2})=\{G\in \operatorname {CL} ({\mathcal {W}}_{2}):\nu _{0}\in G,G{\text{ is a geodesic set s.t. }}\operatorname {dim} (G)\leq k\},\;k\geq 1.$ Let the reference measure $\nu _{0}$ be the Wasserstein mean $\mu _{\oplus }$. Then, a principal geodesic subspace (PGS) of dimension $k$ with respect to $\mu _{\oplus }$ is a set $G_{k}=\operatorname {argmin} _{G\in {\text{CG}}_{\nu _{\oplus },k}({\mathcal {W}}_{2})}K_{W_{2}}(G)$.

Note that the tangent space $T_{\mu _{\oplus }}$ is a subspace of $L_{\mu _{\oplus }}^{2}$, the Hilbert space of ${\mu _{\oplus }}$-square-integrable functions. Obtaining the PGS is equivalent to performing PCA in $L_{\mu _{\oplus }}^{2}$ under constraints to lie in the convex and closed subset. Therefore, a simple approximation of the Wasserstein Geodesic PCA is the Log FPCA by relaxing the geodesicity constraint, while alternative techniques are suggested.

Distributional regression

Fréchet regression

Fréchet regression is a generalization of regression with responses taking values in a metric space and Euclidean predictors. Using the Wasserstein metric $d_{W_{2}}$, Fréchet regression models can be applied to distributional objects. The global Wasserstein-Fréchet regression model is defined as

$${\begin{aligned}m_{\oplus }(x)&=\operatorname {argmin} _{\omega \in {\mathcal {F}}}\mathbb {E} \left[s_{G}(X,x)d_{W_{2}}^{2}(\nu ,\omega )\right],\\s_{G}(X,x)&=1+(X-\mathbb {E} [X])^{\top }{\text{Var}}(X)^{-1}(x-\mathbb {E} [X]),\end{aligned}}$$ (1)

which generalizes the standard linear regression.

For the local Wasserstein-Fréchet regression, consider a scalar predictor $X\in \mathbb {R}$ and introduce a smoothing kernel $K_{h}(\cdot )=h^{-1}K(\cdot /h)$. The local Fréchet regression model, which generalizes the local linear regression model, is defined as $${\begin{aligned}l_{\oplus }(x)&=\operatorname {argmin} _{\omega \in {\mathcal {F}}}\mathbb {E} \left[s_{L}(X,x,h)d_{W_{2}}^{2}(\nu ,\omega )\right],\\s_{L}(X,x,h)&=\sigma _{0}^{-2}\{K_{h}(X-x)[\mu _{2}-\mu _{1}(X-x)]\},\end{aligned}}$$ where $\mu _{j}=\mathbb {E} \left[K_{h}(X-x)(X-x)^{j}\right]$, $j=0,1,2,$ and $\sigma _{0}^{2}=\mu _{0}\mu _{2}-\mu _{1}^{2}$.

Transformation based approaches

Consider the response variable $\nu$ to be probability distributions. With the space of density functions ${\mathcal {F}}_{f}$ and a Hilbert space of functions $\mathbb {H}$, consider continuous and invertible transformations $\Psi :{\mathcal {F}}_{f}\to \mathbb {H}$. Examples of transformations include log hazard transformation, log quantile density transformation, or centered log-ratio transformation. Linear methods such as functional linear models are applied to the transformed variables. The fitted models are interpreted back in the original density space ${\mathcal {F}}$ using the inverse transformation.

Random object approaches

In Wasserstein regression, both predictors $\omega$ and responses $\nu$ can be distributional objects. Let $\omega {\oplus }$ and $\nu _{\oplus }$ be the Wasserstein mean of $\omega$ and $\nu$, respectively. The Wasserstein regression model is defined as $$\mathbb {E} (\log _{\nu _{\oplus }}\nu |\log _{\omega {\oplus }}\omega )=\Gamma (\log _{\omega {\oplus }}\omega ),$$ with a linear regression operator $$\Gamma g(t)=\langle \beta (\cdot ,t),g\rangle _{\omega {\oplus }},\;t\in D,\;g\in T_{\omega {\oplus }},\;\beta :D^{2}\to \mathbb {R} .$$ Estimation of the regression operator is based on empirical estimators obtained from samples. Also, the Fisher-Rao metric $d_{FR}$ can be used in a similar fashion.

Hypothesis testing

Wasserstein F-test

Wasserstein $F$-test has been proposed to test for the effects of the predictors in the Fréchet regression framework with the Wasserstein metric. Consider Euclidean predictors $X\in \mathbb {R} ^{p}$ and distributional responses $\nu \in {\mathcal {W}}_{2}$. Denote the Wasserstein mean of $\nu$ as $\mu _{\oplus }^{*}$, and the sample Wasserstein mean as ${\hat {\mu }}_{\oplus }^{*}$. Consider the global Wasserstein-Fréchet regression model $m_{\oplus }(x)$ defined in (1), which is the conditional Wasserstein mean given $X=x$. The estimator of $m_{\oplus }(x)$, ${\hat {m}}_{\oplus }(x)$ is obtained by minimizing the empirical version of the criterion.

Let $F$, $Q$, $f$, $F_{\oplus }^{*}$, $Q_{\oplus }^{*}$, $f_{\oplus }^{*}$, $F_{\oplus }(x)$, $Q_{\oplus }(x)$, and $f_{\oplus }(x)$ denote the cumulative distribution, quantile, and density functions of $\nu$, $\mu _{\oplus }^{*}$, and $m_{\oplus }(x)$, respectively. For a pair $(X,\nu )$, define $T=Q\circ F_{\oplus }(X)$ be the optimal transport map from $m_{\oplus }(X)$ to $\nu$. Also, define $S=Q_{\oplus }(X)\circ F_{\oplus }^{*}$, the optimal transport map from $\mu _{\oplus }^{*}$ to $m_{\oplus }(x)$. Finally, define the covariance kernel $K(u,v)=\mathbb {E} [{\text{Cov}}((T\circ S)(u),(T\circ S)(v))]$ and by the Mercer decomposition, $K(u,v)=\sum _{j=1}^{\infty }\lambda _{j}\phi _{j}(u)\phi _{j}(v)$.

If there are no regression effects, the conditional Wasserstein mean would equal the Wasserstein mean. That is, hypotheses for the test of no effects are $$H_{0}:m_{\oplus }(x)\equiv \mu _{\oplus }^{*}\quad {\text{v.s.}}\quad H_{1}:{\text{Not }}H_{0}.$$ To test for these hypotheses, the proposed global Wasserstein $F$-statistic and its asymptotic distribution are $$F_{G}=\sum _{i=1}^{n}d_{W_{2}}^{2}({\hat {m}}_{\oplus }(x),{\hat {\mu }}_{\oplus }^{*}),\quad F_{G}|X_{1},\cdots ,X_{n}{\overset {d}{\longrightarrow }}\sum _{j=1}^{\infty }\lambda _{j}V_{j}\;a.s.,$$ where $V_{j}{\overset {iid}{\sim }}\chi _{p}^{2}$. An extension to hypothesis testing for partial regression effects, and alternative testing approximations using the Satterthwaite's approximation or a bootstrap approach are proposed.

Tests for the intrinsic mean

The Hilbert sphere ${\mathcal {S}}^{\infty }$ is defined as ${\mathcal {S}}^{\infty }=\left\{f\in \mathbb {H} :\|f\|_{\mathbb {H} }=1\right\}$, where $\mathbb {H}$ is a separable infinite-dimensional Hilbert space with inner product $\langle \cdot ,\cdot \rangle _{\mathbb {H} }$ and norm $\|\cdot \|_{\mathbb {H} }$. Consider the space of square root densities ${\mathcal {X}}=\left\{x:D\to \mathbb {R} :x={\sqrt {f}},\int _{D}f(t)dt=1\right\}$. Then with the Fisher-Rao metric $d_{FR}$ on $f$, ${\mathcal {X}}$ is the positive orthant of the Hilbert sphere ${\mathcal {S}}^{\infty }$ with $\mathbb {H} =L^{2}(D)$.

Let a chart $\tau :U\subset {\mathcal {S}}^{\infty }\to \mathbb {G}$ as a smooth homeomorphism that maps $U$ onto an open subset $\tau (U)$ of a separable Hilbert space $\mathbb {G}$ for coordinates. For example, $\tau$ can be the logarithm map.

Consider a random element $x={\sqrt {f}}\in {\mathcal {X}}$ equipped with the Fisher-Rao metric, and write its Fréchet mean as $\mu$. Let the empirical estimator of $\mu$ using $n$ samples as ${\hat {\mu }}$. Then central limit theorem for ${\hat {\mu }}_{\tau }=\tau ({\hat {\mu }})$ and $\mu _{\tau }=\tau (\mu )$ holds: ${\sqrt {n}}({\hat {\mu }}_{\tau }-\mu _{\tau }){\overset {L}{\longrightarrow }}Z,\;n\to \infty$, where $Z$ is a Gaussian random element in $\mathbb {G}$ with mean 0 and covariance operator ${\mathcal {T}}$. Let the eigenvalue-eigenfunction pairs of ${\mathcal {T}}$ and the estimated covariance operator ${\hat {\mathcal {T}}}$ as $(\lambda _{k},\phi _{k})_{k=1}^{\infty }$ and $({\hat {\lambda }}_{k},{\hat {\phi }}_{k})_{k=1}^{\infty }$, respectively.

Consider one-sample hypothesis testing $$H_{0}:\mu =\mu _{0}\quad {\text{v.s.}}\quad H_{1}:\mu \neq \mu _{0},$$ with $\mu _{0}\in {\mathcal {S}}^{\infty }$. Denote $\|\cdot \|_{\mathbb {G} }$ and $\langle \cdot ,\cdot \rangle _{\mathbb {G} }$ as the norm and inner product in $\mathbb {G}$. The test statistics and their limiting distributions are $${\begin{aligned}T_{1}&=n\|\tau ({\hat {\mu }})-\tau (\mu _{0})\|_{\mathbb {G} }^{2}{\overset {L}{\longrightarrow }}\lambda _{k}W_{k},\\S_{1}&=n\sum _{k=1}^{K}{\frac {\langle \tau ({\hat {\mu }})-\tau (\mu _{0}),{\hat {\phi }}_{k}\rangle _{\mathbb {G} }^{2}}{{\hat {\lambda }}_{k}}}{\overset {L}{\longrightarrow }}\chi _{K}^{2},\end{aligned}}$$ where $W_{k}{\overset {iid}{\sim }}\chi _{1}^{2}$. The actual testing procedure can be done by employing the limiting distributions with Monte Carlo simulations, or bootstrap tests are possible. An extension to the two-sample test and paired test are also proposed.

Distributional time series

Autoregressive (AR) models for distributional time series are constructed by defining stationarity and utilizing the notion of difference between distributions using $d_{W_{2}}$ and $d_{FR}$.

In Wasserstein autoregressive model (WAR), consider a stationary density time series $f_{t}$ with Wasserstein mean $f_{\oplus }$. Denote the difference between $f_{t}$ and $f_{\oplus }$ using the logarithm map, $f_{t}\ominus f_{\oplus }=\log _{f_{\oplus }}f_{t}=T_{t}-{\text{id}}$, where $T_{t}=Q_{t}\circ F_{\oplus }$ is the optimal transport from $f_{\oplus }$ to $f_{t}$ in which $F_{t}$ and $F_{\oplus }$ are the cdf of $f_{t}$ and $f_{\oplus }$. An $AR(1)$ model on the tangent space $T_{f_{\oplus }}$ is defined as $V_{t}=\beta V_{t-1}+\epsilon _{t},\;t\in \mathbb {Z} ,$ for $V_{t}\in T_{f_{\oplus }}$ with the autoregressive parameter $\beta \in \mathbb {R}$ and mean zero random i.i.d. innovations $\epsilon _{t}$. Under proper conditions, $\mu _{t}=\exp _{f_{\oplus }}(V_{t})$ with densities $f_{t}$ and $V_{t}=\log _{f_{\oplus }}(\mu _{t})$. Accordingly, $WAR(1)$, with a natural extension to order $p$, is defined as $$T_{t}-{\text{id}}=\beta (T_{t-1}-{\text{id}})+\epsilon _{t}.$$

On the other hand, the spherical autoregressive model (SAR) considers the Fisher-Rao metric. Following the settings of ##Tests for the intrinsic mean, let $x_{t}\in {\mathcal {X}}$ with Fréchet mean $\mu _{x}$. Let $\theta =\arccos(\langle x_{t},\mu _{x}\rangle )$, which is the geodesic distance between $x_{t}$ and $\mu _{x}$. Define a rotation operator $Q_{x_{t},\mu _{x}}$ that rotates $x_{t}$ to $\mu _{x}$. The spherical difference between $x_{t}$ and $\mu _{x}$ is represented as $R_{t}=x_{t}\ominus \mu _{x}=\theta Q_{x_{t},\mu _{x}}$. Assume that $R_{t}$ is a stationary sequence with the Fréchet mean $\mu _{R}$, then $SAR(1)$ is defined as $$R_{t}-\mu _{R}=\beta (R_{t-1}-\mu _{R})+\epsilon _{t},$$ where $\mu _{R}=\mathbb {E} R_{t}$ and mean zero random i.i.d innovations $\epsilon _{t}$. An alternative model, the differenced based spherical autoregressive (DSAR) model is defined with $R_{t}=x_{t+1}\ominus x_{t}$, with natural extensions to order $p$. A similar extension to the Wasserstein space was introduced.