Goal

This is the report from my point of view on the paper Functional Data Analysis for Density Functions by Transformation to a Hilbert Space (Petersen and Müller 2016).

1. Functional data that are nonnegative and have a constrained integral can be considered as samples of one-dimensional density function.

2. The density function is nonnegative and has an integral of 1, hence commonly used Hilbert space based methods of functional data analysis are not applicable.

3. To address this issue, the paper proposes a new method which maps the densities to a Hilbert space, where the map is a continuous and invertible function.

4. Basic operations are implemented in the Hilbert space, then apply the inverse transformation to the density space.

Introduction

The functional modeling of density functions is difficult due to the two constrains $\int f(x)dx=1$ and $f\ge 0$ which means that the densities space is not linear.

The key idea is to map probability densities into a linear function space by using a suitably chosen continuous and invertible map $\psi$.

Preliminaries

Density modelling

Assume that data consist of a sample of $n$ (random) density functions $f_1,\dots ,f_n$, where the densities are supported on a common interval $[0,T]$ for some $T>0$.

Denote the space of continuous and strictly positive densities on $[0,1]$ by $\mathcal{G}$.

The sample consists of i.i.d. realizations of an underlying stochastic process, that is, each density is independently distributed as $f\sim \mathfrak{F}$, where $\mathfrak{F}$ is an $L^2$ process on $[0,1]$ taking values in some space $\mathcal{F}\subset \mathcal{G}$

A basic assumption we make on the space $\mathcal{F}$ is:

Densities $f$ can equivalently be represented as:

1. cumulative distribution functions (c.d.f.) $F$ with domain $[0,1]$, hazard functions $h=f/(1-F)$ (possibly on a subdomain of $[0,1]$ where $F(x)<1$).
2. quantile functions $Q=F^{-1}$, with support $[0,1]$.
3. Occasionally of interest is the equivalent notion of the quantile-density function $q(t)=Q^{\prime }(t)=\frac{d}{dt}{F}^{-1}(t)=[f(Q(t)){]}^{-1}$, from which we obtain $f(x)=[q(F(x)){]}^{-1}$.

All of these functions provide equivalent characterizations of distributions.

In many situations, the densities themselves will not be directly observed. Instead, for each $i$, we may observe an i.i.d. sample of data $W_{il}$, $l=1,\dots ,N_i$, that are generated by the random density $f_i$.

Thus, there are two random mechanisms at work that are assumed to be independent: the first generates the sample of densities and the second generates the samples of real-valued random data; one sample for each random density in the sample of densities.

Hence, the probability space can be thought of as a product space $({\mathrm{\Omega }}_1\times {\mathrm{\Omega }}_2,\mathcal{A},P)$, where $P=P_1\otimes P_2$.

Metrics in the space of density functions

In previous applied and methodological work, it was found that a metric $d_Q$ based on quantile functions $d_Q(f,g{)}^2={\int }_0^1(F^{-1}(t)-G^{-1}(t){)}^2\mathrm{d}t$ is particularly promising from a practical point of view.

This quantile metric has connections to the optimal transport problem , and corresponds to the Wasserstein metric between two probability measures,

$$\begin{array}{}\text{(1)}& d_W(f,g{)}^2=\underset{X\sim f,Y\sim g}{inf}E(X-Y{)}^2,\end{array}$$

where the expectation is with respect to the joint distribution of $(X,Y)$. The equivalence $d_Q=d_W$ can be most easily seen by applying a covariance identity, details can be found in the supplementary material.

We will develop our methodology for a general metric, which will be denoted by $d$ in the following, and may stand for any of the above metrics in the space of densities.

Density estimation

The densities themselves must first be estimated.

Consider the estimation of a density $f\in \mathcal{F}$ from an i.i.d. sample (generated by $f$) of size $N$ by an estimator $\check{f}$. Here, $N=N(n)$ will implicitly represent a sequence that depends on $n$, the size of the sample of random densities.

For the theoretical results, a density estimator $\check{f}$ must satisfy the following consistency properties in terms of the $L^2$ and uniform metrics (denoted as $d_2$ and $d_{\mathrm{\infty }}$, resp.):

1. For a sequence $b_N=o(1)$, the density estimator $\check{f}$, based on an i.i.d. sample of size $N$, satisfies $\check{f}\ge 0$, ${\int }_0^1\check{f}(x)dx=1$ and

$$\begin{array}{}\text{(2)}& \underset{f\in \mathcal{F}}{sup}E(d_2(f,\check{f}{)}^2)=O(b_N^2).\end{array}$$

1. For a sequence $a_N=o(1)$ and some $R>0$, the density estimator $\check{f}$, based on an i.i.d. sample of size $N$, satisfies

$$\begin{array}{}\text{(3)}& \underset{f\in \mathcal{F}}{sup}P(d_{\mathrm{\infty }}(f,\check{f})>R{a}_N)\to 0.\end{array}$$

Therefore, we introduce here a modified density estimator of kernel type that is guaranteed to satisfy $\text{(2)}$ and $\text{(3)}$.

Let $$ be a kernel that corresponds to a continuous probability density function and $h<1/2$ be the bandwidth. We define a new kernel density estimator to estimate the density $f\in \mathcal{F}$ on $[0,1]$ from a sample $W_1,\dots ,W_N\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }f$ by

$$\begin{array}{}\text{(4)}& \check{f}(x)=\sum _{l=1}^N\kappa (\frac{x-W_l}{h})w(x,h)/\sum _{l=1}^N{\int }_0^1\kappa (\frac{y-W_l}{h})w(y,h)dy,\end{array}$$

for $x\in [0,1]$ and $0$ elsewhere.

Here, the kernel $\kappa$ is assumed to satisfy the following additional conditions:

1. (K1) The kernel $\kappa$ is of bounded variation and is symmetric about $0$.

2. (K2) The kernel $\kappa$ satisfies ${\int }_0^1\kappa (u)du>0$, and ${\int }_{\mathbb{R}}|u|\kappa (u)du$, ${\int }_{\mathbb{R}}{\kappa }^2(u)du$ and ${\int }_{\mathbb{R}}|u|{\kappa }^2(u)du$ are finite.

The weight function

$$\begin{array}{}\text{(5)}& w(x,h)=\{\begin{array}{ll}({\int }_{-x/h}^1\kappa (u)du{)}^{-1}& \text{for }x\in [0,h),\\ ({\int }_{-1}^{(1-x)/h}\kappa (u)du{)}^{-1},& \text{for }x\in (1-h,1]\\ 1& \text{otherwise}\end{array}\end{array}$$

is designed to remove boundary bias.

The following result demonstrates that this modified kernel estimator indeed satisfies conditions $\text{(2)}$ and $\text{(3)}$. Furthermore, this result provides the rate in $\text{(2)}$ for this estimator as $b_N=N^{-1/3}$, which is known to be the optimal rate under our assumptions, where the class of densities $\mathcal{F}$ is assumed to be continuously differentiable, and it also shows that rates $a_N=N^{-c}$, for any $c\in (0,1/6)$ are possible in $\text{(2)}$.

Functional data analysis for density process

For a generic density function process $f\sim \mathfrak{F}$, denote:

1. the mean function by $\mu (x)=E(f(x))$,

2. the covariance function by $G(x,y)=\mathrm{Cov}(f(x),f(y))$,

3. the orthonormal eigenfunctions and eigenvalues of the linear covariance operator $$\begin{array}{}\text{(6)}& (Af)(t)=\int G(s,t)f(s)\mathrm{d}s\end{array}$$

by $\{{\varphi }_k{\}}_{k=1}^{\mathrm{\infty }}$ and $\{{\lambda }_k{\}}_{k=1}^{\mathrm{\infty }}$, where the latter are positive and in decreasing order.

If $f_1,\dots ,f_n$ are i.i.d. distributed as $f$, then by the Karhunen–Loève expansion, for each $i$,

$$\begin{array}{}\text{(7)}& f_i(x)=\mu (x)+\sum _{k=1}^{\mathrm{\infty }}{\xi }_{ik}{\varphi }_k(x),\end{array}$$

where ${\xi }_{ik}={\int }_0^1(f_i(x)-\mu (x)){\varphi }_k(x)\mathrm{d}x$ are the uncorrelated principal components with zero mean and variance ${\lambda }_k$. The Karhunen–Loève expansion constitutes the foundation for the commonly used FPCA technique.

The next part of this section is on the estimation of the mode of variation, i.e., the eigenfunctions. This part is the routine part of the FPCA, and I will skip it.

Transformation approach

1. The proposed transformation approach is to map the densities into a new space $L^2(\mathcal{T})$ via a functional transformation $\psi$, where $\mathcal{T}\subset \mathbb{R}$ is a compact interval.

2. Then we work with the resulting $L^2$ process $X:=\psi (f)$

3. By performing FPCA in the linear space $L^2(\mathcal{T})$ and then mapping back to density space, this transformation approach can be viewed as an intrinsic analysis, as opposed to ordinary FPCA.

In the new space $L^2(\mathcal{T})$:

• $\nu$ and $H$ denoting the mean and covariance functions, respectively, of the process $X$.
• $\{{\rho }_k{\}}_{k=1}^{\mathrm{\infty }}$ denoting the orthonormal eigenfunctions of the covariance
• $\{{\eta }_k{\}}_{k=1}^{\mathrm{\infty }}$ is the corresponding eigenvalues,

The Karhunen–Loeve expansion for each of the transformed processes $X_i=\psi (f_i)$ is

$$X_i(t)=\nu (t)+\sum _{k=1}^{\mathrm{\infty }}{\eta }_{ik}{\rho }_k(t),t\in \mathcal{T},$$

with principal components ${\eta }_{ik}={\int }_{\mathcal{T}}(X_i(t)-\nu (t)){\rho }_k(t)dt$.

Our goal is to find suitable transformations $\psi :\mathcal{G}\to L^2(\mathcal{T})$ from density space to a linear functional space.

We begin with two specific examples of relevant transformations. For clarity, for functions in the native density space $\mathcal{G}$ we denote the argument by $x$, while for functions in the transformed space $L^2(\mathcal{T})$ the argument is $t$.

The log hazard transformation

Since hazard functions diverge at the right endpoint of the distribution, which is 1, we consider quotient spaces induced by identifying densities which are equal on a subdomain $\mathcal{T}=[0,1_{\delta }]$, where $1_{\delta }=1-\delta$ for some $0<\delta <1$. With a slight abuse of notation, we denote this quotient space as $\mathcal{G}$ as well.

The log hazard transformation ${\psi }_H:\mathcal{G}\to L^2(\mathcal{T})$ is $${\psi }_H(f)(t)=\log (h(t))=\log \{\frac{f(t)}{1-F(t)}\},t\in \mathcal{T}.$$

Since the hazard function is positive but otherwise not constrained on $\mathcal{T}$, it is easy to see that $\psi$ indeed maps density functions to $L^2(\mathcal{T})$. The inverse map can be defined for any continuous function $X$ as

$${\psi }_H^{-1}(X)(x)=\exp \{X(x)-{\int }_0^x{e}^{X(s)}ds\},x\in [0,1_{\delta }].$$

Solving the differential equation $$h(t)=\frac{f(t)}{1-F(t)}=\frac{F^{\prime }(t)}{1-F(t)}$$ for $F(t)$, it can be shown that $$F(t)=1-\exp (-{\int }_0^{t}h(t)dt)$$

Note that for this case one has a strict inverse only modulo the quotient space. However, in order to use metrics such as $d_W$, we must choose a representative. A~straightforward way to do this is to assign the remaining mass uniformly, that is,

$${\psi }_H^{-1}(X)(x)={\delta }^{-1}\exp \{-{\int }_0^{1_{\delta }}{e}^{X(s)}ds\},x\in (1_{\delta },1].$$

The log quantile density transformation

For $\mathcal{T}=[0,1]$, the log quantile density (LQD) transformation ${\psi }_Q:\mathcal{G}\to L^2(\mathcal{T})$ is given by

$${\psi }_Q(f)(t)=\log (q(t))=-\log \{f(Q(t))\},t\in \mathcal{T}.$$

It is then natural to define the inverse of a continuous function $X$ on $\mathcal{T}$ as the density given by $\exp \{-X(F(x))\}$, where $Q(t)=F^{-1}(t)={\int }_0^t{e}^{X(s)}ds$. Since the value $F^{-1}(1)$ is not fixed, the support of the densities is not fixed within the transformed space, and as the inverse transformation should map back into the space of densities with support on $[0,1]$, we make a slight adjustment when defining the inverse by

$${\psi }_Q^{-1}(X)(x)={\theta }_X\exp \{-X(F(x))\},F^{-1}(t)={\theta }_X^{-1}{\int }_0^t{e}^{X(s)}ds,$$

where ${\theta }_X={\int }_0^1{e}^{X(s)}ds$. Since $F^{-1}(1)=1$ whenever $X\in {\psi }_Q(\mathcal{G})$, this definition coincides with the natural definition mentioned above on ${\psi }_Q(\mathcal{G})$.

In the transformation approach we construct modes of variation in the transformed space for processes $X=\psi (f)$ and then map these back into the density space, defining transformation modes of variation

$$\begin{array}{}\text{(8)}& g_k(x,\alpha ,\psi )={\psi }^{-1}(\nu +\alpha \sqrt{{\tau }_k}{\rho }_k)(x).\end{array}$$

Estimation of these modes is done by first estimating the mean function $\nu$ and covariance function $H$ of the process $X$. Letting ${\hat{X}}_i=\psi ({\check{f}}_i)$, the empirical estimators are

$$\begin{array}{rcl}\text{(9)}& \tilde{\nu }(t)& =& \frac{1}{n}\sum _{i=1}^n{X}_i(t)\text{respectively}\hat{\nu }(t)=\frac{1}{n}\sum _{i=1}^n{\hat{X}}_i(t);\tilde{H}(s,t)& =& \frac{1}{n}\sum _{i=1}^n{X}_i(s)X_i(t)-\tilde{\nu }(s)\tilde{\nu }(t)\text{respectively}\\ \text{(10)}& \hat{H}(s,t)& =& \frac{1}{n}\sum _{i=1}^n{\hat{X}}_i(s){\hat{X}}_i(t)-\hat{\nu }(s)\hat{\nu }(t).\end{array}$$

Estimated eigenvalues and eigenfunctions (${\tilde{\tau }}_k$ and ${\tilde{\rho }}_k$, resp., ${\hat{\tau }}_k$ and ${\hat{\rho }}_k$) are then obtained from the mean and covariance estimates as before, yielding the transformation mode of variation estimators $$\begin{array}{rcl}{\tilde{g}}_k(x,\alpha ,\psi )& =& {\psi }^{-1}(\tilde{\nu }+\alpha \sqrt{{\tilde{\tau }}_k}{\tilde{\rho }}_k)(x)\text{respectively}\\ \text{(11)}& {\hat{g}}_k(x,\alpha ,\psi )& =& {\psi }^{-1}(\hat{\nu }+\alpha \sqrt{{\hat{\tau }}_k}{\hat{\rho }}_k)(x).\end{array}$$

In contrast to the modes of variation resulting from ordinary FPCA in Section 3, the transformation modes are bona fide density functions for any value of $\alpha$. Thus, for reasonably chosen transformations, the transformation modes can be expected to provide a more interpretable description of the variability contained in the sample of densities. Indeed, the data application in Section 6.2 shows that this is the case, using the log quantile density transformation as an example.

The truncated representations of the original densities in the sample are then given by

$$\begin{array}{}\text{(12)}& f_i(x,K,\psi )={\psi }^{-1}(\nu +\sum _{k=1}^K{\eta }_{ik}{\rho }_k)(x).\end{array}$$

Utilizing ($\text{9}$) and the ensuing estimates of the eigenfunctions, the (transformation) principal components, for the case of fully observed densities, are obtained in a straightforward manner,

$$\begin{array}{}\text{(13)}& {\tilde{\eta }}_{ik}={\int }_{\mathcal{T}}(X_i(t)-\tilde{\nu }(t)){\tilde{\rho }}_k(t)dt,\end{array}$$

whence $${\tilde{f}}_i(x,K,\psi )={\psi }^{-1}(\tilde{\nu }+\sum _{k=1}^K{\tilde{\eta }}_{ik}{\tilde{\rho }}_k)(x).$$

In practice, the truncation point $K$ can be selected by choosing a cutoff for the fraction of variance explained. This raises the question of how to quantify \emph{total variance}. For the chosen metric $d$, we propose to use the Frechet variance

$$\begin{array}{}\text{(14)}& V_{\mathrm{\infty }}:=E(d(f,f_{\oplus }{)}^2),\end{array}$$ where $f_{\oplus }$ is the Frechet mean of the densities in the sample. which is estimated by its empirical version

$$\begin{array}{}\text{(15)}& {\tilde{V}}_{\mathrm{\infty }}=\frac{1}{n}\sum _{i=1}^{n}d(f_i,{\tilde{f}}_{\oplus }{)}^2,\end{array}$$

using an estimator ${\tilde{f}}_{\oplus }$ of the Frechet mean. Truncating at $K$ included components as in ($\text{12}$) and denoting the truncated versions as $f_{i,K}$, the variance explained by the first $K$ components is

$$\begin{array}{}\text{(16)}& V_K:=V_{\mathrm{\infty }}-E(d(f_1,f_{1,K}{)}^2),\end{array}$$

which is estimated by

$$\begin{array}{}\text{(17)}& {\tilde{V}}_K={\tilde{V}}_{\mathrm{\infty }}-\frac{1}{n}\sum _{i=1}^{n}d(f_i,{\tilde{f}}_{i,K}{)}^2.\end{array}$$

The ratio $V_K/V_{\mathrm{\infty }}$ is called the fraction of variance explained (FVE), and is estimated by ${\tilde{V}}_K/{\tilde{V}}_{\mathrm{\infty }}$. If the truncation level is chosen so that a fraction $p$, $0<p<1$, of total variation is to be explained, the optimal choice of $K$ is

$$\begin{array}{}\text{(18)}& K^{\ast }=min\{K:\frac{V_K}{V_{\mathrm{\infty }}}>p\},\end{array}$$

which is estimated by

$$\begin{array}{}\text{(19)}& {\tilde{K}}^{\ast }=min\{K:\frac{{\tilde{V}}_K}{{\tilde{V}}_{\mathrm{\infty }}}>p\}.\end{array}$$

As will be demonstrated in the data illustrations, this more general notion of variance explained is a useful concept when dealing with densities or other functions that are not in a Hilbert space. Specifically, we will show that density representations in $\text{(12)}$, obtained via transformation, yield higher FVE values than the ordinary representations, thus giving more efficient representations of the sample of densities.

Reference