Compositional splines for representation of density functions

TL;DR

This paper introduces compositional splines, a new class of spline functions designed to accurately approximate probability density functions while respecting their scale and integral properties, enhancing functional data analysis.

Contribution

The paper develops compositional splines based on the Bayes space approach, enabling consistent density approximation and smoothing with potential orthonormalization for statistical analysis.

Findings

Effective density approximation demonstrated on anthropometric data.

Construction of smoothing compositional splines shown to be practical.

Potential for improved functional principal component analysis of densities.

Abstract

In the context of functional data analysis, probability density functions as non-negative functions are characterized by specific properties of scale invariance and relative scale which enable to represent them with the unit integral constraint without loss of information. On the other hand, all these properties are a challenge when the densities need to be approximated with spline functions, including construction of the respective spline basis. The Bayes space methodology of density functions enables to express them as real functions in the standard $L^2$ space using the centered log-ratio transformation. The resulting functions satisfy the zero integral constraint. This is a key to propose a new spline basis, holding the same property, and consequently to build a new class of spline functions, called compositional splines, which can approximate probability density functions in a consistent way. The paper provides also construction of smoothing compositional splines and possible orthonormalization of the spline basis which might be useful in some applications. Finally, statistical processing of densities using the new approximation tool is demonstrated in case of simplicial functional principal component analysis with anthropometric data.

Figures

Tables

Table 1. Table 1 : Histogram data for four age groups: [ 15 , 16 ) , [ 22 , 23 ) , [ 23 , 24 ) 15 16 22 23 23 24 \left[15,16\right),\left[22,23\right),\left[23,24\right) and [ 30 , 31 ) 30 31 \left[30,31\right) . f i subscript f 𝑖 \textbf{f}_{i} are raw density values at midpoints t i = ( t i , 1 , … , t i , q i ) ⊤ subscript t 𝑖 superscript subscript 𝑡 𝑖 1 … subscript 𝑡 𝑖 subscript 𝑞 𝑖 top \textbf{t}_{i}=\left(t_{i,1},\ldots,t_{i,q_{i}}\right)^{\top} of weight classes with proportions p i subscript p 𝑖 \textbf{p}_{i} for i = 1 , 8 , 9 , 16 𝑖 1 8 9 16 i=1,8,9,16 ; q i subscript 𝑞 𝑖 q_{i} indicates the number of the weight classes.

$[15,16)$	${\text{p}}_1$	0.0656	0.2625	0.3375	0.2156	0.0750	0.0281	0.0094	0.0062			$q_1=8$
	${\text{f}}_1$	0.0075	0.0300	0.0386	0.0246	0.0086	0.0032	0.0011	0.0007
	${\text{t}}_1$	44.375	53.125	61.875	70.625	79.375	88.125	96.875	105.625
$[22,23)$	${\text{p}}_8$	0.0156	0.0869	0.1804	0.2138	0.2272	0.1514	0.0935	0.0200	0.0067	0.0045	$q_8=10$
	${\text{f}}_8$	0.0022	0.0124	0.0258	0.0305	0.0325	0.0216	0.0134	0.0029	0.0010	0.0006
	${\text{t}}_8$	43.5	50.5	57.5	64.5	71.5	78.5	85.5	92.5	99.5	106.5
$[23,24)$	${\text{p}}_9$	0.0078	0.0908	0.2100	0.1971	0.1659	0.1659	0.1011	0.0259	0.0337	0.0017	$q_9=10$
	${\text{f}}_9$	0.0011	0.0130	0.0300	0.0282	0.0237	0.0237	0.0144	0.0037	0.0048	0.0002
	${\text{t}}_9$	43.5	50.5	57.5	64.5	71.5	78.5	85.5	92.5	99.5	106.5
$[30,31)$	${\text{p}}_{16}$	0.0568	0.1023	0.2045	0.2386	0.2159	0.1364	0.0455				$q_{16}=7$
	${\text{f}}_{16}$	0.0057	0.0102	0.0205	0.0239	0.0216	0.0136	0.0045
	${\text{t}}_{16}$	45.0	55.0	65.0	75.0	85.0	95.0	105.0

Table 2. Table 2 : Input data for smoothing procedure: clr ​ ( f i ) clr subscript f 𝑖 \mbox{clr}(\textbf{f}_{i}) are raw clr density values at midpoints t i = ( t i , 1 , … , t i , q i ) ⊤ subscript t 𝑖 superscript subscript 𝑡 𝑖 1 … subscript 𝑡 𝑖 subscript 𝑞 𝑖 top \textbf{t}_{i}=\left(t_{i,1},\ldots,t_{i,q_{i}}\right)^{\top} of the weight classes for i = 1 , 2 , … , N 𝑖 1 2 … 𝑁 i=1,2,\ldots,N ; q i subscript 𝑞 𝑖 q_{i} indicates the number of the weight classes.

$[15,16)$	$\text{clr}\left[{\text{f}}_1\right]$	0.100	1.486	1.737	1.289	0.233	-0.748	-1.846	-2.252			$q_1=8$
	${\text{t}}_1$	44.375	53.125	61.875	70.625	79.375	88.125	96.875	105.625
$[16,17)$	$\text{clr}\left[{\text{f}}_2\right]$	-0.210	1.217	1.760	1.636	0.396	-0.392	-2.001	-2.407			$q_2=8$
	${\text{t}}_2$	44.375	53.125	61.875	70.625	79.375	88.125	96.875	105.625
$[17,18)$	$\text{clr}\left[{\text{f}}_3\right]$	-1.375	0.570	1.316	1.669	1.381	0.534	-0.364	-2.069	-1.663		$q_3=9$
	${\text{t}}_3$	43.889	51.667	59.444	67.222	75.000	82.778	90.556	98.333	106.111
$[18,19)$	$\text{clr}\left[{\text{f}}_4\right]$	-1.354	0.592	1.419	1.443	1.406	1.131	0.563	-0.661	-1.171	-3.369	$q_4=10$
	${\text{t}}_4$	43.5	50.5	57.5	64.5	71.5	78.5	85.5	92.5	99.5	106.5
$[19,20)$	$\text{clr}\left[{\text{f}}_5\right]$	-1.536	0.628	1.408	1.555	1.535	1.209	0.302	-0.774	-1.536	-2.789	$q_5=10$
	${\text{t}}_5$	43.5	50.5	57.5	64.5	71.5	78.5	85.5	92.5	99.5	106.5
$[20,21)$	$\text{clr}\left[{\text{f}}_6\right]$	-1.341	0.674	1.333	1.558	1.638	1.452	0.422	-0.568	-2.034	-3.133	$q_6=10$
	${\text{t}}_6$	43.5	50.5	57.5	64.5	71.5	78.5	85.5	92.5	99.5	106.5
$[21,22)$	$\text{clr}\left[{\text{f}}_7\right]$	-1.746	0.451	1.185	1.463	1.411	1.131	0.531	-0.242	-1.746	-2.439	$q_7=10$
	${\text{t}}_7$	43.5	50.5	57.5	64.5	71.5	78.5	85.5	92.5	99.5	106.5
$[22,23)$	$\text{clr}\left[{\text{f}}_8\right]$	-1.168	0.550	1.281	1.450	1.511	1.106	0.624	-0.917	-2.015	-2.421	$q_8=10$
	${\text{t}}_8$	43.5	50.5	57.5	64.5	71.5	78.5	85.5	92.5	99.5	106.5
$[23,24)$	$\text{clr}\left[{\text{f}}_9\right]$	-1.884	0.573	1.412	1.348	1.177	1.177	0.681	-0.680	-0.417	-3.388	$q_9=10$
	${\text{t}}_9$	43.5	50.5	57.5	64.5	71.5	78.5	85.5	92.5	99.5	106.5
$[24,25)$	$\text{clr}\left[{\text{f}}_{10}\right]$	-1.602	0.595	1.186	1.274	1.106	0.796	0.056	-0.423	-2.988		$q_{10}=9$
	${\text{t}}_{10}$	43.889	51.667	59.444	67.222	75.000	82.778	90.556	98.333	106.111
$[25,26)$	$\text{clr}\left[{\text{f}}_{11}\right]$	-1.401	0.471	0.768	0.824	1.145	0.850	0.209	-1.178	-1.688		$q_{11}=9$
	${\text{t}}_{11}$	43.889	51.667	59.444	67.222	75.000	82.778	90.556	98.333	106.111
$[26,27)$	$\text{clr}\left[{\text{f}}_{12}\right]$	-1.045	0.513	0.901	1.180	1.258	0.513	-0.485	-2.836			$q_{12}=8$
	${\text{t}}_{12}$	44.375	53.125	61.875	70.625	79.375	88.125	96.875	105.625
$[27,28)$	$\text{clr}\left[{\text{f}}_{13}\right]$	-0.816	0.570	0.742	0.742	1.056	0.570	-0.256	-2.608			$q_{13}=8$
	${\text{t}}_{13}$	44.375	53.125	61.875	70.625	79.375	88.125	96.875	105.625
$[28,29)$	$\text{clr}\left[{\text{f}}_{14}\right]$	-1.155	0.579	0.790	0.965	0.690	-0.308	-1.561				$q_{14}=7$
	${\text{t}}_{14}$	45.0	55.0	65.0	75.0	85.0	95.0	105.0
$[29,30)$	$\text{clr}\left[{\text{f}}_{15}\right]$	-1.060	0.480	0.837	0.674	0.614	-0.773	-0.773				$q_{15}=7$
	${\text{t}}_{15}$	45.0	55.0	65.0	75.0	85.0	95.0	105.0
$[30,31)$	$\text{clr}\left[{\text{f}}_{16}\right]$	–0.756	-0.168	0.525	0.679	0.579	0.120	-0.979				$q_{16}=7$
	${\text{t}}_{16}$	45.0	55.0	65.0	75.0	85.0	95.0	105.0

Table 3. Table 3 : Z ​ B 𝑍 𝐵 ZB -spline coefficients for clr transformed density functions of N = 16 𝑁 16 N=16 age groups.

age group	spline coefficients, ${\text{z}}_i={(z_{i,-3},\dots ,z_{i,1})}^{\top },i=1,\dots ,N$
$[15,16)$	-6.950	6.647	46.536	40.973	13.163
$[16,17)$	-7.806	-0.596	41.616	45.181	14.083
$[17,18)$	-16.677	-11.292	18.284	43.917	9.102
$[18,19)$	-17.067	-8.988	21.373	33.533	20.188
$[19,20)$	-18.483	-9.902	22.408	38.249	16.447
$[20,21)$	-17.242	-7.010	18.199	46.788	18.682
$[21,22)$	-20.452	-10.875	11.653	36.887	14.797
$[22,23)$	-15.236	-5.368	16.735	46.421	14.071
$[23,24)$	-22.485	-12.348	17.033	23.450	20.153
$[24,25)$	-19.873	-14.176	13.567	20.115	19.448
$[25,26)$	-19.011	-5.949	-4.623	30.860	9.973
$[26,27)$	-14.997	-10.545	2.638	28.225	19.143
$[27,28)$	-14.461	-4.455	-0.689	21.892	18.070
$[28,29)$	-18.518	-11.045	-2.723	21.744	10.395
$[29,30)$	-16.445	-9.417	-1.814	23.562	2.889
$[30,31)$	-5.077	-15.534	-4.171	8.220	7.618

Table 4. Table 4 : B 𝐵 B -spline coefficients for clr transformed density functions of N = 16 𝑁 16 N=16 age groups.

age group	spline coefficients, ${\text{b}}_i={(b_{i,-3},\dots ,b_{i,2})}^{\top },i=1,\dots ,N$
$[15,16)$	-1.264	1.236	2.381	-0.332	-2.472	-2.289
$[16,17)$	-1.419	0.655	2.520	0.213	-2.764	-2.449
$[17,18)$	-3.032	0.490	1.766	1.530	-3.095	-1.583
$[18,19)$	-3.103	0.734	1.813	0.726	-1.186	-3.511
$[19,20)$	-3.361	0.780	1.929	0.946	-1.938	-2.860
$[20,21)$	-3.135	0.930	1.505	1.707	-2.498	-3.249
$[21,22)$	-3.719	0.871	1.345	1.507	-1.964	-2.573
$[22,23)$	-2.770	0.897	1.320	1.772	-2.876	-2.447
$[23,24)$	-4.088	0.922	1.754	0.383	-0.293	-3.505
$[24,25)$	-3.613	0.518	1.656	0.391	-0.059	-3.382
$[25,26)$	-3.456	1.187	0.079	2.118	-1.857	-1.734
$[26,27)$	-2.727	0.405	0.787	1.528	-0.807	-3.329
$[27,28)$	-2.629	0.910	0.225	1.348	-0.340	-3.143
$[28,29)$	-3.367	0.679	0.497	1.461	-1.009	-1.808
$[29,30)$	-2.990	0.639	0.454	1.515	-1.838	-0.502
$[30,31)$	-0.923	-0.951	0.678	0.740	-0.053	-1.325

Equations

$$(f\oplus g)(x)=\frac{f(x)g(x)}{{\int }_{I}f(y)g(y)\mathrm{d}y},(\alpha \odot f)(x)=\frac{f(x{)}^{\alpha }}{{\int }_{I}f(y{)}^{\alpha }\mathrm{d}y},x\in I,$$

$$⟨f,g{⟩}_{\mathcal{B}}=\frac{1}{2\eta }{\int }_I{\int }_I\ln \frac{f(x)}{f(y)}\ln \frac{g(x)}{g(y)}\mathrm{d}x\mathrm{d}y,f,g\in {\mathcal{B}}^2(I).$$

$$\text{mbox}clr(f)(x)\equiv f_c(x)=\ln f(x)-\frac{1}{\eta }{\int }_I\ln f(y)\mathrm{d}y.$$

$$\text{mbox}clr(f\oplus g)(x)=f_c(x)+g_c(x),\text{mbox}clr(\alpha \odot f)(x)=\alpha \cdot f_c(x)$$

$$⟨f,g{⟩}_{\mathcal{B}}=⟨f_c,g_c{⟩}_2={\int }_I{f}_c(x)g_c(x)\mathrm{d}x.$$

$${\int }_I\text{mbox}clr(f)(x)\mathrm{d}x={\int }_I\ln f(x)\mathrm{d}x-{\int }_I\frac{1}{\eta }{\int }_I\ln f(y)\mathrm{d}y\mathrm{d}x=0,$$

$$\text{mbox}{clr}^{-1}[f_c](x)=\frac{\exp (f_c(x))}{{\int }_I\exp (f_c(y))\text{mbox}dy};$$

$$f(x)=c_1\odot {\phi }_1(x)\oplus c_2\odot {\phi }_2(x)\oplus \dots \oplus c_r\odot {\phi }_r(x)\oplus \dots =\underset{i=1}{\overset{\infty }{⨁}}{c}_i\odot {\phi }_i(x),x\in I.$$

$${\lambda }_0=a<{\lambda }_1<\cdots <{\lambda }_g<b={\lambda }_{g+1}$$

$$B_{i}^{1}(x)\,=\,\left\{\begin{array}[]{rl}1&\quad\mbox{if}\;x\in[\lambda_{i},\lambda_{i+1})\\ 0&\quad\mbox{otherwise}\end{array}\right.$$

$$B_i^{k+1}(x)=\frac{x-{\lambda }_i}{{\lambda }_{i+k}-{\lambda }_i}{B}_i^k(x)+\frac{{\lambda }_{i+k+1}-x}{{\lambda }_{i+k+1}-{\lambda }_{i+1}}{B}_{i+1}^k(x).$$

$$Z_i^{k+1}(x):=\frac{\text{mbox}d}{\text{mbox}dx}{B}_i^{k+2}(x),$$

$$Z_{i}^{1}(x)=\left\{\begin{array}[]{rl}\dfrac{1}{\lambda_{i+1}-\lambda_{i}}&\quad\mbox{if}\;x\in[\lambda_{i},\lambda_{i+1})\\ &\\ \dfrac{-1}{\lambda_{i+2}-\lambda_{i+1}}&\quad\mbox{if}\;x\in(\lambda_{i+1},\lambda_{i+2}]\end{array}\right.$$

$$Z_i^{k+1}(x)=(k+1)\left(\frac{B_i^{k+1}(x)}{{\lambda }_{i+k+1}-{\lambda }_i}-\frac{B_{i+1}^{k+1}(x)}{{\lambda }_{i+k+2}-{\lambda }_{i+1}}\right).$$

$$\text{mbox}supp{Z}_i^{k+1}(x)=\text{mbox}supp{B}_i^{k+2}(x)=[{\lambda }_i,{\lambda }_{i+k+2}],$$

$$Z_i^{k+1}(x)=0\text{mbox}ifx\notin [{\lambda }_i,{\lambda }_{i+k+2}].$$

$$M_i^{k+1}(x):=\frac{k+1}{{\lambda }_{i+k+1}-{\lambda }_i}{B}_i^{k+1}(x)$$

$$\underset{\mathbb{R}}{\int }{M}_i^{k+1}(x)\text{mbox}dx=1,$$

$$Z_i^{k+1}(x)=M_i^{k+1}(x)-M_{i+1}^{k+1}(x)$$

$$\underset{\mathbb{R}}{\int }{Z}_i^{k+1}(x)\text{mbox}dx=0.$$

$$\dim ({\mathcal{S}}_k^{\Delta \lambda }[a,b])=g+k+1.$$

$${\lambda }_{-k}=\cdots ={\lambda }_{-1}={\lambda }_0=a,b={\lambda }_{g+1}={\lambda }_{g+2}=\cdots ={\lambda }_{g+k+1}.$$

$$s_k\left(x\right)=\sum _{i=-k}^g{b}_i{B}_i^{k+1}\left(x\right).$$

$${\mathcal{Z}}_k^{\Delta \lambda }[a,b]:=\{s_k(x)\in {\mathcal{S}}_k^{\Delta \lambda }[a,b]:\underset{I}{\int }{s}_k(x)\text{mbox}dx=0\}.$$

$$\underset{I}{\int }{s}_k(x)\text{mbox}dx=\frac{1}{k+1}\sum _{i=-k}^g{b}_i({\lambda }_{i+k+1}-{\lambda }_i).$$

$$\dim ({\mathcal{Z}}_k^{\Delta \lambda }[a,b])=\dim ({\mathcal{S}}_k^{\Delta \lambda }[a,b])-\text{mbox}codim({\mathcal{Z}}_k^{\Delta \lambda }[a,b])=g+k.$$

$$s_k\left(x\right)=\sum _{i=-k}^{g-1}{z}_i{Z}_i^{k+1}\left(x\right).$$

$$Z_{i}^{k+1}\left(x\right)=(k+1)\left(B_{i}^{k+1}\left(x\right),B_{i+1}^{k+1}\left(x\right)\right)\left(\begin{array}[]{cc}\dfrac{1}{\lambda_{i+k+1}-\lambda_{i}}&0\\ 0&\dfrac{1}{\lambda_{i+k+2}-\lambda_{i+1}}\end{array}\right)\left(\begin{array}[]{r}1\\ -1\end{array}\right).$$

$$(Z_{-k}^{k+1}\left(x\right),\dots ,Z_{g-1}^{k+1}\left(x\right))=(B_{-k}^{k+1}\left(x\right),\dots ,B_g^{k+1}\left(x\right))\mathbf{DK}={\mathbf{B}}_{k+1}(x)\mathbf{DK},$$

$$\mathbf{D}=(k+1)\text{mbox}diag\left(\frac{1}{{\lambda }_1-{\lambda }_{-k}},\dots ,\frac{1}{{\lambda }_{g+k+1}-{\lambda }_g}\right)=(k+1)\text{mbox}diag\left(\frac{1}{l_1},\dots ,\frac{1}{l_{g+k+1}}\right)$$

Full text

Keywords:

spline representation, spline with zero integral, compositional spline, smoothing spline, simplicial functional principal component analysis

1 Introduction

Probability density functions are popularly known as non-negative functions satisfying the unit integral constraint. This clearly inhibits their direct processing using standard methods of functional data analysis [18] since unconstrained functions are assumed there. The same holds also for approximation of the raw input data using splines which is commonly considered to be a key step in functional data analysis. But more severely, in addition to the apparent unit integral constraint of densities which might seem to represent just a kind of numerical obstruction, density functions are rather characterized by deeper geometrical properties that need to be taken into account for any reliable analysis [6, 21, 22]. Specifically, in contrast to functions in the standard $L^{2}$ space, densities obey the scale invariance and relative scale properties [9]. Scale invariance means that not just the representation of densities with the unit integral constraint, but any its positive multiple conveys the same information about relative contributions of Borel sets on the whole probability mass. Relative scale can be explained directly with an example: the relative increase of a probability over a Borel set from 0.05 to 0.1 (2 multiple) differs from the increase 0.5 to 0.55 (1.1 multiple), although the absolute differences are the same in both cases. If we restrict to a bounded support $I=[a,b]\subset\mathbb{R}$ that is mostly used in practical applications [3, 9, 14, 15], density functions can be represented with respect to Lebesgue reference measure using the Bayes space $\mathcal{B}^{2}(I)$ of functions with square-integrable logarithm [6, 22].

The Bayes space $\mathcal{B}^{2}(I)$ has structure of separable Hilbert space that enables construction of an isometric isomorphism between $\mathcal{B}^{2}(I)$ and $L^{2}(I)$, the $L^{2}$ space restricted to $I$. Accordingly, analogies of summing two functions and multiplication of a function by a real scalar in the $L^{2}$ space together with an inner product between two densities are required. Given two absolutely integrable density functions $f,g\in\mathcal{B}^{2}(I)$ and a real number $\alpha\in\mathbb{R}$ we indicate with $f\oplus g$ and $\alpha\odot f$ the perturbation and powering operations, defined as

[TABLE]

respectively. The resulting functions are readily seen to be probability density functions, though, note that the unit integral constraint representation was chosen just for the sake of convenient interpretation. In [6], it is proven that $\mathcal{B}^{2}(I)$ endowed with the operations $(\oplus,\odot)$ is a vector space. Note that the neutral elements of perturbation and powering are $e(x)=1/\eta$, with $\eta=b-a$ (i.e., the uniform density), and 1, respectively. The difference between two elements $f,g\in\mathcal{B}^{2}(I)$, denoted by $f\ominus g$, is obtained as perturbation of $f$ with the reciprocal of $g$, i.e., $(f\ominus g)(x)=(f\oplus[(-1)\odot g])(x),x\in I$. Finally, to complete the Hilbert space structure, the inner product is defined as

[TABLE]

Form of the inner product clearly indicates that the relevant information in densities is contained in (log-)ratios between elements from the support $I$.

Density functions can be considered as functional counterparts to compositional data, positive vectors carrying relative information [1, 17] that are driven by the Aitchison geometry [16]. In order to enable their statistical processing using standard multivariate methods in real space [5], the preferred strategy is to express them either in centered log-ratio (clr) coefficients [1] with respect to a generating system, or in logratio coordinates, preferably with respect to an orthonormal basis [7]. The latter coordinates (called also isometric log-ratio coordinates), as well as the clr coefficients, provide isometry between the Aitchison geometry and the real Euclidean space. A similar strategy is used also for densities in the Bayes space [22]. An isometric isomorphism between $\mathcal{B}^{2}(I)$ and $L^{2}(I)$ is represented by the centered log-ratio (clr) transformation [14, 22], defined for $f\in\mathcal{B}^{2}(I)$ as

[TABLE]

We remark that such an isometry allows to compute operations and inner products among the elements in $\mathcal{B}^{2}(I)$ in terms of their counterpart in $L^{2}(I)$ among the clr-transforms, i.e.

[TABLE]

and

[TABLE]

However, the clr transformation induces an additional constraint,

[TABLE]

that needs to be taken into account for computation and analysis on clr-transformed density functions. As the clr space is clearly a subspace of $L^{2}(I)$, hereafter it is denoted as $L_{0}^{2}(I)$. The inverse clr transformation is obtained as

[TABLE]

again as before, the denominator is used just to achieve the unit integral constraint representation of the resulting density (without loss of relative information, carried by the density function).

According to [22], it is not necessary to restrict ourselves to the constrained clr space, because a basis in $\mathcal{B}^{2}(I)$ can be easily constructed. Specifically, let $\psi_{0}(x),\psi_{1}(x),\psi_{2}(x),\dots$ is a basis in $L^{2}(I)$ and assume that $\psi_{0}(x)$ is a constant function, then $\varphi_{1}(x):=\exp(\psi_{1}(x)),$ $\varphi_{2}(x):=\exp(\psi_{2}(x)),\dots$ form a basis in $\mathcal{B}^{2}(I)$. Of course, also here an orthonormal basis is preferable, but it is not always possible in applications. Nevertheless, if this would be so, then a function $f\in\mathcal{B}^{2}(I)$ can be projected orthogonally to the space spanned, e.g., by the first $r$ functions $\varphi_{1}(x),\varphi_{2}(x),\dots,\varphi_{r}(x)$. This is done through the respective coefficients $c_{1},\dots,c_{r}$ in the basis expansion

[TABLE]

Functional data analysis relies strongly on approximation of the input functions using splines [18]. However, splines are mostly utilized purely as an approximation tool, without considering further methodological consequences. Because statistical processing of density functions requires a deeper geometrical background, provided by the Bayes spaces, this should be followed also by the respective spline representation, performed preferably in the clr space $L_{0}^{2}(I)$. In [12], a first attempt of constructing a spline representation that would honor the zero integral constraint (4) was performed. The problem is that $B$-splines that form basis for the spline expansion in [12] come from $L^{2}(I)$, but not from $L_{0}^{2}(I)$. This paper presents an important step ahead – such splines are constructed that form basis functions in the clr space $L_{0}^{2}(I)$. Consequently, the splines can be expressed also directly in $\mathcal{B}^{2}(I)$ and the spline representation formulated in terms of the Bayes space which can be used for interpretation purposes; hereafter we refer to compositional splines. Apart from methodological advantages, using compositional splines simplifies construction and interpretation of spline coefficients that can be considered as coefficients of a (possibly orthonormal) basis in $\mathcal{B}^{2}(I)$.

The paper is organized as follows. In the next section the construction of splines basis in $L_{0}^{2}(I)$ is presented together with a comparison to spline functions introduced in [12]. Section 3 is devoted to smoothing splines in $L_{0}^{2}(I)$ and Section 4 discusses orthogonalization of basis functions (that form, by construction, an oblique basis). Section 5 introduces a new class of splines that reflect the Bayes spaces methodology, compositional splines. Section 6 demonstrates usefulness of the new approximation tool in context of simplicial functional principal component analysis with anthropometric data and the final Section 7 concludes.

2 Construction of spline in $L_{0}^{2}(I)$

Because the clr transformation enables to process density functions in the standard $L^{2}$ space, just restricted according to zero integral constraint (4), it is natural that also construction of compositional splines should start in $L_{0}^{2}(I)$. Nevertheless, before doing so, some basic facts about $B$-spline representation of splines are recalled, see [2, 4, Schum07] for details. Let the sequence of knots

[TABLE]

be given. The (normalized) $B$-spline of degree [math] (order $1$) is defined as

[TABLE]

and the (normalized) $B$-spline of degree $k$, $k\in\mathbb{N}$, (order $k+1$) is defined by

[TABLE]

Now let the functions $Z_{i}^{k+1}(x)$ for $k\geq 0$, $k\in\mathbb{N}$, be defined

[TABLE]

i.e., more precisely for $k=0$

[TABLE]

and for $k\geq 1$

[TABLE]

Noteworthy, functions $Z_{i}^{k+1}(x)$ have similar properties as $B$-splines $B_{i}^{k+1}(x)$.

They are piecewise polynomials of degree $k$. Particularly, $Z_{i}^{1}(x)$ is a piecewise constant polynomial, $Z_{i}^{2}(x)$ is a piecewise linear polynomial, see Figure 1, $Z_{i}^{3}(x)$ is a piecewise quadratic polynomial, see Figure 2. For other examples see Figures 3-5. 2. 2.

It is evident that for $k\geq 1$ the function $Z_{i}^{k+1}(x)$ and its derivatives up to order $k-1$ are all continuous. 3. 3.

It is easy to check that for $k\geq 0$

[TABLE]

and of course

[TABLE] 4. 4.

From the perspective of $L_{0}^{2}(I)$, a crucial point is that the integral of $Z_{i}^{k+1}(x)$ equals to zero. If we consider Curry-Schoenberg $B$-spline $M_{i}^{k+1}(x)$ [2], which are defined as

[TABLE]

with property

[TABLE]

than it is clear that

[TABLE]

and

[TABLE]

It is known that for the vector space ${\cal S}_{k}^{\Delta\lambda}[a,b]$ of polynomial splines of degree $k>0$, $k\in\mathbb{N}$, defined on a finite interval $I=[a,b]$ with the sequence of knots $\Delta\lambda=\left\{\lambda_{i}\right\}_{i=0}^{g+1}$, $\lambda_{0}=a<\lambda_{1}<\ldots<\lambda_{g}<b=\lambda_{g+1}$, the dimension is

[TABLE]

For the construction of all basis functions $B_{i}^{k+1}(x)$, it is necessary to consider some additional knots. Without loss of generality we can add coincident knots

[TABLE]

Then every spline $s_{k}(x)\in{\cal S}_{k}^{\Delta\lambda}[a,b]$ in $L^{2}(I)$ has a unique representation

[TABLE]

In [12, 20], the splines with zero integral are studied. There is given the necessary and sufficient condition for $B$-splines coefficients of these splines. However, typical $B$-splines $B_{i}^{k+1}(x)$, thus ignoring the constraint (4) in $L_{0}^{2}(I)$ for construction of the $B$-spline basis, were used there.

Now, regarding the definition (7), we are able to use spline functions $Z_{i}^{k+1}(x)$ which have zero integral on $I$ (denoted also as $ZB$-splines in the sequel). In the following, ${\cal Z}_{k}^{\Delta\lambda}[a,b]$ denotes the vector space of polynomial splines of degree $k>0$, defined on a finite interval $[a,b]$ with the sequence of knots $\Delta\lambda$ and having zero integral on $[a,b]$, it means

[TABLE]

Theorem 2.1**.**

The dimension of the vector space ${\cal Z}_{k}^{\Delta\lambda}[a,b]$ defined by the formula (12) is $g+k$.

Proof.

For spline $s_{k}(x)\in{\cal S}_{k}^{\Delta\lambda}[a,b]$, $s_{k}\left(x\right)=\sum\limits_{i=-k}^{g}b_{i}B_{i}^{k+1}\left(x\right)$, with the coincident additional knots it is known, [4], that

[TABLE]

It means that $B$-spline coefficients of $s_{k}(x)\in{\cal Z}_{k}^{\Delta\lambda}[a,b]\subset{\cal S}_{k}^{\Delta\lambda}[a,b]$ satisfy condition $0=\sum\limits_{i=-k}^{g}\,b_{i}(\lambda_{i+k+1}-\lambda_{i})\,=\,\mathbf{A}\mathbf{b}$ with $\mathbf{A}=(\lambda_{1}-\lambda_{-k},\cdots,\lambda_{g+k+1}-\lambda_{g})$, $\mathbf{b}=(b_{-k},\cdots,b_{g})^{\top}$. And it is obvious that $\mbox{codim}({\cal Z}_{k}^{\Delta\lambda}[a,b])=1$, thus

[TABLE]

∎

Theorem 2.2**.**

For the coincident additional knots (10), the functions $Z_{-k}^{k+1}(x),$ $\cdots,Z_{g-1}^{k+1}(x)$ form a basis for the space ${\cal Z}_{k}^{\Delta\lambda}[a,b]$.

Proof.

Since $M_{i}^{k+1}(x)$ form a basis for the spline space ${\cal S}_{k}^{\Delta\lambda}[a,b]$ and $Z_{i}^{k+1}(x)\,=\,M_{i}^{k+1}(x)-M_{i+1}^{k+1}(x)$, the functions $Z_{i}^{k+1}(x)$, $i=-k,\ldots,g-1$, are linearly independent and lie in ${\cal Z}_{k}^{\Delta\lambda}[a,b]$ with $\dim({\cal Z}_{k}^{\Delta\lambda}[a,b])=g+k$. Therefore $Z_{i}^{k+1}(x)$, $i=-k,\ldots,g-1$, form a basis for the ${\cal Z}_{k}^{\Delta\lambda}[a,b]$. ∎

With regard to this theorem, every spline $s_{k}(x)\in{\cal Z}_{k}^{\Delta\lambda}[a,b]$ has a unique representation

[TABLE]

Now we can proceed to matrix notation of $s_{k}(x)\in{\cal Z}_{k}^{\Delta\lambda}[a,b]$. With respect to (8) and (9), we are able to write the functions $Z_{i}^{k+1}\left(x\right)$ in matrix notation as

[TABLE]

Then it is clear that

[TABLE]

where

[TABLE]

and

[TABLE]

Therefore the spline $s_{k}(x)\in{\cal Z}_{k}^{\Delta\lambda}[a,b]$, $s_{k}\left(x\right)=\sum\limits_{i=-k}^{g-1}b_{i}Z_{i}^{k+1}\left(x\right)$ can be written in matrix notation as

[TABLE]

where $\mathbf{Z}_{k+1}(x)=(Z_{-k}^{k+1}\left(x\right),\ldots,Z_{g-1}^{k+1}\left(x\right))$ and $\mathbf{z}=\left(z_{-k},\ldots,z_{g-1}\right)^{\top}$.

Remark 1**.**

The formula (16) is very useful, because we can use the standard $B$-spline basis for working with splines honoring the zero integral constraint, which is very convenient from a computational point of view.

Example 2.1**.**

We consider knots $\Delta\lambda=\left\{\lambda_{i}\right\}_{i=0}^{g+1}$, $\lambda_{0}=0=a<2<5<9<14<b=20=\lambda_{5}$. The task is to find a cubic spline with the given sequence of knots and which has zero integral on the interval $[0,20]$. It is evident that $k=3$, $g=4$. We consider the additional knots

[TABLE]

The basis functions of the space ${\cal Z}_{3}^{\Delta\lambda}[0,20]$ are plotted in Figure 7. Every spline $s_{3}(x)\in{\cal Z}_{3}^{\Delta\lambda}[0,20]$ can be written as

[TABLE]

Thus, e.g., for $\mathbf{z}=(z_{-3},\ldots,z_{3})^{\top}=(0.5,-1,2,3,-8,9,1)^{\top}$ the cubic spline $s_{3}(x)$ with zero integral is plotted in Figure 6.

3 Smoothing spline in $L_{0}^{2}(I)$

In [12], the construction of smoothing splines in the space $L_{0}^{2}(I)$ was studied, however using standard $B$-spline basis functions $B_{i}^{k+1}(x)$. Now we are able to construct smoothing splines in this space with new basis functions $Z_{i}^{k+1}(x)$. For this purpose, let data $(x_{i},y_{i})$, $a\leq x_{i}\leq b$, weights $w_{i}>0$, $i=1,\ldots,n$, sequence of knots $\Delta\lambda=\left\{\lambda_{i}\right\}_{i=0}^{g+1}$, $\lambda_{0}=a<\lambda_{1}<\ldots<\lambda_{g}<b=\lambda_{g+1}$, $n\geq g+1$ and a parameter $\alpha\in(0,1)$ be given. For arbitrary $l\in\left\{1,\ldots,k-1\right\}$ our task is to find a spline $s_{k}(x)\in{\cal Z}_{k}^{\Delta\lambda}[a,b]\subset L_{0}^{2}(I)$, which minimizes the functional

[TABLE]

Note that the choice of parameter $\alpha$ and $l$, where $l$ stands for $l$th derivation, affects smoothness of the resulting spline. Let us denote $\mathbf{x}=\left(x_{1},\ldots,x_{n}\right)^{\top}$, $\mathbf{y}=\left(y_{1},\ldots,y_{n}\right)^{\top}$, $\mathbf{w}=\left(w_{1},\ldots,w_{n}\right)^{\top}$ and $\mathbf{W}=diag\left(\mathbf{w}\right)$. Regarding the representation (13) and matrix notation (16), the functional $J_{l}(s_{k})$ can be written as a quadratic function

[TABLE]

see [10, 11, 12] for details. In fact the matrix

[TABLE]

is positive definite, because $B_{i}^{k+1-l}(x)\geq 0$, $i=-k+l,\ldots,g$ are basis functions. Upper triangular matrix ${\bf S}_{l}={\bf D}_{l}{\bf L}_{l}\ldots{\bf D}_{1}{\bf L}_{1}\in\mathbb{R}^{g+k+1-l,g+k+1}$ has full row rank. ${\bf D}_{j}\in\mathbb{R}^{g+k+1-j,g+k+1-j}$ is a diagonal matrix such that

[TABLE]

with

[TABLE]

and

[TABLE]

Finally, ${\bf B}_{k+1}({\bf x})\in\mathbb{R}^{n,g+k+1}$ stands for the collocation matrix, i.e.

[TABLE]

Using the notation $\mathbf{U}:={\bf D}{\bf K}$,

[TABLE]

and

[TABLE]

it is possible to rewrite the quadratic function $J_{l}(\bf z)$ as

[TABLE]

Our task is to find a spline $s_{k}(x)\in{\cal Z}_{k}^{\Delta\lambda}[a,b]$ which minimizes the functional $J_{l}(s_{k})$, in other words, we want to find a minimum of the function (21). It is obvious that this function has just one minimum if and only if the matrix ${\bf G}$ is positive definite (p.d.). From (20) it can be easily seen that

[TABLE]

From Schoenberg-Whitney theorem and its generalization, see [2] and [10], it is known that matrix ${\bf B}_{k+1}(\bf x)$ is of full column rank if and only if there exists $\{u_{-k},\ldots,u_{g}\}\subset\{x_{1},\ldots,x_{n}\}$ with $u_{i}<u_{i+1}$, $i=-k,\ldots,g-1$, such that $\lambda_{i}<u_{i}<\lambda_{i+k+1}$, $i=-k,\ldots,g$. In this case from the necessary and sufficient condition for a unique minimum of quadratic function, i.e.

[TABLE]

we get a system of linear equations ${\bf Gz}={\bf g}$ and then the unique solution of this system is given by

[TABLE]

Consequently, the resulting smoothing spline is obtained by the formula

[TABLE]

in matrix notation using standard $B$-splines $B_{i}^{k+1}(x)$ as

[TABLE]

where the vector $\mathbf{z}^{*}=\left(z_{-k}^{*},\ldots,z_{g-1}^{*}\right)^{{}^{\prime}}$ is obtained as

.

4 Orthogonalization of basis functions

A further step is to orthogonalize the basis

[TABLE]

of the space ${\cal Z}_{k}^{\Delta\lambda}[a,b]$ that is by construction obligue with respect to the $L^{2}$ space metric. For this purpose the idea presented in [19] is used. We search for a linear transformation $\boldsymbol{\Phi}$ such that

[TABLE]

forms an orthogonal set of basis functions of the space ${\cal Z}_{k}^{\Delta\lambda}[a,b]$, i.e.

[TABLE]

Regarding the lemma presented in [19] and notation used here, we can formulate the following statement.

Lemma 4.1**.**

An invertible transformation $\boldsymbol{\Phi}$ orthogonalizes the basis functions $\mathbf{Z}_{k+1}(x)$ if and only if it satisfies the condition that

[TABLE]

where $\boldsymbol{\Sigma}$ represents the positive definite matrix

[TABLE]

With respect to the definition of basis functions $\mathbf{Z}_{k+1}(x)=\left(Z_{-k}^{k+1}\left(x\right),\ldots,\right.$ $\left.Z_{g-1}^{k+1}\left(x\right)\right)^{\top}$ the matrix $\boldsymbol{\Sigma}$ can be expressed as

[TABLE]

where ${\bf M}\,:=\,{\bf M}_{k0}$. The linear transformation $\boldsymbol{\Phi}$ is not unique and can be computed for example by the Cholesky decomposition. The basis functions

[TABLE]

are orthogonal and have a zero integral. The linear and quadratic $ZB$-splines with zero integral and their orthogonalization are plotted in Figures 8 and 9.

To sum up, the spline $s_{k}(x)$ with zero integral can be constructed as a linear combination of orthogonal splines $O_{i}^{k+1}(x)$ having zero integral in a form

[TABLE]

On the other hand, the standard and well-known $B$-splines $B_{i}^{k+1}(x)$ can be used to represent $s_{k}(x)\in{\cal Z}_{k}^{\Delta\lambda}[a,b]$ in matrix form

[TABLE]

This formulation seems to be very useful because it allows us to use existing $B$-spline codes in software R or Matlab, for example, the collocation matrix or computation of integrals in (19).

5 Compositional splines in the Bayes space

Construction of splines directly in $L_{0}^{2}(I)$ has important practical consequences, however, it is crucial also from the theoretical perspective. Expressing splines as functions in $L_{0}^{2}(I)$ enables to back-transform them to the original Bayes space $\mathcal{B}^{2}(I)$ using (5). It results in compositional $CB$-splines, obtained from (8) as

[TABLE]

Accordingly, for instance $ZB$-splines from Figures 1-5 can be now expressed directly in the Bayes space $\mathcal{B}^{2}(I)$ as $CB$-splines, see Figures 10-12. Note that $CB$-splines $\zeta_{i}^{k+1}(x)$ fulfill the unit integral constraint which is, however, not necessary for further considerations. As a consequence, it is immediate to define vector space ${\cal C}_{k}^{\Delta\lambda}[a,b]$ of compositional polynomial splines of degree $k>0$, defined on a finite interval $[a,b]$ with the sequence of knots $\Delta\lambda$. From isomorphism between ${\cal C}_{k}^{\Delta\lambda}[a,b]$ and ${\cal Z}_{k}^{\Delta\lambda}[a,b]$ it holds that

[TABLE]

Moreover, from isometric properties of clr transformation it follows that every compositional spline $\xi_{k}(x)\in{\cal C}_{k}^{\Delta\lambda}[a,b]$ in $\mathcal{B}^{2}(I)$ can be uniquely represented as

[TABLE]

$CB$-splines $\zeta_{i}^{k+1}(x)$ forming the basis are by the default setting (8) not orthogonal. Their orthogonalization can be done as described in Section 4, i.e. by employing $L_{0}^{2}(I)$ with the back-transformation to $\mathcal{B}^{2}(I)$.

The resulting compositional splines (with either orthogonal, or non-orthogonal $CB$-spline basis) can be used for representation of densities directly in $\mathcal{B}^{2}(I)$. This is an important step in construction of methods of functional data analysis involving density functions, like for ANOVA modeling [22] or for the SFPCA method introduced in [9] and demonstrated in the next section. In the latter case $CB$-splines were indeed already used for construction of the procedure, although at the respective level of development the authors were not aware of that. With $CB$-splines one has a guarantee that methods are developed consistently in the Bayes space. Moreover, the possibility of having an orthogonal basis enables to gain additional features resulting from orthogonality of finite dimensional projection in combination with approximate properties of spline functions.

As usual, compositional splines can be tuned according to concrete problem, with the advantage of their direct formulation in the Bayes space sense.

Example 5.1**.**

To illustrate smoothing of concrete data with a compositional spline, 1000 values from standard normal distribution were simulated and the support was determined by minimum and maximum simulated values, $I=[x_{min},x_{max}]$. Data were collected in a form of histogram, where the breakpoints are equidistantly spaced. The representative points of the histogram cells are denoted as asterisks in Figure 13 (left) and these points form the input data $(x_{i},y_{i})$, $i=1,\ldots,n$ for smoothing purpose. The $y$-values stand for discretized relative contributions to the overall probability mass, therefore discrete clr transformation [1] is needed to obtain a real vector with zero sum constraint (Figure 13, right). These data points were smoothed using the procedure from Section 3 and back-transformed to the original space. In the concrete setting $k=2$, $l=1$, $\alpha=0.5$, $\Delta\lambda=\{x_{1},-2,-1,0,1,2,x_{n}\}$, $w_{i}=1$, $\forall i=1,\ldots,n$ were considered. The resulting spline $s_{2}(x)$ with zero integral on interval $[x_{min},x_{max}]$ is plotted in Figure 13 (right). In the left plot the compositional spline $\xi_{2}(x)$ with unit integral by using (5) is depicted.

6 Application to anthropometric data

For the purpose of illustrating the smoothing procedure outlined in Section 3, a real-world data set dealing with the most commonly used anthropometric measure relating to body weight is presented. The data set we consider collects the body weight of apparently healthy Czech adolescents and young adults aged 15-31 years (the total of 4436 records) which were recruited non-randomly by offering free body composition assessment. Body weight was measured by the InBody 720 device (Biospace Co., Ltd, Seoul, Korea), recorded as the total body mass rounded to the nearest 0.1 kg.

The raw data for each of $N=16$ age groups, i.e. $\left[15,16\right),\left[16,17\right),\ldots$, $\left[30,31\right)$, were turned into a form of histogram data as follows. The sampled values of the body weight in each age group were divided into equally-spaced classes of the united range 40-110 kg and the optimal number of classes, denoted by $q_{i}$, $i=1,\ldots,16$, was set according to the well-known Sturge’s rule separately across the age groups. Because of the insufficient number of sampled data for males and females in each age group, gender information was ignored. Although there might be some differences between male and female samples due to sex dimorphism, they are not that dramatic (e.g., contrary to body height) because the weight is influenced also by external factors (nutrition, physical activity), and still allow for a reasonable aggregation of data. Subsequently, the proportions in classes within each age group were computed and present zero-values (zero counts in the respective classes) were imputed by values $(2/3)\cdot(1/n_{i}),i=1,\ldots,N$ according to [13], where $n_{i}$ stands for the number of observations in $i$th age group. Finally, the raw discretized density data $f_{i,j},i=1,\ldots,N,j=1,\ldots,q_{i}$ which correspond to the midpoints $t_{i,j}$ of classes, $i=1,\ldots,N,j=1,\ldots,q_{i}$ (i.e., $f_{i,j}=f(t_{i,j}))$, were obtained by dividing (not necessary normalized) proportions $\textbf{p}_{i}=\left(p_{i,1},\ldots,p_{i,q_{i}}\right)^{\top}$ of counts in classes by the length of the respective intervals resulting from partition of the weight range in each of age groups. Figure 14 shows four examples of histograms with different number of classes together with raw data (Table 1) to be smoothed. To do so, their transformation into real vectors is needed. We note that if the histogram data are constructed on subintervals of the same length, i.e. with equally-spaced breakpoints, it enables to use the discrete version of the clr transformation [1] directly on the vector of proportions $\textbf{p}_{i},i=1,\ldots,N$ by considering the scale invariance property. Otherwise, the input of the clr transformation must be vectors with raw density data $\textbf{f}_{i}=\left(f_{i,1},\ldots,f_{i,q_{i}}\right)^{\top},i=1,\ldots,N$ in order to avoid misleading results which would not reflect the actual behavior of data. Vectors of clr transforms are hereafter denoted as $\mbox{clr}(\textbf{f}_{i})=\left(\mbox{clr}(f_{i,1}),\ldots,\mbox{clr}(f_{i,q_{i}})\right)^{\top}$ for $i=1,\ldots,N$ and clr values are listed in Table 2.

Having the collected data $(t_{i,j},\mbox{clr}(f_{i,j}))$, we proceed to smooth them with the compositional smoothing splines using a system of $ZB$-spline basis functions from the $L^{2}_{0}(I)$. They are considered on domain $I=\left[40,107\right]$ which has been modified in order to avoid undesired artifacts in densities at their right-hand side. For all $N$ observations, the same strategy was followed to set the values of the input parameters for the smoothing procedure. We employed cubic smoothing splines ($k=3,l=2$) with the given sequence of knots $\Delta\lambda=\{40,62,84,107\}$, the vectors of weights $\textbf{w}_{i}$ for all input data equal to vectors of ones and the smoothing parameter $\alpha$ was set to $0.5$. That is, when minimizing the penalized functional (18), the same importance is assigned to both smoothness of the smoothing splines as well as to their approximative properties. The resulting compositional smoothing splines are obtained via their clr representation

[TABLE]

the corresponding $ZB$-spline coefficients are reported in Table 3.

Figure 15(a) displays an example of three raw density data from Figure 14 together with smoothed curves in the $L^{2}_{0}$ space (right) and after the inverse transformation (5) in the $\mathcal{B}^{2}$ space (left). The whole sample of smoothed density functions (Figure 15(b)) is plotted on blue scale distinguishing the age groups. Two trends are apparent – in the younger age groups, the estimated density functions are right skewed and exhibit lower variability while with increasing age they become more symmetric followed by higher variability. Nevertheless, density function in age group $\left[23,24\right)$ does not fully respect this behavior: the variability trend holds, but the distribution of weights is more similar to those in younger age groups as it is skewed more to the left. In general, it seems that adolescents appear to be predominantly of a lower body weight than the older persons whose weight is more spread over the weight classes and more pronounced in the middle part of the distribution. Accordingly, there is also a higher incidence of higher weights in comparison with the younger adolescents.

Importantly, the quality of smoothing is the same irrespective of the shape of the distribution. It might just happen that the smoothed densities exhibit heavier tails although they are not indicated by the data $f_{i,j}$, even with some artifacts typical for overfitting. This is obviously due to keeping the zero integral constraint, which is more sensitive to deviations from monotonic character of densities at their tails. A possible way out, applied here, was to reduce slightly the range from $I=\left[40,110\right]$ to $I=\left[40,107\right]$ in order to keep predominantly the monotonic behavior of (normalized) counts $f_{i,j}$ at the right tails.

Of course, the smoothed data can be further analyzed using methods of functional data analysis [18], adapted in order to respect specific properties of densities. It is demonstrated here for case of the compositional functional principal component analysis (SFPCA) [9]. This statistical tool has been recently designed based on the Bayes space methodology, so it enables to capture the main modes of relative variability in a data set consisting of sampled density functions. Given a data set of $N$ zero-mean functional observations $X_{1},\ldots,X_{N}$ in $\mathcal{B}^{2}(I)$, SFPCA aims to find the (normalized and orthogonal) directions of maximum variability in dataset, i.e., a collection of density functions $\left\{\theta_{\kappa}\right\}_{\kappa\geq 1}$ called simplicial functional principal components (SPFCs) which maximize the following objective function over $\theta\in\mathcal{B}^{2}(I)$,

[TABLE]

where the orthogonality condition is assumed to be fulfilled for $K\geq 2$. Since $\left\langle X_{i},\theta\right\rangle_{\mathcal{B}}$ represents a projection of $X_{i}$ along the direction $\theta\in\mathcal{B}^{2}(I)$, in fact we look for orthogonal basis functions in $\mathcal{B}^{2}(I)$ maximizing the relative variability of these projections. The maximization task (27) is efficiently implemented by applying the clr transformation (3) and the output can be back-transformed from the $L^{2}(I)$ to the $\mathcal{B}^{2}(I)$ space, as detailed in [9].

The interpretation of SFPCs can be performed by displaying:

• •

individual SFPCs (as clr transformed density functions, SFPCs always represent contrasts between the parts of domain $I$);

• •

overall mean density function $\bar{X}$ along with its perturbation by SFPCs powered by a suitable constant,

[TABLE]

where ${\rho}_{\kappa}$ is an amount of the variability of the dataset along the direction $\textit{x}_{\kappa}$ and it holds ${\rho}_{1}\geq{\rho}_{2}\geq\ldots$. This is a natural choice because SFPCs represent variation around the overall mean density function;

• •

the projection of dataset along the directions $\theta_{\kappa}$,

[TABLE]

where $\textit{x}_{i\kappa}=\left\langle X_{i},\theta_{\kappa}\right\rangle_{\mathcal{B}},i=1,\ldots,N$ are so called principal component scores associated with the $\kappa$th SFPCs $\theta_{\kappa}$. The scores can be plotted for pairs of the first SFPCs to assess the relationship among sampled density functions or to reveal presence of outlying observations;

• •

to complete the above interpretation it is important to note that the original functions $X_{1},\ldots,X_{N}$ and the eigenfunctions $\theta_{\kappa},\kappa\geq 1$ are based on a $CB$-spline basis expansion

[TABLE]

respectively, by considering (26).

SFPCA is also a statistical method for reducing dimensionality of dataset. The number of SFPCs can be determined from the scree plot which displays cumulative percentage of the total variance explained by each subsequent SFPC. That is, the dimensionality is identified by a point in scree plot at which explained variability drops off.

For the actual computation, $CB$-splines (29) of the input functional observations, represented by the corresponding $ZB$-splines (26), were expressed by $B$-splines with $B$-spline coefficients (listed in Table 4) using formula (16). The output of SFPCA for body weight density functions is reported in Figure 16. According to scree plot (Figure 16(b)), two or three SFPCs should be taken, but we resort to use only first two of them which capture together almost 85$\%$ of the total variability of the data set. The first SFPC (Figure 16(c)) represents the contrast between the weight below and above 78 kg, which could be considered as a reference (average) weight. Hence, higher scores along the SFPC1 are expected for age groups with higher incidence of individuals with higher body weight than the average, and, conversely, lower scores are associated with age groups with higher incidence of individuals with lower body weight then the average. The interpretation of SFPC1 can be obviously linked with age, see Figure 16(b). The scree plot more or less separates rather right skewed weight density functions of younger age groups (located on the left in the scree plot) from those more symmetric ones associated with older age groups (located on the right in the scree plot).

The second SFPC (Figure 16(d)) characterizes the variability within the tails of density functions, i.e. the main contribution to the variability along SFPC2 is provided by the lowest and highest weight values ($\leq$ 51 kg and $\geq$ 98 kg respectively). It contrasts low and high weights (associated with high scores along the SFPC2) against middle weight values (associated with low scores along the SFPC2), see Figure 16(b). The consistent interpretation can be also observed from Figure 16(e) which displays the variation along the first two directions – SFPC1 and SFPC2 – with respect to sample mean $\bar{f}(t),t\in I$ (i.e. $\bar{f}\oplus/\ominus 2\sqrt{\rho_{\kappa}}\odot\text{SFPC}_{\kappa},{\kappa}=1,2$).

Figures 16(f) and 16(g), respectively, represent two main modes of variability in the data set ($\left\langle f_{i},\theta_{\kappa}\right\rangle_{\mathcal{B}}\odot\theta_{\kappa},\kappa=1,2,i=1,\ldots,N$). For instance, the variation along SFPC2 is confirmed to be exhibited in tails of density functions and the observations with lowest (gold curve) and highest score (red curve) further support the conclusions made so far. The high scores along the second direction thus reflect heavier tails and, conversely, the low scores along the second direction reflects low incidence of individuals with extreme (both small and high) weights. Nevertheless, the relationship of scores (Figure 16(b)) is apparent: at the beginning, they continue to fall, reach a bottom and then continue to grow. The relationship might be partially explained by unequal representation of men and women in age groups and unequal number of observations in these age groups. Another reason might be that data corresponding to age groups with low SFPC2 scores were collected mostly from students of the Faculty of Physical Culture at Palacký University in Olomouc, Czech Republic, which form more homogeneous population than an average one. In any case, the second SFPC reveals an interesting feature which is worth to be further investigated.

7 Conclusions

The compositional splines, which enable to construct a spline basis in the clr space of density functions ($ZB$-spline basis) and consequently also in the original space of densities ($CB$-spline basis), might become an important contribution within the Bayes space methodology for processing of functional data carrying relative information. They provide a solid theoretical base for further developments of the approximation theory in context of the Bayes spaces, but even more importantly, compositional splines can be used also for adaptation of popular methods of functional data analysis for density functions. Here the case of compositional functional principal component analysis was presented, but similarly, e.g., regression analysis or classification methods could be developed. Also further tuning of the compositional splines is possible, here represented by the smoothing compositional splines or by orthonormalization of the $ZB$-basis. The latter case be used for an orthogonal projection of a density function on a subset of $CB$-splines, to further applications within the approximation theory or also for development of the theoretical framework of functional data analysis.

The pending challenge is to generalize the methodology introduced above also to $p$ dimensional density functions, $p>1$, which can be formally extended from any univariate density $f(x),\,x\in I=[a,b]$ to $f(\mathbf{x})$, where $\mathbf{x}=(x_{1},\ldots,x_{p})^{\top}\in\mathcal{I}=I_{1}\times\ldots\times I_{p}=[a_{1},b_{1}]\times\ldots\times[a_{p},b_{p}]$, in Equations (1) to (6); $\eta=b-a$ would be replaced by $H=\prod_{i=1}^{p}(b_{i}-a_{i})$. Currently an approach which focuses on keeping the zero integral constraint of the clr transformed densities was developed in [8] as a generalization of [12], which, however, does not lead to a compositional counterpart of the $B$-spline basis. A consistent approach in this direction is currently under development.

Acknowledgements

The authors gratefully acknowledge both the support by Czech Science Foundation GA18-09188S, the grant IGA_PrF_2018_024 Mathematical Models of the Internal Grant Agency of the Palacký University in Olomouc, and the grant COST Action CRoNoS IC1408.