# Principal component analysis in Bayes spaces for sparsely sampled density
functions

Lisa Steyer & Sonja Greven

###### Abstract

This paper presents a novel approach to functional principal component
analysis (FPCA) in Bayes spaces in the setting where densities are the object
of analysis, but only few individual samples from each density are observed.
We use the observed data directly to account for all sources of uncertainty,
instead of relying on prior estimation of the underlying densities in a two-
step approach, which can be inaccurate if small or heterogeneous numbers of
samples per density are available. To account for the constrained nature of
densities, we base our approach on Bayes spaces, which extend the Aitchison
geometry for compositional data to density functions. For modeling, we exploit
the isometric isomorphism between the Bayes space and the $\mathbb{L}^{2}$
subspace $\mathbb{L}_{0}^{2}$ with integration-to-zero constraint through the
centered log-ratio transformation. As only discrete draws from each density
are observed, we treat the underlying functional densities as latent variables
within a maximum likelihood framework and employ a Monte Carlo Expectation
Maximization (MCEM) algorithm for model estimation. Resulting estimates are
useful for exploratory analyses of density data, for dimension reduction in
subsequent analyses, as well as for improved preprocessing of sparsely sampled
density data compared to existing methods. The proposed method is applied to
analyze the distribution of maximum daily temperatures in Berlin during the
summer months for the last 70 years, as well as the distribution of rental
prices in the districts of Munich.

##  1 Introduction

A classic task in statistics is to estimate the underlying density from sample
data, since density functions can be used to describe the distribution of
real-valued random variables. However, when not all observations are
identically distributed, but constitute repeated draws from a set of density
functions $f_{1},\dots,f_{n}$, e.g. for $n$ different individuals,
distributional properties of these density functions themselves may be the
actual target of a statistical analysis. Examples in which densities or
distributions are considered the observational units exist in many different
fields. These include the size distributions of different zooplankton in
oceanology (Nerini and Ghattas, 2007), the distributions of firm size in
econometrics (Huynh and Jacho-Chávez, 2010), mortality densities at different
locations in epidemiology (Scimone et al., 2021), and the densities of glucose
levels among several diabetes patients in medical research (Matabuena et al.,
2021).

Probability density functions can be seen as a special case of functional
data, but considering them as the unit of observation poses two major
challenges. First, any density function must be non-negative and integrate to
one to be valid. Second, in practice, density functions are often unobserved
and accessible only through discrete samples. That is, for each density
function $f_{i}$, $i=1,\dots,n$ there is usually only an independent and
identically distributed sample $x_{ij}\sim f_{i}$, $j=1,\dots,m_{i}$
available. Our goal is to develop a Principal Component Analysis (PCA) for
densities that can handle both of these challenges. A PCA for density data is
of interest for several reasons. First, resulting estimates are useful for
exploratory analyses to better understand the main modes of variation in
density data. Second, the resulting dimension reduction allows to succinctly
describe differences and trends in densities and the corresponding principal
components (PCs) can be used as a parsimonious data-driven basis in subsequent
analyses, as common in functional data analysis (Yao et al., 2005; Chiou and
Li, 2007; Scheipl et al., 2015). Third, the reconstructed densities resulting
from the PCA can also be used as improved preprocessing of sparsely sampled
density data compared to existing methods, if subsequent analysis methods need
observed or reconstructed density data as input (Scimone et al., 2021; Maier
et al., 2022).

Existing research has primarily addressed one of the two challenges associated
with studying densities as functions, while overlooking the other. Initially,
Kneip and Utikal (2001) used functional principal component analysis (FPCA) in
the unbounded space $\mathbb{L}^{2}$ of quadratic integrable functions without
considering the density constraints. They did, however, account for discretely
observed data by estimating the covariance surface based on the combined
observations from all densities. In recent years, researchers have directed
their attention towards incorporating the intrinsic geometric constraints in
the space of density functions. The following discussion provides an overview
of their work. However, when dealing with discretely observed data, their
primary approach has been to estimate the observed densities through
preprocessing steps, such as aggregating the data using histograms or kernel
density estimates, and ignoring the reconstruction uncertainty in the further
analysis.

To address the density constraints, several metrics have been considered for
the space of probability density functions. In particular, the Wasserstein
distance (Panaretos and Zemel, 2019) for probability measures is widely used,
while the Fisher-Rao metric (Srivastava et al., 2007) imposes a manifold
structure on the space of density functions. Although statistical analysis can
be performed directly on manifolds (e.g. geodesic PCA in Wasserstein space
(Bigot et al., 2013)), it is often more convenient to map the densities to a
space with a simpler structure, perform statistical analysis there, and then
back-transform the results to the original density space. Various
transformations have been considered, such as the log-hazard and log-quantile
density transformations (Petersen and Müller, 2016) for mapping the density
functions to a Hilbert space. In particular, Hron et al. (2016) used the
centered log-ratio (clr) transformation to obtain FPCA for densities. The clr
transformation is particularly useful here because it defines a one-to-one
mapping between the squared-log integrable (proper and improper) density
functions and the separable Hilbert space $\mathbb{L}^{2}_{0}$, which
represents the space of square-integrable functions that integrate to zero.
This means that the clr transformation also induces a Hilbert space structure
on the space of squared-log integrable (proper and improper) density
functions, which is called Bayes Hilbert space (Egozcue et al., 2006; van den
Boogaart et al., 2014).

This paper develops Functional Principal Component Analysis (FPCA) in the
Bayes spaces for the setting where only samples from the set of densities of
interest are available. Our approach utilizes the observed data directly for
calculations instead of estimating the underlying densities beforehand in a
two-step approach. Estimating density $f_{i}$ becomes particularly challenging
when dealing with small sample sizes $m_{i}$ in each observation unit. To
address this, a two-step approach is proposed by Qiu et al. (2022) for cases
with heterogeneous sample sizes. Firstly, the underlying process is estimated
based on a subsample where each density is densely observed. Then, the
remaining densities are estimated using the process estimated in the first
step as a prior. However, this approach is only feasible when there exists a
representative subset of densities that have been densely observed.

The advantage of our approach is its applicability even when all densities are
sparsely observed. To achieve this, we incorporate the observed data of all
densities in a maximum likelihood framework, treating the underlying densities
as latent variables. To estimate the model, we utilize the Expectation-
Maximization (EM) algorithm (Dempster et al., 1977). As the expectation step
in our framework is not analytically solvable, requiring a Monte Carlo
approximation, we base the estimation of our model on the so-called Monte
Carlo Expectation-Maximization (MCEM) adaptation of the EM algorithm (Wei and
Tanner, 1990).

We organize our contributions as follows. First, in Section 2, we introduce
the latent density model and develop an MCEM algorithm suitable for this
scenario. Then, in Section 3, we apply our methodology in two different
applications. First, in the context of climate change and changing (extreme)
temperatures in particular, we study how the distribution of maximum
temperatures during the summer months in Berlin has evolved over the last 70
years. Second, we look at the distribution of rental prices in the different
districts of Munich. In Section 4, we use a simulation to demonstrate that our
method is particularly effective for analyzing densities when few observations
drawn from them are available. Finally, we conclude the paper with a
discussion in Section 5.

##  2 Principal component analysis (PCA) for densities based on individual
samples

We consider densities as elements of the Bayes Hilbert space, which has proven
to be a valuable framework for modeling densities. In order to perform
Principal Component Analysis (PCA) in the Bayes Hilbert space, which is a
separable Hilbert space, we employ the Karhunen-Loève decomposition. While PCA
was originally developed as a dimension reduction tool for finite-dimensional
data, the Karhunen-Loève decomposition extends this concept to infinite-
dimensional Hilbert spaces (Hsing and Eubank, 2015).

###  2.1 PCA in Bayes Hilbert spaces

We start with reviewing the structure of the Bayes Hilbert space. For
simplicity and as it seems natural in most applications, we restrict ourselves
in this work mainly to densities with respect to the Lebesgue measure
$\lambda$ defined on a compact interval $I\subset\mathbb{R}$ although the
construction can be done for general measures (van den Boogaart et al., 2014).
In Subsection 2.4, we briefly discuss how the case of compositional data can
be treated in a similar way.

######  Theorem 2.1 (Bayes Hilbert space (Egozcue et al., 2006)).

Let $B=\\{f=\exp(g)|g\in\mathbb{L}^{2}(I)\\}$ and consider the equivalence
relation $f_{1}\sim f_{2}\Leftrightarrow\exists\alpha>0:f_{1}=\alpha f_{2}$
for $f_{1},f_{2}\in B$. Denote by $\mathcal{B}=B/_{\sim}$ the set of
equivalence classes $[f]$ with $f\in B$. Then $\mathcal{B}$ equipped with the
operations

  * $\oplus$

(addition) given by the pertubation operator $[f_{1}]\oplus[f_{2}]=[f_{1}\cdot
f_{2}]$ for all $[f_{1}],[f_{2}]\in\mathcal{B}$,

  * $\odot$

(scalar multiplication) given by the powering operation
$\alpha\odot[f]=[f^{\alpha}]$ for all $[f]\in\mathcal{B},\alpha\in\mathbb{R}$
and

  * $\langle\cdot,\cdot\rangle_{\mathcal{B}}$

the scalar product given via
$\langle[f_{1}],[f_{2}]\rangle_{\mathcal{B}}=\frac{1}{2|I|}\int_{I}\int_{I}\log\left(\frac{f_{1}(x)}{f_{1}(y)}\right)\log\left(\frac{f_{2}(x)}{f_{2}(y)}\right)dx\
dy$ for all $[f_{1}],[f_{2}]\in\mathcal{B}$

is a separable Hilbert space.

As a separable Hilbert space the Bayes Hilbert space is isometrically
isomorphic to any other infinite dimensional separable Hilbert space, in
particular $\mathcal{B}$ is isometrically isomorphic to $\mathbb{L}^{2}_{0}$,
the space of square-integrable functions integrating to zero, via the centered
log-ratio transformation.

######  Lemma 2.2 (Centered log-ratio transformation).

The centered log-ratio (clr) transformation

$\displaystyle\operatorname{clr}:\mathcal{B}\to\mathbb{L}^{2}_{0},\quad[f]\mapsto\log(f)-\frac{1}{|I|}\int_{I}\log(f(x))dx$
(1)

is a bijective isometry with the inverse given via
$\operatorname{clr}^{-1}(g)=[\exp(g)]$ for all $g\in\mathbb{L}_{0}^{2}$.

In particular, the clr transformation is well defined as
$\operatorname{clr}([\alpha f])=\operatorname{clr}([f])$ for all $f\in B$ and
$\alpha\in\mathbb{R}$. A detailed proof for Lemma 2.2 can be found in Appendix
A.1 (compare with van den Boogaart et al. (2014) for densities on general
measure spaces). This identification of the Bayes space $\mathcal{B}$ with
$\mathbb{L}_{0}^{2}$ allows statistical modeling to be performed in
$\mathbb{L}_{0}^{2}$ instead of directly in $\mathcal{B}$, allowing the
application of techniques developed for functional data. In particular, Hron
et al. (2016) used this correspondence to introduce PCA in Bayes space, while
Scimone et al. (2021) and Maier et al. (2022) exploit its use for regression
purposes.

Similarly, we will achieve a principal component decomposition of observed
densities $f_{1},\dots,f_{n}$, $n\in\mathbb{N}$ via considering them being the
back-transforms of realizations $g_{1},\dots,g_{n}$ of a stochastic process
$\mathcal{G}=\\{G(x)\\}_{x\in I}\subset\mathbb{L}_{0}^{2}$ characterized by
its mean function $\mu(x)=\mathbb{E}(G(x))$ and covariance kernel
$K(x_{1},x_{2})=\text{Cov}(G(x_{1}),G(x_{2}))$ for all $x,x_{1},x_{2}\in I$.
The Karhunen-Loève decomposition (Karhunen, 1946; Loève, 1946) then yields the
functional principal component representation

$\displaystyle G(x)=\mu(x)+\sum_{k=1}^{\infty}Z_{k}\varphi_{k}(x)$ (2)

where $\varphi_{k}$, $k\in\mathbb{N}$ are the orthonormal eigenfunctions of
the covariance operator $\mathbb{L}_{0}^{2}\to\mathbb{L}_{0}^{2}$,
$g\mapsto\int_{I}K(x_{1},\cdot)g(x_{1})dx_{1}$ associated with the covariance
kernel $K$ and uncorrelated principal component scores $Z_{k}$ of decreasing
importance, with $\mathbb{E}(Z_{k})=0$ and $\text{Var}(Z_{k})=\sigma_{k}^{2}$
the corresponding eigenvalues, $\sigma^{2}_{1}\geq\sigma^{2}_{2}\dots\geq 0$.
For more details on this decomposition for second-order stochastic processes
refer to Hsing and Eubank (2015).

For a given sample of (fully observed) functions $g_{1},\dots,g_{n}$,
$n\in\mathbb{N}$ the unknown parameters $\varphi_{k},Z_{k}$ of this process
could then be estimated via computing the eigendecomposition of the sample
covariance operator associated with the sample covariance kernel
$\hat{K}_{n}(x_{1},x_{2})=\frac{1}{n}\sum_{i=1}^{n}(g_{i}(x_{1})-\hat{\mu}(x_{1}))(g_{i}(x_{2})-\hat{\mu}(x_{2}))$
with $\hat{\mu}(x)=\frac{1}{n}\sum_{i=1}^{n}g_{i}(x)$ for all
$x,x_{1},x_{2}\in I$. Thus, the eigenfunctions $\varphi_{k}$ could be
estimated as the eigenfunctions $\hat{\varphi}_{k}$ of the sample covariance
and the distribution of $Z_{k}$ as the empirical distribution of the factor
loadings $z_{ik}=\int_{I}(g_{i}(x)-\hat{\mu}(x))\hat{\varphi}_{k}dx$,
$i=1,\dots,n$ for all $k=1,\dots,N$ where $N$ is the number of non-zero
eigenvalues of the sample covariance.

The correspondence of $\mathbb{L}_{0}^{2}$ and the Bayes Hilbert space via the
clr transformation (Lemma 2.2) gives an analog principal component
decomposition to Equation (2) for densities. The process
$\\{\operatorname{clr}^{-1}(G)(x)\\}_{x\in I}$ is given as

$\displaystyle\operatorname{clr}^{-1}(G)(x)=\operatorname{clr}^{-1}(\mu)(x)\oplus\bigoplus_{k=1}^{\infty}\exp(Z_{k})\odot\operatorname{clr}^{-1}(\varphi_{k})(x)$
(3)

inheriting the properties of the decomposition of $\mathcal{G}$ in
$\mathbb{L}_{0}^{2}$ via the clr transformation (Equation (1)). Namely, we
obtain orthonormal eigenfunctions since
$\langle\operatorname{clr}^{-1}(\varphi_{k}),\operatorname{clr}^{-1}(\varphi_{l})\rangle_{\mathcal{B}}=\langle\varphi_{k},\varphi_{l}\rangle_{\mathbb{L}_{2}}=0$
for all $k\neq l$ and
$\|\varphi\|_{\mathcal{B}}=\|\varphi\|_{\mathbb{L}_{2}}=1$ for all
$k\in\mathbb{N}$ and principal components scores $\exp(Z_{k}),k\in\mathbb{N}$
with $Z_{k}\perp Z_{l}$ and $\mathbb{E}(Z_{k})=0$ for all $k,l\in\mathbb{N}$.
If we additional assume that $G$ is Gaussian, this also implies that the
principal component scores $\exp(Z_{k}),k\in\mathbb{N}$ are independent (and
therefore also uncorrelated) and we can compute their expectation as the
evaluation of the moment generating function, i.e.
$\mathbb{E}(\exp(Z_{k}))=\exp(0.5\sigma_{k}^{2})$, where
$\sigma^{2}_{k}=\text{Var}(Z_{k})$.

Note that elements in $\mathcal{B}$ are equivalence classes with respect to
scalar multiplication. As a result, they either possess a unique
representative which is a proper density function, or solely consist of
improper densities. In practical applications, the focus is often on modeling
proper density functions. Consequently, it becomes necessary to impose
restrictions on $\mathbb{L}_{0}^{2}$ accordingly. Specifically, only those
functions $g\in\mathbb{L}_{0}^{2}$ that satisfy $\int_{I}\exp(g(x))dx<\infty$
are suitable for representing proper densities. This requirement can be
effectively met, for example, by using spline representations for elements of
$\mathbb{L}_{0}^{2}$ (Machalová et al., 2021; Maier et al., 2022) since the
exponential transformation of a spline is bounded, guaranteeing integrability
over $I$. Furthermore, if we consider a finite set of bounded density
functions $f_{1},\dots,f_{n}$, such as those obtained from kernel density
estimates or histograms, and apply principal component decomposition in the
Bayes space (3), we are modeling only functions in the span of the data, which
are also bounded. Consequently, in such cases, the eigenfunctions (and any
linear combinations of them) are inherently proper densities.

Each equivalence class $[f]$ in the subset of $\mathcal{B}$ constructed such
that it contains only proper densities can be uniquely identified with the
element $\tilde{f}\in[f]$ that satisfies $\int_{I}\tilde{f}(x)dx=1$. To
simplify the notation, albeit with a slight misuse, we will refer to the
equivalence class by this particular element, $\tilde{f}=[f]$ in the following
discussion. Hence, if $[\exp(g)]$ contains proper density functions, i.e. if
$\int_{I}\exp(g(x))dx$ is finite, we denote by

$\displaystyle\operatorname{clr}^{-1}(g)=\frac{\exp(g)}{\int_{I}\exp(g(x))dx}$
(4)

the back-transformed element in $\mathcal{B}$ under the inverse clr
transformation.

###  2.2 Likelihood formulation assuming latent densities

We have seen in the previous section that the correspondence of
$\mathbb{L}_{0}^{2}$ and the Bayes Hilbert space via the clr transformation
provides a convenient approach for conducting PCA on fully observed densities,
since it allows estimation via first transforming the density functions to the
Hilbert space $\mathbb{L}^{2}$ using the clr transformation, and then
performing the principal component decomposition in this well-known function
space. This procedure for fully observed density functions was previously
suggested by Hron et al. (2016). However, their approach motivates the so-
called simplicial principal component analysis solely as a maximization
problem. In other words, they seek to find the projections of the observed
densities that maximize the variance along their directions. While this leads
to the same decomposition for fully observed density functions when the
covariance kernel is estimated using the sample covariance kernel, the
stochastic process perspective becomes especially useful in the more common
scenario we have in mind.

We focus on analyzing densities that are neither directly observable nor can
be satisfactorily estimated from observed data as a preprocessing step.
Instead, we assume that we have access to samples $x_{i1},\dots,x_{im_{i}}$,
where $m_{i}\in\mathbb{N}$, drawn from each probability distribution with
density $f_{i}$, where $i=1,\dots,n$. Our objective is then to conduct
Maximum-Likelihood estimation for the parameters $\mu$ and $K$ of the
underlying process $\mathcal{G}$ in $\mathbb{L}_{0}^{2}$ based on these
samples. By estimating these parameters, we can compute the eigenvalues and
eigenfunctions of the estimated covariance operator, which allows us to obtain
the principal component decomposition which directly yields the principal
component decomposition in the Bayes space via the inverse clr transformation.
To accomplish this, we need to make a distributional assumption for $G$. For
simplicity, we assume that $G$ follows a Gaussian process with a finite
Karhunen-Loève decomposition. More precisely, we assume the following model.

######  Definition 2.3 (Latent density model).

Let $GP(\mu,K)$ be a Gaussian process with mean function $\mu$ and covariance
kernel $K$ taking values in a finite dimensional subspace
$\mathcal{H}\subset\mathbb{L}_{0}^{2}$ consisting of bounded functions. Then
we assume the following data generating process

$\displaystyle
X_{ij}\overset{i.i.d.}{\sim}\operatorname{clr}^{-1}(G_{i})=\frac{\exp(G_{i})}{\int_{I}\exp(G_{i}(x))dx}$

with $G_{i}$ being independent replicates of $GP(\mu,K)$ for all
$i=1,\dots,n$, $n\in\mathbb{N}$, $j=1,\dots,m_{i}$.

Assuming a finite dimensional subspace $\mathcal{H}$ is not restrictive in
practice, where observed data always lies in a finite dimensional subspace.
This assumption allows us to apply maximum likelihood theory, since in this
case the parameters $\mu$ and $K$ of the underlying Gaussian process can be
described by a finite set of real-valued parameters and the likelihood has a
similar form as the marginal likelihood of a mixed model. More precisely, we
can formulate the following corollary (proof in Appendix A.2).

######  Corollary 2.4.

For any orthonormal basis $e_{1},\dots,e_{N}\subset\mathbb{L}^{2}$ with
$\mathcal{H}\subseteq\text{span}\\{e_{1},\dots,e_{N}\\}$ we have that
$G\overset{i.i.d.}{\sim}GP(\mu,K)$ is equivalent to
$G=\sum_{k=1}^{N}\theta_{k}e_{k}$ with
$\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{N})^{T}\overset{i.i.d.}{\sim}\mathcal{N}(\boldsymbol{\nu},\boldsymbol{\Sigma})$
for $\mu=\sum_{k=1}^{N}\nu_{k}e_{k}$ and
$K(x_{1},x_{2})=\sum_{k=1}^{N}\sum_{l=1}^{N}e_{k}(x_{1})e_{l}(x_{2})\Sigma_{kl}$,
where $\boldsymbol{\Sigma}=\left(\Sigma_{kl}\right)_{k,l=1,\dots,N}$ and
$\boldsymbol{\nu}=(\nu_{1},\dots,\nu_{N})^{T}$.
If $\boldsymbol{v}_{1},\dots,\boldsymbol{v}_{N}$ are the eigenvectors of
$\boldsymbol{\Sigma}$ with corresponding eigenvalues
$\sigma^{2}_{1},\dots,\sigma^{2}_{N}$ then
$\varphi_{l}=\sum_{k=1}^{N}v_{lk}e_{k}$, $l=1,\dots,N$ are the eigenfunctions
of the covariance operator given by the covariance function $K$ with the same
eigenvalues $\sigma^{2}_{1},\dots,\sigma^{2}_{N}$, where
${\boldsymbol{v}}_{l}=({v}_{l1},\dots,{v}_{lN})$ for all $l=1,\dots,N$.

With this equivalence, the latent density model 2.3 can also be written as

$\displaystyle
X_{ij}\overset{i.i.d.}{\sim}\operatorname{clr}^{-1}(G_{i})=\frac{\exp(G_{i})}{\int_{I}\exp(G_{i}(x))dx}\quad\text{with}\quad
G_{i}=\sum_{k=1}^{N}\theta_{ik}e_{k}\text{ and
}\boldsymbol{\theta}_{i}=(\theta_{i1},\dots,\theta_{iN})\overset{i.i.d.}{\sim}\mathcal{N}(\boldsymbol{\nu},\boldsymbol{\Sigma})$
(5)

and estimation of the parameters $\mu$ and $K$ is equivalent to estimation of
$\boldsymbol{\nu}$ and $\boldsymbol{\Sigma}$. Note that we do not assume
span{$e_{1},\dots,e_{N}\\}=\mathcal{H}$, we only need to cover $\mathcal{H}$,
since for any finite basis $\\{e_{1},\dots,e_{N}\\}$ in $\mathbb{L}^{2}$ the
sum to zero constrain carries over to the parameters $\boldsymbol{\nu}$ and
$\boldsymbol{\Sigma}$. That means if $\int_{I}e_{k}(x)dx=\int_{I}e_{l}(x)dx$
for all $k,l=1,\dots,N$, one just needs to ensure that the entries in
$\boldsymbol{\nu}$ as well as the entries of all eigenvectors
$\boldsymbol{v}_{1},\dots,\boldsymbol{v}_{N}$ of $\boldsymbol{\Sigma}$ sum to
zero.

Corollary 2.4 also implies that the maximum likelihood estimators will be
asymptotically consistent, if model 2.3 is correctly specified, that is the
image $\mathcal{H}$ of the the process $G$ is actually in
$\text{span}\\{e_{1},\dots,e_{N}\\}$. An interesting question for future
research is whether in the case of misspecification, in particular when the
image space $\mathcal{H}$ is not finite, a sequence of finite-dimensional
subspaces of $\mathbb{L}^{2}$ can be constructed such that the corresponding
estimators are asymptotically consistent in this case as well.

To estimate the finite dimensional parameters $\boldsymbol{\nu}$ and
$\boldsymbol{\Sigma}$ of the latent density model via maximum likelihood, we
need to maximise the likelihood function given the realizations
$\boldsymbol{x}_{i}=(x_{i1},\dots,x_{im_{i}})^{T}$ from the random sample
$\boldsymbol{X}_{i}=(X_{i1},\dots,X_{im_{i}})^{T}$, $i=1,\dots,n$. Let
$G_{i}=\sum_{k=1}^{N}\theta_{ik}e_{k}$ and
$\boldsymbol{\theta}_{i}=(\theta_{i1},\dots\theta_{iN})^{T}$ for all
$i=1,\dots,n$. Then the marginal likelihood for the parameters
$\boldsymbol{\nu}$ and $\boldsymbol{\Sigma}$ is given as

$\displaystyle
L(\boldsymbol{\nu},\boldsymbol{\Sigma}|\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n})=\prod_{i=1}^{n}\int_{\mathbb{R}^{N}}\frac{\exp(\sum_{j=1}^{m_{i}}\sum_{k=1}^{N}\theta_{ik}e_{k}(x_{ij}))p(\boldsymbol{\theta}_{i}|\boldsymbol{\nu},\boldsymbol{\Sigma})}{\left(\int_{I}\exp(\sum_{k=1}^{N}\theta_{ik}e_{k}(x))dx\right)^{m_{i}}}d\boldsymbol{\theta}_{i}.$
(6)

For a detailed derivation, please refer to Appendix A.3. Maximizing this
marginal likelihood can be seen as an empirical Bayes approach, where the
prior for $\boldsymbol{\theta}_{i}$ is a multivariate normal distribution with
mean $\boldsymbol{\nu}$ and covariance $\boldsymbol{\Sigma}$. Note that by $p$
we denote a general density function, for example here
$p(\boldsymbol{\theta}_{i}|\boldsymbol{\nu},\boldsymbol{\Sigma})$ denotes the
density of a multivariate normal distribution with parameters
$\boldsymbol{\nu}$ and $\boldsymbol{\Sigma}$.

Due to the complicated nature of the likelihood function and the potential
abundance of parameters in $\boldsymbol{\nu}$ and $\boldsymbol{\Sigma}$,
numerical optimization of (6) is challenging. Therefore, in the following
section, we use the Monte Carlo Expectation Maximization (MCEM) algorithm as a
numerical method to effectively tackle this maximization problem.

###  2.3 Model estimation using an MCEM algorithm

The EM algorithm, developed by Dempster et al. (1977), addresses the challenge
of maximum likelihood estimation in the presence of incomplete or missing
data. This algorithm provides a framework for estimating the parameters of
statistical models that involve unobserved or latent variables. It iteratively
updates parameter estimates by incorporating both observed data and estimates
of the missing data. In our specific context, we want to use the EM algorithm
to handle latent, unobserved densities $f_{1},\dots,f_{n}$, along with
observed values $\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}$ sampled from
these densities. Central to the EM algorithm is the notion of the expected
complete-data log-likelihood, which in our case becomes

$\displaystyle
Q(\boldsymbol{\nu},\boldsymbol{\Sigma}|\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$
$\displaystyle=\mathbb{E}(\log(p(\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n},\boldsymbol{\theta}_{1}\dots,\boldsymbol{\theta}_{n}|\boldsymbol{\nu},\boldsymbol{\Sigma})))$
$\displaystyle=\mathbb{E}(\log(p(\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}|\boldsymbol{\theta}_{1}\dots,\boldsymbol{\theta}_{n})))+\mathbb{E}(\log(p(\boldsymbol{\theta}_{1}\dots,\boldsymbol{\theta}_{n}|\boldsymbol{\nu},\boldsymbol{\Sigma})))+\text{const.}$
$\displaystyle=\sum_{i=1}^{n}\mathbb{E}(\log(p(\boldsymbol{\theta}_{i}|\boldsymbol{\nu},\boldsymbol{\Sigma})))+\text{const.}$
(7)

where the expectation is taken with respect to the conditional distribution
$\boldsymbol{\theta}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}$
of the parameters
$\boldsymbol{\theta}_{i}=(\theta_{i1},\dots,\theta_{iN})^{T}$ of the latent
densities for all $i=1,\dots,n$ given the current estimates
$\boldsymbol{\nu}^{(h)}$ and $\boldsymbol{\Sigma}^{(h)}$ for
$\boldsymbol{\nu}$ and $\boldsymbol{\Sigma}$. However, in the given setting,
the conditional distribution needed to compute the expectation of the
complete-data log-likelihood is not directly available, making the computation
intractable. To address this challenge, Wei and Tanner (1990) introduce the
Monte Carlo Expectation Maximization (MCEM) algorithm, which uses a Monte
Carlo approach to approximate the expected value in $Q$. Thus, for our
particular use case, we need to generate samples of
$\boldsymbol{\theta}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}$
for all $i=1,\dots,n$. In the following, we outline the procedure for
obtaining these samples and implementing the MCEM algorithm in our use case.
Subsections 2.3.1 and 2.3.2 detail the E- and the M-steps, respectively, 2.3.3
the selection of model space and initial values, and 2.3.4 summarizes the
complete algorithm.

####  2.3.1 E-step

For all $i=1,\dots,n$ we approximate the conditional expectation
$\mathbb{E}(\log(p(\boldsymbol{\theta}_{i}|\boldsymbol{\nu},\boldsymbol{\Sigma})))$
where the expectation is taken with respect to
$\boldsymbol{\theta}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}$
using importance sampling, which is a method for estimating properties of a
target distribution by sampling from another, auxiliary distribution. In our
case we sample for all $i=1,\dots,n$, $r$ replicates of the parameters
$\boldsymbol{\theta}_{i}$ of the latent densities from an auxiliary
distribution with density $p_{i}^{*}(\boldsymbol{\theta}_{i})$ instead of
$\boldsymbol{\theta}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}$.
This yields replicates
$\boldsymbol{\theta}_{i}^{(1)},\dots,\boldsymbol{\theta}_{i}^{(r)}$,
$r\in\mathbb{N}$, and we approximate

$\displaystyle\mathbb{E}(\log(p(\boldsymbol{\theta}_{i}|\boldsymbol{\nu},\boldsymbol{\Sigma})))\approx\sum_{t=1}^{r}\frac{\omega_{it}}{\sum_{t=1}^{r}\omega_{it}}\log(p(\boldsymbol{\theta}_{i}^{(t)}|\boldsymbol{\nu},\boldsymbol{\Sigma}))$
(8)

with weights $\omega_{it},t=1,\dots,r$ given as
$\omega_{it}=\frac{p(\boldsymbol{\theta}_{i}^{(t)}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})}{p_{i}^{*}(\boldsymbol{\theta}_{i}^{(t)})}$
for all $i=1,\dots,n$. For details on this method and a comprehensive
treatment of several related Monte Carlo methods, see the book by Hammersley
and Handscomb (1964).

The key to this method lies in selecting an appropriate auxiliary
distribution. To achieve this, we use the eigen decompostion of
$\boldsymbol{\Sigma}^{(h)}$ with sorted eigenvalues
${\sigma_{1}^{2}}^{(h)}\geq\dots\geq{\sigma_{N}^{2}}^{(h)}$ and corresponding
eigenvectors $\boldsymbol{v}_{1}^{(h)},\dots,\boldsymbol{v}_{N}^{(h)}$. Then
$\boldsymbol{\theta}_{i}\sim\mathcal{N}(\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$
is equivalent to
$\boldsymbol{z}_{i}=\boldsymbol{V}^{(h)}(\boldsymbol{\theta}_{i}-\boldsymbol{\nu}^{(h)})\sim\mathcal{N}(\boldsymbol{0},\text{diag}({\sigma_{1}^{2}}^{(h)},\dots,{\sigma_{N}^{2}}^{(h)}))$,
where
$\boldsymbol{V}^{(h)}=(\boldsymbol{v}_{1}^{(h)},\dots,\boldsymbol{v}_{N}^{(h)})$
is a matrix whose columns are the eigenvectors of $\boldsymbol{\Sigma}^{(h)}$.
Thus sampling from
$\boldsymbol{\theta}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}$
is equivalent to sampling from the conditional distribution of the scores
$\boldsymbol{z}_{i}=(z_{i1},\dots,z_{iN})$ given as

$\displaystyle
p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$
$\displaystyle\propto
p(\mathbf{x}_{i}|\mathbf{z}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})p(\mathbf{z}_{i}|\boldsymbol{\Sigma}^{(h)})=p(\mathbf{x}_{i}|\boldsymbol{\theta}_{i}=\boldsymbol{V}^{(h)T}\mathbf{z}_{i}+\boldsymbol{\nu}^{(h)})\prod_{k=1}^{N}p(z_{ik}|{\sigma_{k}^{2}}^{(h)})$
$\displaystyle=\prod_{j=1}^{m_{i}}\operatorname{clr}^{-1}\left(\sum_{k=1}^{N}\nu_{k}^{(h)}e_{k}+\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}\right)(x_{ij})\prod_{k=1}^{N}p(z_{ik}|{\sigma_{k}^{2}}^{(h)})$
$\displaystyle=\frac{\exp\left(\sum_{j=1}^{m_{i}}\left(\mu^{(h)}(x_{ij})+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x_{ij})\right)\right)}{\left(\int_{I}\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)dx\right)^{m_{i}}}\prod_{k=1}^{N}p(z_{ik}|{\sigma_{k}^{2}}^{(h)})$
(9)

where $z_{ik}|{\sigma_{k}^{2}}^{(h)}\sim\mathcal{N}(0,{\sigma_{k}^{2}}^{(h)})$
for all $i=1,\dots,n$ and $k=1,\dots,N$. Here,
$\mu^{(h)}=\sum_{k=1}^{N}\nu_{k}^{(h)}e_{k}$ is the current estimate for the
mean function and
$g_{i}=\mu^{(h)}+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}$,
$i=1,\dots,n$ are the current predictions for the latent clr transformed
densities. Here we again take the Bayesian perspective, where
$p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$
is the posterior distribution for the prior
$\mathcal{N}(\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$. Note that
this posterior is a proper distribution if the prior
$\mathcal{N}(\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$ is proper,
that is if all eigenvalues of $\boldsymbol{\Sigma}^{(h)}$ are finite.

######  Lemma 2.5.

Let ${\sigma_{1}^{2}}^{(h)}<\infty$. Then
$\int_{\mathbb{R^{N}}}p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})d\mathbf{z}_{i}<\infty$
for all $i$.

A proof for this statement can be found in Appendix A.4. This also implies
that
$p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$
is decreasing as $\|\mathbf{z}_{i}\|\to\infty$ and since it is also
continuous, it attains its mode $\mathbf{z}_{i}^{*}\in\mathbb{R}^{n}$. This is
not necessarily the case if
$\mathcal{N}(\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$ is improper
(see Appendix A.5 for a counterexample).

Since we assume here that the prior distribution
$\mathcal{N}(\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$ is proper, the
mode of the posterior distribution
$\mathbf{z}_{i}^{*}=\operatorname{argmax}_{\mathbf{z}_{i}\in\mathbb{R}^{N}}p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$
is attained. Hence we can choose a multivariate normal distribution centered
around the mode as an auxiliary distribution for the scores $\mathbf{z}_{i}$.
We further choose the variances to be proportional to the prior variances
${\sigma_{1}^{2}}^{(h)},\dots,{\sigma_{N}^{2}}^{(h)}$. This means that for a
tuning parameter $\lambda>0$, we choose the auxiliary distribution
$p_{i}^{*}(\boldsymbol{z}_{i})$ to be
$\mathcal{N}(\boldsymbol{z}_{i}^{*},\lambda\text{diag}({\sigma_{1}^{2}}^{(h)},\dots,{\sigma_{N}^{2}}^{(h)}))$.
Usually we set $\lambda=1$, but if one wants to explore the parameter space
for $\mathbf{z}$ more or less, one can also set $\lambda$ larger or smaller.
Consequently, once we compute the mode $\mathbf{z}_{i}^{*}$, sampling from the
auxiliary distribution $p_{i}^{*}$ reduces to independently sampling each
element of the vector $\mathbf{z}_{i}$ from a univariate normal distribution.
To numerically compute the mode, i.e., the maximizer of (2.3.1) with respect
to $\mathbf{z}_{i}$, it is useful to derive the gradient of its log
transformation to apply a gradient descent algorithm.

######  Lemma 2.6.

The gradient of the log conditional density of the scores
$\mathbb{R}^{N}\to\mathbb{R}$,
$\mathbf{z}_{i}\mapsto\log(p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}))$
is given as

$\displaystyle\nabla\log(p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}))$
$\displaystyle=\sum_{k=1}^{N}\boldsymbol{v}^{(h)}_{k}\left(\sum_{j=1}^{m_{i}}e_{k}(x_{ij})-m_{i}\langle
f_{\mathbf{z}_{i}},e_{k}\rangle_{\mathbb{L}_{2}}\right)-\left(\frac{z_{il}}{{\sigma_{l}^{2}}^{(h)}}\right)_{l=1,\dots,N}$

where
$f_{\mathbf{z}_{i}}=\operatorname{clr}^{-1}\left(\mu^{(h)}+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}\right)$
for all $\mathbf{z}_{i}=(z_{i1},\dots,z_{iN})^{T}\in\mathbb{R}^{N}$.

A detailed derivation can be found in Appendix A.6. With this readily
available gradient, finding the mode becomes numerically feasible and we can
obtain i.i.d. samples for the scores $\mathbf{z}_{it}$ and corresponding
weights $\omega_{it}\in\mathbb{N}$ for all $t=1,\dots,r$ using the importance
sampling described above.

The equivalence of conditionally sampling $\boldsymbol{\theta}_{i}$ or
$\mathbf{z}_{i}$ also yields samples $\boldsymbol{\theta}_{i}^{(t)}$ from
$\boldsymbol{\theta}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}$
for all $t=1,\dots,r$ via
$\boldsymbol{\theta}_{i}^{(t)}=\boldsymbol{\nu}^{(h)}+\boldsymbol{V}^{(h)T}\mathbf{z}_{it}$.
Hence, we approximate the expected complete-data log-likelihood given in (2.3)
by

$\displaystyle
Q(\boldsymbol{\nu},\boldsymbol{\Sigma}|\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})\approx\sum_{i=1}^{n}\sum_{t=1}^{r}\frac{\omega_{it}}{\sum_{t=1}^{r}\omega_{it}}\log(p(\boldsymbol{\theta}_{i}^{(t)}|\boldsymbol{\nu},\boldsymbol{\Sigma}))+\text{const.}$
(10)

using the Monte-Carlo approximation given in (8) for the conditional
expectation.

####  2.3.2 M-step

The M-step updates the parameters $\boldsymbol{\nu}$ and $\boldsymbol{\Sigma}$
by maximizing this approximation (10) of the $Q$-function. The new estimates
are then given by

$\displaystyle\left(\boldsymbol{\nu}^{(h+1)},\boldsymbol{\Sigma}^{(h+1)}\right)=\operatorname{argmax}_{\boldsymbol{\nu},\boldsymbol{\Sigma}}Q(\boldsymbol{\nu},\boldsymbol{\Sigma}|\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})\approx\operatorname{argmax}_{\boldsymbol{\nu},\boldsymbol{\Sigma}}\sum_{i=1}^{n}\sum_{t=1}^{r}\frac{\omega_{it}}{\sum_{t=1}^{r}\omega_{it}}\log(p(\boldsymbol{\theta}_{i}^{(t)}|\boldsymbol{\nu},\boldsymbol{\Sigma})).$

Since for the latent density model (5) we assume that
$\boldsymbol{\theta}_{i}^{(t)}|\boldsymbol{\nu},\boldsymbol{\Sigma}$ follows a
multivariate normal distribution with mean $\boldsymbol{\nu}$ and covariance
matrix $\boldsymbol{\Sigma}$, this optimization problem corresponds to a
weighted maximum likelihood estimation of the mean and the covariance matrix.
Remarkably, this maximization problem also arises when the EM algorithm is
used to estimate a Gaussian mixture distribution, allowing us to derive the
solution based on the computations performed for this problem (e.g. Bishop and
Nasrabadi, 2006). Hence, we compute the updates for the parameters of our
model as

$\displaystyle\boldsymbol{\nu}^{(h+1)}$
$\displaystyle=\frac{1}{\sum_{i=1}^{n}\sum_{t=1}^{r}\omega_{it}}\sum_{i=1}^{n}\sum_{t=1}^{r}\omega_{it}\boldsymbol{\theta}_{i}^{(t)}$
$\displaystyle\boldsymbol{\Sigma}^{(h+1)}$
$\displaystyle=\frac{1}{\sum_{i=1}^{n}\sum_{t=1}^{r}\omega_{it}}\sum_{i=1}^{n}\sum_{t=1}^{r}\omega_{it}(\boldsymbol{\theta}_{i}^{(t)}-\boldsymbol{\nu}^{(h+1)})(\boldsymbol{\theta}_{i}^{(t)}-\boldsymbol{\nu}^{(h+1)})^{T}.$

These are the weighted mean and weighted covariance matrix of the samples of
the principal component scores $\boldsymbol{\theta}_{i}^{(t)}$, $i=1,\dots,n$,
$t=1,\dots,r$.

####  2.3.3 Selection of model space and initial values

In order to apply our method to real-world problems, we first need to find a
suitable model space as well as suitable initial values for the MCEM
algorithm. We suggest using piecewise constant spline functions for modeling.
Since the piecewise constant functions are dense in $\mathbb{L}^{2}$, they
allow us to approximate any function in $\mathbb{L}_{0}^{2}$, and thus any
density in $\mathcal{B}$, with arbitrary accuracy if the nodes are chosen to
be on a fine grid. Therefore, we fix a fine grid for the knots
$\kappa_{1},\dots,\kappa_{N+1}$ and choose as a basis the indicator functions
which are one between two neighboring knots and zero elsewhere. That is
$e_{k}=\mathds{1}_{[\kappa_{k},\kappa_{k+1}]}$ for all $k=1,\dots,N$.

To obtain suitable initial values for $\boldsymbol{\nu}^{(0)}$ and
$\boldsymbol{\Sigma}^{(0)}$, we propose to first estimate the latent densities
$\hat{f}_{1},\dots,\hat{f}_{n}$ by kernel density estimation. We then develop
their clr transformations $\hat{g}_{1},\dots,\hat{g}_{n}$ in our basis
$e_{1},\dots,e_{N}$. Subsequently, we estimate $\boldsymbol{\nu}^{(0)}$ and
$\boldsymbol{\Sigma}^{(0)}$ as the empirical mean and covariance of the
coefficients $\boldsymbol{\theta}_{1}\dots,\boldsymbol{\theta}_{n}$ of
$\hat{g}_{1},\dots,\hat{g}_{n}$, respectively. This approach effectively
restricts the model space to the span of the kernel density estimates. This
could be a problem, for example, if only a small sample of densities is
available, or if the kernel density estimates are close to zero in some parts
of the support.

In this case, an alternative would be to initially select a lower dimensional,
smooth spline space for modeling. If we choose an orthonormal basis
$e_{1},\dots,e_{N}\in\mathbb{L}_{0}^{2}$, such as normalized versions of the
orthogonal compositional splines suggested by Machalová et al. (2021), we can
choose arbitrary values for the initial mean and covariance of the
coefficients, for instance,
$\boldsymbol{\nu}^{(0)}=\boldsymbol{0}\in\mathbb{R}^{N}$ and
$\boldsymbol{\Sigma}^{(0)}=\mathbb{I}_{N}\in\mathbb{R}^{N\times N}$, the
identity matrix.

On the other hand, when the number of densities $n$ is large, which results
also in a $N=n$ basis functions by the procedure given above, not only does
computing the mode become a high-dimensional optimization problem, which is
computationally demanding, also calculating the weights $\omega_{it}$ becomes
unstable as in this case typically many of the variances
${\sigma_{k}^{2}}^{(h)}$ will be very small, thus the product
$\prod_{k=1}^{N}p(z_{k}|{\sigma_{k}^{2}}^{(h)})$ will be close to zero. In
this case, we suggest reducing the dimensionality of
$GP(\hat{\mu}^{(h)},\hat{K}^{(h)})$ in each step $h$ by setting the variances
${\sigma_{k}^{2}}^{(h)}=0$ for $k>N^{\prime}$, where $N^{\prime}$ is a chosen
value such that $N^{\prime}\ll N$. The specific value of $N^{\prime}$ can be
determined based on the desired proportion of variance explained, while the
proportion of variance explained should be considerably greater than desired
for the final PCA.

####  2.3.4 Estimation of the density PCA using the MCEM algorithm

Subsections 2.3.1-2.3.3 derived all necessary ingredients to now use the MCEM
algorithm to estimate the PCA. First, we initialize $\boldsymbol{\nu}^{(0)}$
and $\boldsymbol{\Sigma}^{(0)}$ according to 2.3.3. Then we iterate the
following two steps until the convergence criteria
$\|\boldsymbol{\nu}^{(h+1)}-\boldsymbol{\nu}^{(h)}\|_{\mathbb{L}^{2}}<\epsilon$
and
$\|\boldsymbol{\Sigma}^{(h+1)}-\boldsymbol{\Sigma}^{(h)}\|_{\mathbb{L}^{2}}<\epsilon$
for a threshold $\epsilon>0$ are reached.

\- E-step following 2.3.1,

\- M-step following 2.3.2.

Estimates at convergence then approximate the maximum likelihood estimates of
the observed data likelihood, that is
$\hat{\boldsymbol{\nu}}\approx\boldsymbol{\nu}^{(h+1)}$ and
$\hat{\boldsymbol{\Sigma}}\approx\boldsymbol{\Sigma}^{(h+1)}$ . Finally, the
equivalence given in Corollary 2.4 yields estimates for $\mu$, $\varphi_{k}$
and $\sigma^{2}_{k}$ in (2) using the eigendecomposition of
$\hat{\boldsymbol{\Sigma}}$, where
$\hat{\boldsymbol{v}}_{1},\dots,\hat{\boldsymbol{v}}_{N}$ are the eigenvectors
of $\boldsymbol{\Sigma}$ with corresponding eigenvalues
$\hat{\sigma}^{2}_{1},\dots,\hat{\sigma}^{2}_{N}$. Then, the estimate for the
mean is obtained as $\hat{\mu}=\sum_{k=1}^{N}\hat{\nu}_{k}e_{k}$, where
$\hat{\boldsymbol{\nu}}=(\hat{\nu}_{1},\dots,\hat{\nu}_{N})^{T}$. The
estimates for the eigenfunctions are given as
$\hat{\varphi}_{l}=\sum_{k=1}^{N}\hat{v}_{lk}e_{k}$, where
$\hat{\boldsymbol{v}}_{l}=(\hat{v}_{l1},\dots,\hat{v}_{lN})$ for all
$l=1,\dots,N$ with corresponding eigenvalues
$\hat{\sigma}^{2}_{1},\dots,\hat{\sigma}^{2}_{N}$.

The scores $Z_{ik}$, $k=1,\dots,N$ for each density $f_{i},i=1,\dots,n$ are
then predicted as the posterior mode given the estimates
$\hat{\boldsymbol{\nu}}$ for the mean and $\hat{\boldsymbol{\Sigma}}$ for the
covariance of the coefficients, that is
$\hat{\mathbf{z}}_{i}=\operatorname{argmax}_{\mathbf{z}_{i}\in\mathbb{R}^{N}}p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\hat{\boldsymbol{\nu}},\hat{\boldsymbol{\Sigma}})$
with $\hat{\mathbf{z}}_{i}=(\hat{z}_{i1},\dots\hat{z}_{iN})$. This also yields
predictions for the latent densities as
$\hat{f}_{i}=\operatorname{clr}^{-1}(\hat{g}_{i})$ with
$\hat{g}_{i}=\hat{\mu}+\sum_{k=1}^{N}\hat{z}_{ik}\hat{\varphi}_{k}$ for all
$i=1,\dots,n$.

###  2.4 PCA for compositional data as a special case

Although we have limited our considerations so far to densities with respect
to the Lebesgue measure, it’s important to recognize that our methodology can
be seamlessly extended to densities with respect to arbitrary measures. In the
following, we illustrate this via showing how our approach can be used to
derive Principal Component Analysis (PCA) for compositional data, namely
densities with respect to the discrete measure. Notably, our method has the
advantage of being applicable to “sparsely observed” compositional data, known
as count compositions in the compositional data literature (Filzmoser et al.,
2018), even when some of the categories have no observations and without
imputations.

The discrete measure on the power set $\mathcal{P}(\\{A_{1},\dots,A_{k}\\})$
of a finite set of disjoint outcomes $A_{1},\dots,A_{k}$ is given as
$\eta=\sum_{k=1}^{N}\delta_{A_{k}}$, where $\delta_{A_{k}}(B)=1$, if $B=A_{k}$
and $\delta_{A_{k}}(B)=0$ else, for all
$B\in\mathcal{P}(\\{A_{1},\dots,A_{k}\\})$. Hence, every probability measure
on $\mathcal{P}(\\{A_{1},\dots,A_{k}\\})$ is given by a discrete density, i.e.
probability mass function $f:\\{A_{1},\dots,A_{k}\\}\to\mathbb{R}$ with
respect to $\eta$ by the Radon-Nikodym Theorem. Here, the density $f$ is
characterized solely by the values $\pi_{k}=f(A_{k})$ for all $k=1,\dots,N$,
and since we consider only probability measures, it must hold that
$\sum_{k=1}^{N}\pi_{k}=1$. That means $\mathcal{B}$, the set of densities with
respect to the discrete measure on $\mathcal{P}(\\{A_{1},\dots,A_{k}\\})$ can
be identified with the simplex
$\\{\boldsymbol{\pi}\in\mathbb{R}^{N}|\sum_{k=1}^{N}\pi_{k}=1,\pi_{k}\geq 0\
\forall k=1,\dots,N\\}$ and via the discrete centered log-ratio transformation
$\boldsymbol{\rho}=\operatorname{clr}(\boldsymbol{\pi})=(\log(\pi_{1})-\frac{1}{N}\sum_{k=1}^{N}\log(\pi_{k}),\dots,\log(\pi_{N})-\frac{1}{N}\sum_{k=1}^{N}\log(\pi_{k}))$
with the $N-1$ dimensional Hilbert space
$\mathcal{H}=\mathbb{R}^{N}_{0}=\\{\boldsymbol{\rho}\in\mathbb{R}^{N}|\sum_{k=1}^{N}\rho_{k}=0\\}$
(Aitchison, 1982). If we equip this space with the standard scalar product on
$\mathbb{R}^{N}$, we obtain the Aitchison geometry on $\mathcal{B}$, which
defines a Hilbert space structure for compositional data, i.e. for densities
with respect to the discrete measure.

In order to use our method to perform principal component analysis for
observed count compositions $\boldsymbol{\pi_{1}},\dots,\boldsymbol{\pi_{n}}$,
we need to find a suitable basis representation of $\mathcal{H}$, as well as
appropriate initial values for the mean and covariance of the corresponding
basis coefficients. For densities with respect to the discrete measure, kernel
density estimation cannot be employed due to the usual lack of an order
relation and thus neighborhood structure to smooth over. We therefore rely on
the alternative described in Subsubsection 2.3.3 and use an orthonormal basis
of $\mathcal{H}=\mathbb{R}^{N}_{0}$.

There are many options since an orthonormal basis can be obtained from any
basis of $\mathbb{R}^{N}_{0}$ using Gram-Schmidt orthogonalization. For
example, Egozcue et al. (2003) suggest to use the following such basis

$\displaystyle
e_{k}=\sqrt{\frac{k}{k+1}}(\overbrace{k^{-1},\dots,k^{-1}}^{k\text{
times}},-1,0,\dots,0)^{T}$

with $k=1,\dots,N$ as an orthonormal basis of $\mathbb{R}^{N}_{0}$, which
yields an orthonormal basis of $\mathcal{B}$ via the inverse clr transform.
With this basis choice we propose to set
$\boldsymbol{\nu}^{(0)}=\boldsymbol{0}\in\mathbb{R}^{N-1}$ and
$\boldsymbol{\Sigma}^{(0)}=\mathbb{I}_{N}\in\mathbb{R}^{(N-1)\times(N-1)}$ as
initial values for the mean and covariance of the basis coefficients and
proceed with the MCEM algorithm as in the continuous case described above.

##  3 Applications

In this section, we demonstrate the applicability and advantages of our latent
density-based PCA. We consider two different applications. In the first
application, we analyze densities describing distributions of summer daily
maximum temperatures in Berlin. In the second application, the latent
densities describe the distribution of rental prices for each district in
Munich.

###  3.1 Distributions of daily maximum temperature in summer months per year

According to the Copernicus Climate Change Service (C3S) report
(https://climate.copernicus.eu/esotc/2022/temperature), 2022 was the second
warmest year on record in Europe, with temperatures 0.9°C above the long-term
average. In particular, the summer of that year set a new mark as the hottest
on record, with temperatures 1.4°C above average, topping the previous warmest
summer in 2021 by 0.3°C.

These findings are consistent with the observed trend of increasing
temperatures, an indicator of climate change. To refine the analysis of summer
temperature trends, we focus not on average temperatures, but on trends in the
entire distribution of daily maximum temperatures during the summer months of
June, July, and August each year. That is, we consider as observational units
densities that describe the distributions of daily maximum temperatures in
summer. Thus, for every year, we treat the daily maximum temperature measured
for the 92 days in June, July, and August as observations from these densities
per year. The daily maximum temperature data we use in this application have
been collected from 1951 to 2022 at a single weather observatory, which is
Berlin Tempelhof. Data and metadata are available at https://www.ecad.eu,
provided by the ECA&D project (Klein Tank et al., 2002).

Figure 1: Latent density PCA for daily maximum temperature (°C). Top: Effect
of adding/subtracting $\hat{\sigma}_{k}\hat{\varphi}_{k}$ to the clr
transformed mean density $\mu$, where $\hat{\varphi}_{k}$ is the $k$th
principal component, with corresponding eigenvalue $\hat{\sigma}_{k}^{2}$,
$k=1,2,3,4$. Middle: Effect on the density level, i.e.
$\operatorname{clr}^{-1}$ transformations of the functions in the top row.
Bottom: Temporal trend of the corresponding predicted scores per year, with
scatterplot smoother and pointwise confidence bands overlaid.

To visualize and describe how the distribution of daily maximum temperature
has changed over the period from 1951 to 2022, we proceed as follows. First,
we obtain a low-dimensional representation of the latent daily maximum
temperature densities using our latent density PCA and with kernel density
estimates as initial estimates (Figure 6 in the appendix) For technical
details of the estimation please refer to Appendix B.1). Then, in Figure 1, we
visualize the first four principal components on clr level (top row), and
transformed back to density level (middle row). In the bottom row, we plot the
temporal trend of the corresponding predicted scores, overlaying a scatterplot
smoother and pointwise confidence bands based on Wood (2017).

This shows that adding a multiple of the first principal component, which
explains 40% of the total variability in latent densities, to the estimated
mean mainly causes a rightward shift of the density. In particular, a positive
value of the first principal component implies that high temperatures are more
likely than in average years. Looking at the effect on clr level, we notice
that this is especially true for temperatures above 35°C, which are more
likely to occur in years with high first principal component scores. Notably,
these first principal component scores show a clear increase over the time
course from 1951 to 2022 (Figure 1, bottom left), meaning that hot and also
very hot daily maximal temperatures in summer became more likely.

In contrast, the scores associated with the second and fourth principal
components show no or almost no visible temporal trend. Adding these principal
components to the estimated mean results in smaller changes in the shape of
the density, with subtler shifts to the right in certain areas of the density.
However, adding the third principal component to the estimated mean shifts the
density towards experiencing moderate to hot temperatures (25°C-35°C) more
frequently, while decreasing the likelihood of milder temperatures
(15°C-25°C). The scores associated with this third principal component also
show an increase over time.

This means that the trends of both the first and third principal component
scores indicate that hot and very hot days are becoming more likely. When we
plot both scores together (Figure 2) we see that all early years,
corresponding to (dark) blue points, tend to lie in the bottom left corner,
while recent years (yellow and orange points) are predominantly in the top
right corner. Thus, recent years have high first principal component scores
and/or high third principal component scores, which means the likelihood for
hot and/or very hot days has been higher in these years than in earlier years.

Figure 2: Scores associated with the first and third principal components for
all years.

This application shows that the estimation of the latent density model (5) is
suitable for identifying a small number of principal directions of variation
in the data for densities when only discrete observations of each density are
available. These principal components can then be used to visualize the data
and/or for further analysis, such as to relate them to other scalar variables
(in our case, the year of observation).

However, in this homogeneous and only mildly sparse setting with 92
observations per density, differences to the simpler two-step approach with
pre-smoothing (our starting values), a PCA obtained from the clr
transformations of the kernel density estimates (see Figure 7), are relatively
small. One notable difference, though, is that in the latent density model,
the trend of the score associated with the first principal component changes
more rapidly around 1990, with a steeper increase since that time, than that
of the first principal component scores obtained in the two-step approach. The
trend in the third PC score is somewhat more pronounced as well, which may be
due to unstable kernel density estimates in areas where the densities are
close to zero. In contrast to this first application, our second application
has more heterogeneous sample sizes and we will observe how this unbalanced
setting affects the estimation.

###  3.2 Distributions of rental prices in the districts of Munich

We consider data that was collected in Munich in 2019 (the most recent year)
to construct an official rent index by estimating the rent given certain
covariates of the apartments in a regression model. For details on the data
collection, the available variables, estimation and interpretation of the rent
index, see (Windmann and Kauermann, 2019, in German). A sub-dataset, which
contains data on all 3255 apartments but only a subset of the covariates used
for the official Munich rent index, is provided in the supplemental material
of Fahrmeir et al. (2023) and is also used here.

The goal of the analysis here is to describe the differences in the
distributions of rental prices across the districts of Munich. To do this, we
model the distribution of net rents per square meter, assuming a latent
density for each district. Figure 8 in Appendix B.2 displays the histogram
estimates and kernel density estimates (Gaussian kernel, bandwidth = 2) for
the densities in each district. Table 1 shows that the number of observations
$m_{i}$, $i=1,\dots,25$ per district varies considerably in this dataset. It
ranges from 29 observations in district 23-Allach-Untermenzing to 261
observations in district 9-Neuhausen-Nymphenburg. Below, we will compare the
estimation of our latent density model (5) for this heterogeneous sampling
scheme to a two-step approach, where the densities are first estimated and
then PCA is performed.

district $i$ |  name | $m_{i}$
---|---|---
1 |  Altstadt-Lehel | 79
2 |  Ludwigsvorstadt-Isarvorstadt | 217
3 |  Maxvorstadt | 219
4 |  Schwabing-West | 241
5 |  Au-Haidhausen | 230
6 |  Sendling | 136
7 |  Sendling-Westpark | 95
8 |  Schwanthalerhöhe | 110
9 |  Neuhausen-Nymphenburg | 261
10 |  Moosach | 81
11 |  Milbertshofen-Am Hart | 101
12 |  Schwabing-Freimann | 198
13 |  Bogenhausen | 165

district $i$ |  name | $m_{i}$
---|---|---
14 |  Berg am Laim | 75
15 |  Trudering-Riem | 92
16 |  Ramersdorf-Perlach | 111
17 |  Obergiesing-Fasangarten | 125
18 |  Untergiesing-Harlaching | 151
19 |  Thalkirchen-Obersendling-Forstenried-Fürstenried-Solln | 160
20 |  Hadern | 64
21 |  Pasing-Obermenzing | 107
22 |  Aubing-Lochhausen-Langwied | 35
23 |  Allach-Untermenzing | 29
24 |  Feldmoching-Hasenbergl | 38
25 |  Laim | 135

Table 1: Munich districts: number $i$, name and number of observations per
distric $m_{i}$

We first consider in Figure 3(a) the results of estimating the latent density
model with kernel density estimates as initial estimates. Details for the
estimation can be found in Appendix B.2. In this model, the first principal
component causes a shift towards more expensive apartments, and looking at the
effect on clr level, we see that this principal component also primarily
describes whether the occurrence of very expensive apartments (more than 25€
per square meter) is likely. Looking at which districts have a positive score
for the first principal component, we see that these districts are mainly in
the city center and correspond closely to the districts that have an
increasing influence on the expected net rent in the regression model
estimated in Fahrmeir et al. (2023).

The second principal component describes the presence of very cheap and
expensive apartments. However, this principal component has only a small
influence on the shape of the rent per square meter density. In the latent
density model the variance of the principal components, and thus their
importance, is measured at the clr level. Since in this case the largest
effect is on parts of the clr transformed densities with negative values, the
multiplicative effect on the actual densities is in parts where the density is
close to zero and thus the effect, corresponding to heavier tails, is hardly
visible. It is therefore not surprising that the spatial distribution of the
corresponding scores has little structure, since the second principal
component describes the occurrence of a few extreme observations in the
districts.

The effect of the third principal component on the densities appears to be
larger than the effect of the second principal component, although the
variability explained is smaller, as it affects areas around the mode of the
distribution. For this component, a negative score mostly describes a larger
share of affordable housing (5€ to 10€ per square meter). These negative
scores are mainly predicted for districts in the south of Munich, i.e.
6-Sendling and neighboring districts.

(a) PCA based on the latent density model.

(b) Two-step approach using kernel density estimates and then applying PCA
after clr transformation.

Figure 3: Comparison of the latent density model with a two-step approach
using kernel density estimates as preprocessing on the Munich rent dataset.
For each method, the first row shows the effect of adding a principal
component times the corresponding standard deviation to the mean on clr level,
the second row shows the same effect on density level, and the third row shows
a map of Munich districts, where the color represents the predicted scores.

Comparing the effects of the first three principal components estimated with
the latent density model with the estimates obtained using a two-step approach
based on pre-smoothing with kernel density estimates (Gaussian kernel,
bandwidth $=2$), we see that the estimates differ considerably. In particular,
for the two-step approach, the main variation in the scores is caused by
density estimates of districts with only a small number of observations
$m_{i}$. That is, the most extreme scores are estimated for district 22 with
$m_{22}=35$, district 23 with $m_{23}=29$, and district 24 with $m_{24}=38$,
for the first three principal components, respectively. These districts
benefit from the shrinkage effect of the latent density model, which divides
the total variance into a part which is due to the underlying stochastic
process for the latent densities and a part due to sampling from them. This
shrinkage causes the predictions for the densities in these districts to be
closer to the overall mean. Note also that the percentages of variance
explained for the PCs shown in figures 3(a) and 3(b), respectively, are
correspondingly relative to different total variances.

Figure 8 in the appendix shows that for these densities, the latent density
predictions appear more plausible than the kernel density estimates,
especially in areas where there are few observations, i.e. for very cheap or
expensive housing. However, nearly all predicted densities appear to better
reflect the underlying data using our approach, as the kernel density
estimates generally underestimate the modes. The reason for this behavior is
that we had to choose a relatively large bandwidth for the kernels to avoid
estimating close to zero densities in other parts of the domain, especially
for small samples. While we use the same kernel density estimates for our
approach as starting values, the influence of such problems and the choice of
the bandwidth in general is much smaller due to the later updates of the
model. In order to systematically investigate the influence of the number of
observations per density on the model estimation, a simulation is carried out
in the following subsection.

##  4 Simulation

This subsection aims to evaluate how well our latent density model can recover
the mean and covariance structure of the latent process for a varying number
of observations per density, ranging from very few observations to a moderate
number. We also compare the performance of the model with two-step approaches,
where the density estimates are obtained first, and then the PCA is performed
after applying the clr transformation to each density. This corresponds to the
simplicial PCA proposed by Hron et al. (2016). For the first comparison, we
obtain kernel density estimates with a Gaussian kernel and then perform PCA on
the clr transformed densities, which is also used as the initial estimate for
our method. Second, we use a compositional spline estimate for the clr
transformed densities, as suggested by Machalová et al. (2021), before
performing the PCA.

To this end we choose the following simulation setting. In each simulation run
we simulate $n=30$ densities on the interval $I=[0,1]$ from a Gaussian process
with true clr transformed mean function
$\mu(x)=-20(x-\frac{1}{2})^{2}+\frac{5}{3}$ and only two principal components.
These are given on clr level as $g_{1}(x)=\frac{1}{5}\sin(10(x-\frac{1}{2}))$
with corresponding factor $Z_{1}\sim\mathcal{N}(0,0.5)$ and
$g_{2}(x)=\frac{1}{10}\cos(2\pi(x-\frac{1}{2}))$ with corresponding factor
$Z_{2}\sim\mathcal{N}(0,0.2)$. Note that these functions satisfy
$\mu,g_{1},g_{2}\in\mathbb{L}^{2}_{0}$ and $g_{1}\perp g_{2}$. The samples for
the densities $f_{i}$, $i=1,\dots,30$ are then obtained as
$f_{i}=\operatorname{clr}^{-1}(\mu+z_{i1}g_{1}+z_{i2}g_{2})$ where
$z_{i1}\overset{i.i.d.}{\sim}Z_{1}$ and $z_{i2}\overset{i.i.d.}{\sim}Z_{2}$.
The resulting densities are shown in Figure 5 in the top row on the left for a
simulation run with $m_{i}=40$. Finally, we sample observations
$x_{ij}\overset{i.i.d}{\sim}f_{i}$ with $j=1,\dots,m_{i}$ from each density
$f_{i}$, $i=1,\dots,n$. For the number of observations per density we consider
$m_{i}\in\\{20,40,80,160\\}$ and repeat the simulation 100 times for each
$m_{i}$.

For the two-step approaches we then estimate the densities, either with kernel
density estimates using Gaussian kernels with bandwidths $0.12,0.09,0.08,0.07$
for the different setting with $m_{i}=20,40,80,160$, respectively, or using
cubic compositional splines with five knots. The kernel density estimates are
also used as initial estimates for our latent density model. The number of
Monte Carlo samples in the E step 2.3.1 is chosen to be $r=10h$, where $h$ is
the iteration index, i.e., the number of samples increases over the
iterations. The parameter $\lambda$ for the proposal density is set to $1$. In
each iteration, the dimension reduction is at most $0.0001$, i.e. we keep as
many principal components as necessary to explain at least 99.999 % of the
variance.

The performance of the different methods is evaluated as the distance to the
oracle estimates, that are the pointwise estimates of the mean and the
covariance function based on the true underlying densities $f_{1},\dots,f_{n}$
(Figure 4) evaluated on a equidistant grid with 200 grid points. More
precisely, we compute the distance of the mean estimates as
$\sqrt{\int_{0}^{1}(\tilde{\mu}(x)-\hat{\mu}(x))^{2}dx}$, where $\tilde{\mu}$
is the oracle estimate for the mean and $\hat{\mu}$ is the estimate of each
method. Analogously, the distance of the covariance functions is obtained as
$\sqrt{\int_{0}^{1}\int_{0}^{1}(\tilde{C}(x_{1},x_{2})-\hat{C}(x_{1},x_{2}))^{2}dx_{1}dx_{2}}$,
where $\tilde{C}$ is the oracle estimate for the covariance and $\hat{C}$ is
the estimate of each method.

Figure 4: $\mathbb{L}^{2}$ distance of the estimated mean and covariance
functions for each method to the oracle estimates based on the true underlying
densities $f_{1},\dots,f_{n}$.

As expected, the performance of all three methods improves as the number of
observations per density increases. Still, over all values of
$m_{i}\in\\{20,40,80,160\\}$, and for both the mean and the covariance
function, our latent density model has the smallest average distance to the
oracle estimate over the 100 replicates. This shows that our method
outperforms both two-step approaches in this scenario. However, when comparing
the two-step approaches, it is worth noting that kernel density estimates seem
to be better for estimating the mean, while compositional splines excel at
capturing the covariance structure. For a more concrete picture of how the
mean and covariance estimates for the three methods behave relative to the
oracle estimate, Figure 5 shows each for an example with $m_{i}=40$
observations per density.

Figure 5: One randomly selected simulation run with $m_{i}=40$ observations
per density. Here ’oracle’ refers to estimation based on the true underlying
densities. Top: Density estimates/predictions. Middle: Mean and first two
principal components. Bottom: Estimated covariance structure.

In this example, the estimates for the mean function, the first two principal
components, and the covariance function of the latent density model (second
column) appear similar to the corresponding oracle estimates (first column).
Also, the pattern of predicted latent densities (first row, third column) is
similar to the pattern of true underlying densities. For the two-step approach
based on kernel density estimates (third row), the modes of the estimated
densities appear to be too low, and although the mean has a similar shape as
the oracle mean, the first two principal components, as well as the covariance
function, differ substantially from the oracle estimates.

The compositional splines used for the two-stage approach in the fourth row
appear to behave similarly to the true densities near the center of the
observed interval (in $\approx[0.2,0.8]$), but show implausible
characteristics near the boundary, i.e., near $0$ and $1$. This is also
evident in the estimates for the mean and the first two principal components.
Similarly, the estimate for the covariance function seems to be close to the
oracle estimate in the center of the $I\times I$ domain, but not at the
boundaries.

##  5 Discussion

We have proposed improvements to existing PCA methods for densities in the
Bayes Hilbert space, by explicitly incorporating in a maximum likelihood
approach that there are usually only discrete samples of the densities
available. This differs from two-step approaches, where the densities are
first estimated in a preprocessing step and then PCA is performed for these
estimates in the Bayes Hilbert space, ignoring uncertainty arising from the
first density estimation step. We confirmed in applications and in a
simulation that our latent density model can be successfully estimated
employing an MCEM algorithm, and that the estimation of the mean and
covariance structure for the densities is superior to the estimation using
two-step approaches. While improvements are particularly pronounced for small
samples per density, differences persist even for moderately large samples.

Consequently, resulting estimates for the principal components using our
approach are better suited for understanding the variation in the underlying
densities, and the predicted scores can be better used for dimension reduction
and for further analyses, such as to describe differences and trends in the
densities. In addition, given the importance of principal components for
dimension reduction in functional data analysis (Ramsay and Silverman, 2005;
Chiou and Li, 2007), the principal components could also be used as the basis
for further subsequent functional data analysis methods, such as as a basis in
which to expand model terms in a functional regression model (as for example
in Yao et al., 2005; Scheipl et al., 2015; Volkmann et al., 2023) in the Bayes
Hilbert space. Finally, the predicted latent densities could also be used for
subsequent analysis, providing more reliable reconstructions of the underlying
true densities than the density estimates obtained using preprocessing with
kernel density estimation or compositional splines.

To allow maximum likelihood estimation of the principal components, we assume
a normal distribution of the scores, i.e. a Gaussian process (prior) for the
latent densities, as common in many statistical methods ranging from mixed
models (Guo, 2002) to Gaussian process regression (Rasmussen and Williams,
2005). This distributional assumption may be too restrictive in some
applications, however, in such cases leading to unjustified shrinkage of the
latent densities. In future research, it would thus be appealing to extend our
approach further by including a choice of different distributions for the
scores in the latent density model, e.g. to account for heavy tails.

In this paper, our primary focus was on continuous densities with respect to
the Lebesgue measure on a bounded interval. Furthermore, we provided a brief
outline of a potential extension to discrete measures, which incorporates
compositional data. In addition to an application of our model to the
compositional data case often encountered in practical data scenarios, an
interesting direction for future research would be a generalization of the
approach and implementation to other measures. This could include mixed
discrete and continuous measures as in Maier et al. (2022), the inclusion of
unbounded domains as well as the multivariate case.

## Acknowledgements

We gratefully acknowledge funding by grants GR 3793/3-1 ‘Flexible regression
methods for curve and shape data’ and GR 3793/8-1 ‘Flexible density regression
methods’ from the German Research Foundation (DFG). We would also like to
thank the data providers in the ECA&D project for providing the temperature
data for the first application and Michael Windmann for sharing the Munich
rent data for the second application.

## References

  * Aitchison (1982) J. Aitchison.  The statistical analysis of compositional data.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 44(2):139–160, 1982.  doi: https://doi.org/10.1111/j.2517-6161.1982.tb01195.x.
  * Bigot et al. (2013) Jérémie Bigot, Raúl Gouet, Thierry Klein, and Alfredo Lopez.  Geodesic PCA in the Wasserstein space.  _Annales de l’Institut Henri Poincaré, Probabilités et Statistiques_ , 53, 07 2013.  doi: 10.1214/15-AIHP706.
  * Bishop and Nasrabadi (2006) Christopher M Bishop and Nasser M Nasrabadi.  _Pattern recognition and machine learning_ , volume 4.  Springer, 2006.
  * Chiou and Li (2007) Jeng-Min Chiou and Pai-Ling Li.  Functional clustering and identifying substructures of longitudinal data.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 69(4):679–699, 2007.
  * Dempster et al. (1977) A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 39(1):1–38, 1977.  ISSN 00359246.  URL http://www.jstor.org/stable/2984875.
  * Egozcue et al. (2003) Juan José Egozcue, Vera Pawlowsky-Glahn, Glòria Mateu-Figueras, and Carles Barcelo-Vidal.  Isometric logratio transformations for compositional data analysis.  _Mathematical geology_ , 35(3):279–300, 2003\.
  * Egozcue et al. (2006) Juan José Egozcue, José Luis Díaz-Barrero, and Vera Pawlowsky-Glahn.  Hilbert space of probability density functions based on Aitchison geometry.  _Acta Mathematica Sinica_ , 22:1175–1182, 2006.
  * Fahrmeir et al. (2023) Ludwig Fahrmeir, Göran Kauermann, Gerhard Tutz, and Michael Windmann.  Spatial smoothing revisited: An application to rental data in Munich.  _Statistical Modelling_ , 0(0):online before print, 2023.
  * Filzmoser et al. (2018) Peter Filzmoser, Karel Hron, Matthias Templ, Peter Filzmoser, Karel Hron, and Matthias Templ.  Analyzing compositional data using R.  _Applied Compositional Data Analysis: With Worked Examples in R_ , pages 17–34, 2018.
  * Guo (2002) Wensheng Guo.  Functional mixed effects models.  _Biometrics_ , 58(1):121–128, 2002.
  * Hammersley and Handscomb (1964) J.M. Hammersley and D.C. Handscomb.  _Monte Carlo Methods_.  Methuen’s monographs on applied probability and statistics. Methuen, 1964\.  ISBN 9780416523409.  URL https://books.google.de/books?id=Kk4OAAAAQAAJ.
  * Hron et al. (2016) Karel Hron, Alessandra Menafoglio, Matthias Templ, Klara Hruzová, and Peter Filzmoser.  Simplicial principal component analysis for density functions in Bayes spaces.  _Comput. Stat. Data Anal._ , 94:330–350, 2016.
  * Hsing and Eubank (2015) Tailen Hsing and Randall Eubank.  _Theoretical foundations of functional data analysis, with an introduction to linear operators_ , volume 997.  John Wiley & Sons, 2015.
  * Huynh and Jacho-Chávez (2010) Kim Huynh and David Jacho-Chávez.  Firm size distributions through the lens of functional principal components analysis.  _Journal of Applied Econometrics_ , 25:1211–1214, 11 2010\.  doi: 10.1002/jae.1200.
  * Karhunen (1946) Kari Karhunen.  Zur Spektraltheorie stochasticher Prozesse.  In _Annales Academiae Scientiarum Fennicae Series A_ , volume 1, page 34, 1946.
  * Klein Tank et al. (2002) A. M. G. Klein Tank, J. B. Wijngaard, G. P. Können, R. Böhm, G. Demarée, A. Gocheva, M. Mileta, S. Pashiardis, L. Hejkrlik, C. Kern-Hansen, R. Heino, P. Bessemoulin, G. Müller-Westermeier, M. Tzanakou, S. Szalai, T. Pálsdóttir, D. Fitzgerald, S. Rubin, M. Capaldo, M. Maugeri, A. Leitass, A. Bukantis, R. Aberfeld, A. F. V. van Engelen, E. Forland, M. Mietus, F. Coelho, C. Mares, V. Razuvaev, E. Nieplova, T. Cegnar, J. Antonio López, B. Dahlström, A. Moberg, W. Kirchhofer, A. Ceylan, O. Pachaliuk, L. V. Alexander, and P. Petrovic.  Daily dataset of 20th-century surface air temperature and precipitation series for the European Climate Assessment.  _International Journal of Climatology_ , 22(12):1441–1453, 2002.
  * Kneip and Utikal (2001) Alois Kneip and Klaus Utikal.  Inference for density families using functional principal component analysis.  _Journal of the American Statistical Association_ , 96:519–542, 02 2001.  doi: 10.1198/016214501753168235.
  * Loève (1946) M Loève.  Fonctions aléatoires à décomposition orthogonale exponentielle.  _La Revenue Scientifique_ , 84:159–162, 1946.
  * Machalová et al. (2021) Jitka Machalová, Renáta Talská, Karel Hron, and Aleš Gába.  Compositional splines for representation of density functions.  _Computational Statistics_ , 36:1–34, 06 2021.  doi: 10.1007/s00180-020-01042-7.
  * Maier et al. (2022) Eva-Maria Maier, Almond Stöcker, Bernd Fitzenberger, and Sonja Greven.  Additive density-on-scalar regression in Bayes hilbert spaces with an application to gender economics, 2022.  arXiv:2110.11771.
  * Matabuena et al. (2021) Marcos Matabuena, Alexander Petersen, Juan C Vidal, and Francisco Gude.  Glucodensities: a new representation of glucose profiles using distributional data analysis.  _Statistical methods in medical research_ , 30(6):1445–1464, 2021.
  * Nerini and Ghattas (2007) David Nerini and Badih Ghattas.  Classifying densities using functional regression trees: Applications in oceanology.  _Computational Statistics & Data Analysis_, 51(10):4984–4993, 2007.  ISSN 0167-9473.  doi: https://doi.org/10.1016/j.csda.2006.09.028.  URL https://www.sciencedirect.com/science/article/pii/S0167947306003550.
  * Panaretos and Zemel (2019) Victor M. Panaretos and Yoav Zemel.  Statistical aspects of Wasserstein distances.  _Annual Review of Statistics and Its Application_ , 6(1):405–431, 2019.  doi: 10.1146/annurev-statistics-030718-104938.  URL https://doi.org/10.1146/annurev-statistics-030718-104938.
  * Petersen and Müller (2016) Alexander Petersen and Hans-Georg Müller.  Functional data analysis for density functions by transformation to a Hilbert space.  _Annals of Statistics_ , 44:183–218, 2016.
  * Qiu et al. (2022) Jiaming Qiu, Xiongtao Dai, and Zhengyuan Zhu.  Nonparametric estimation of repeated densities with heterogeneous sample sizes.  _Journal of the American Statistical Association_ , 0(0):1–13, 2022.  doi: 10.1080/01621459.2022.2104728.  URL https://doi.org/10.1080/01621459.2022.2104728.
  * Ramsay and Silverman (2005) James O. Ramsay and Bernhard W. Silverman.  _Functional Data Analysis_.  Springer New York, 2005.
  * Rasmussen and Williams (2005) C.E. Rasmussen and C.K.I. Williams.  _Gaussian Processes for Machine Learning_.  Adaptive Computation and Machine Learning series. MIT Press, 2005.  ISBN 9780262182539.
  * Scheipl et al. (2015) Fabian Scheipl, Ana-Maria Staicu, and Sonja Greven.  Functional additive mixed models.  _Journal of Computational and Graphical Statistics_ , 24(2):477–501, 2015.  doi: 10.1080/10618600.2014.901914.  PMID: 26347592.
  * Scimone et al. (2021) Riccardo Scimone, Alessandra Menafoglio, Laura Sangalli, and Piercesare Secchi.  A look at the spatio-temporal mortality patterns in Italy during the COVID-19 pandemic through the lens of mortality densities.  _Spatial Statistics_ , 49:100541, 09 2021.  doi: 10.1016/j.spasta.2021.100541.
  * Srivastava et al. (2007) Anuj Srivastava, Ian Jermyn, and Shantanu Joshi.  Riemannian analysis of probability density functions with applications in vision.  In _2007 IEEE Conference on Computer Vision and Pattern Recognition_ , pages 1–8, 2007.  doi: 10.1109/CVPR.2007.383188.
  * van den Boogaart et al. (2014) Karl Gerald van den Boogaart, Juan José Egozcue, and Vera Pawlowsky-Glahn.  Bayes hilbert spaces.  _Australian & New Zealand Journal of Statistics_, 56(2):171–194, 2014.  doi: https://doi.org/10.1111/anzs.12074.
  * Volkmann et al. (2023) Alexander Volkmann, Almond Stöcker, Fabian Scheipl, and Sonja Greven.  Multivariate functional additive mixed models.  _Statistical Modelling_ , 23(4):303–326, 2023\.  doi: 10.1177/1471082X211056158.
  * Wei and Tanner (1990) Greg C. G. Wei and Martin A. Tanner.  A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms.  _Journal of the American Statistical Association_ , 85:699–704, 1990.
  * Windmann and Kauermann (2019) M Windmann and G Kauermann.  Statistische Analyse der Nettomieten.  In _Mietspiegel für München 2019. Statistik, Dokumentation und Analysen_ , pages 19–70. Landeshauptstadt München, Sozialreferat, Amt für Wohnen und Migration, 2019.  Chapter 2.
  * Wood (2017) S.N. Wood.  _Generalized Additive Models: An Introduction with R, Second Edition_.  Chapman & Hall/CRC Texts in Statistical Science. CRC Press, 2017.  ISBN 9781498728348.  URL https://books.google.de/books?id=HL-PDwAAQBAJ.
  * Yao et al. (2005) Fang Yao, Hans-Georg Müller, and Jane-Ling Wang.  Functional linear regression analysis for longitudinal data.  _The Annals of Statistics_ , 33(6):2873–2903, 2005.

##  Appendix A Proofs and Computations

###  A.1 Proof of Lemma 2.2

###### Proof.

The clr transformation is well defined as

$\displaystyle\left|\int_{I}\log(f(x))dx\right|\leq|I|\left\|\log(f)\right\|_{\mathbb{L}^{2}}<\infty,$

due to the Cauchy-Schwarz inequality and $\log(f)\in\mathbb{L}^{2}$. We
further compute

$\displaystyle\operatorname{clr}([\alpha f])$ $\displaystyle=\log(\alpha
f)-\frac{1}{|I|}\int_{I}\log(\alpha f(x))dx$
$\displaystyle=\log(\alpha)+\log(f)-\frac{1}{|I|}\int_{I}\log(\alpha)+\log(f(x))dx$
$\displaystyle=\log(f)-\frac{1}{|I|}\int_{I}\log(f(x))dx=\operatorname{clr}([f])$

for all $\alpha\in\mathbb{R}$ and $f\in B$. Since $\operatorname{clr}^{-1}$ is
clearly well defined as well, we show that $\operatorname{clr}$ is bijective
via showing

$\displaystyle\operatorname{clr}(\operatorname{clr}^{-1}(g))=\operatorname{clr}([\exp(g)])=g-\frac{1}{|I|}\int_{I}g(x)dx=g$

for all $g\in\mathbb{L}_{0}^{2}$, and

$\displaystyle\operatorname{clr}^{-1}(\operatorname{clr}(f))=\left[\exp\left(\log(f)-\frac{1}{|I|}\log(f(x))dx\right)\right]=\left[f\exp\left(-\frac{1}{|I|}\log(f(x))dx\right)\right]=[f]$

for all $f\in B$. The clr transformation is also an isometry as for all
$g_{1},g_{2}\in\mathbb{L}_{0}^{2}$ holds

$\displaystyle\langle[\operatorname{clr}^{-1}(g_{1})],[\operatorname{clr}^{-1}(g_{2})]\rangle_{\mathcal{B}}=$
$\displaystyle\frac{1}{2|I|}\int_{I}\int_{I}\log\left(\frac{\exp(g_{1}(x))}{\exp(g_{1}(y))}\right)\log\left(\frac{\exp(g_{2}(x))}{\exp(g_{2}(y))}\right)dx\
dy$ $\displaystyle=$
$\displaystyle\frac{1}{2|I|}\int_{I}\int_{I}\left(g_{1}(x)-g_{1}(y)\right)\left(g_{2}(x)-g_{2}(y)\right)dx\
dy$ $\displaystyle=$
$\displaystyle\frac{1}{2|I|}\int_{I}\int_{I}\left(g_{1}(x)g_{2}(x)-g_{1}(y)g_{2}(x)-g_{1}(x)g_{2}(y)+g_{1}(y)g_{2}(y)\right)dx\
dy$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{I}g_{1}(x)g_{2}(x)dx-\frac{1}{2|I|}\int_{I}g_{1}(y)dy\int_{I}g_{2}(x)dx-\frac{1}{2|I|}\int_{I}g_{1}(x)dx\int_{I}g_{2}(y)dy$
$\displaystyle+\frac{1}{2}\int_{I}g_{1}(y)g_{2}(y)dy$ $\displaystyle=$
$\displaystyle\int_{I}g_{1}(x)g_{2}(x)dx=\langle
g_{1},g_{2}\rangle_{\mathbb{L}^{2}}$

since $\int_{I}g_{1}(x)dx=\int_{I}g_{2}(x)dx=0$. ∎

###  A.2 Proof of Corollary 2.4

###### Proof.

It is trivial to see that if $G=\sum_{k=1}^{N}\theta_{k}e_{k}$ with
$\boldsymbol{\theta}\overset{i.i.d.}{\sim}\mathcal{N}(\boldsymbol{\nu},\boldsymbol{\Sigma})$,
then $G\overset{i.i.d.}{\sim}GP(\mu,K)$ with $\mu=\sum_{k=1}^{N}\nu_{k}e_{k}$
and
$K(x_{1},x_{2})=\sum_{k=1}^{N}\sum_{l=1}^{N}e_{k}(x_{1})e_{l}(x_{2})\boldsymbol{\Sigma}_{kl}$
since

$\displaystyle\mathbb{E}(G(x))$
$\displaystyle=\mathbb{E}\left(\sum_{k=1}^{N}\theta_{k}e_{k}(x)\right)=\sum_{k=1}^{N}\nu_{k}e_{k}\quad\text{and}$
$\displaystyle\text{Cov}(G(x_{1}),G(x_{2}))$
$\displaystyle=\text{Cov}\left(\sum_{k=1}^{N}\theta_{k}e_{k}(x_{1}),\sum_{l=1}^{N}\theta_{l}e_{l}(x_{2})\right)=\sum_{k=1}^{N}\sum_{l=1}^{N}e_{k}(x_{1})e_{l}(x_{2})\text{Cov}(\theta_{k},\theta_{l})$
$\displaystyle=\sum_{k=1}^{N}\sum_{l=1}^{N}e_{k}(x_{1})e_{l}(x_{2})\boldsymbol{\Sigma}_{kl}.$
(11)

The other direction is an implication of the Karhunen-Loève decomposition. If
$\mathcal{H}$ is $N^{\prime}$-dimensional, $G\overset{i.i.d.}{\sim}GP(\mu,K)$
can be decomposed as $G=\mu+\sum_{k=1}^{N^{\prime}}Z_{k}\varphi_{k}$ with
$\varphi_{k}$, $k=1,\dots,{N^{\prime}}$ being the orthonormal eigenfunctions
and uncorrelated scores $Z_{k}$ with $\mathbb{E}(Z_{k})=0$ and
Var$(Z_{k})=\sigma_{k}^{2}$.

Choose $\nu_{k}=\langle\mu,e_{k}\rangle_{\mathbb{L}_{2}}$ to be the
orthonormal projection on $e_{k}$ and
$\boldsymbol{\Sigma}_{kl}=\sum_{j=1}^{N^{\prime}}\sigma_{j}^{2}\langle\varphi_{j},e_{k}\rangle_{\mathbb{L}_{2}}\langle\varphi_{j},e_{l}\rangle_{\mathbb{L}_{2}}$
for all $k,l=1,\dots,N$. Then we compute

$\displaystyle\sum_{k=1}^{N}\nu_{k}e_{k}$
$\displaystyle=\sum_{k=1}^{N}\langle\mu,e_{k}\rangle_{\mathbb{L}_{2}}e_{k}=\mu,$

since this gives the orthogonal projection of $\mu$ on
span$\\{e_{1},\dots,e_{N}\\}$ and
$\mu\in\mathcal{H}\subseteq\text{span}\\{e_{1},\dots,e_{N}\\}$ . We further
compute using the same identity

$\displaystyle\text{Cov}(G(x_{1}),G(x_{2}))$
$\displaystyle=\text{Cov}\left(\sum_{k=1}^{N^{\prime}}Z_{k}\varphi_{k}(x_{1}),\sum_{l=1}^{N^{\prime}}Z_{l}\varphi_{l}(x_{2})\right)=\sum_{k=1}^{N^{\prime}}\sum_{l=1}^{N^{\prime}}\varphi_{k}(x_{1})\varphi_{l}(x_{2})\text{Cov}(Z_{k},Z_{l})$
$\displaystyle=\sum_{k=1}^{N^{\prime}}\sigma_{k}^{2}\varphi_{k}(x_{1})\varphi_{k}(x_{2})$

which is identical to

$\displaystyle\sum_{k=1}^{N}\sum_{l=1}^{N}e_{k}(x_{1})e_{l}(x_{2})\boldsymbol{\Sigma}_{kl}$
$\displaystyle=\sum_{k=1}^{N}e_{k}(x_{1})\sum_{l=1}^{N}e_{l}(x_{2})\sum_{j=1}^{N^{\prime}}\sigma_{j}^{2}\langle\varphi_{j},e_{k}\rangle_{\mathbb{L}_{2}}\langle\varphi_{j},e_{l}\rangle_{\mathbb{L}_{2}}$
$\displaystyle=\sum_{j=1}^{N^{\prime}}\sigma_{j}^{2}\sum_{k=1}^{N}e_{k}(x_{1})\langle\varphi_{j},e_{k}\rangle_{\mathbb{L}_{2}}\sum_{l=1}^{N}e_{l}(x_{2})\langle\varphi_{j},e_{l}\rangle_{\mathbb{L}_{2}}$
$\displaystyle=\sum_{j=1}^{N^{\prime}}\sigma_{j}^{2}\varphi_{j}(x_{1})\varphi_{j}(x_{2}).$

To show the correspondence of the eigenvalue decompositions we need to show
that if $\boldsymbol{v}_{l}=(v_{l1},\dots,v_{lN})$ is an eigenvector of
$\boldsymbol{\Sigma}$ with corresponding eigenvalue $\sigma^{2}_{l}$ then
$\varphi_{l}=\sum_{m=1}^{N}v_{lm}e_{m}$ is an eigenfunction of $K$ with the
same eigenvalue $\sigma^{2}_{l}$. This is true since if we plug in the formula
obtained in (A.2) for $K$ we obtain

$\displaystyle\int_{I}K(x_{1},\cdot)\varphi_{l}(x_{1})dx_{1}$
$\displaystyle=\int_{I}K(x_{1},\cdot)\sum_{m=1}^{N}v_{lm}e_{m}(x_{1})dx_{1}=\sum_{k=1}^{N}e_{k}\int_{I}\left(\sum_{j=1}^{N}e_{j}(x_{1})\boldsymbol{\Sigma}_{kj}\right)\sum_{m=1}^{N}v_{lm}e_{m}(x_{1})dx_{1}$
$\displaystyle=\sum_{k=1}^{N}e_{k}\sum_{j=1}^{N}\boldsymbol{\Sigma}_{kj}\sum_{m=1}^{N}v_{lm}\int_{I}e_{j}(x_{1})e_{m}(x_{1})dx_{1}=\sum_{k=1}^{N}e_{k}\sum_{j=1}^{N}\boldsymbol{\Sigma}_{kj}v_{lj}$
$\displaystyle=\sum_{k=1}^{N}e_{k}\left(\boldsymbol{\Sigma}\boldsymbol{v}_{l}\right)_{k}=\sum_{k=1}^{N}e_{k}\left(\sigma^{2}_{l}\boldsymbol{v}_{l}\right)_{k}=\sigma^{2}_{l}\sum_{k=1}^{N}e_{k}v_{lk}=\sigma^{2}_{l}\varphi_{l}.$

###  A.3 Derivation of the likelihood

$\displaystyle L(\mu,K|\boldsymbol{x}_{i},\dots,\boldsymbol{x}_{n})$
$\displaystyle=\prod_{i=1}^{n}p(\boldsymbol{x}_{i}|\mu,K)=\prod_{i=1}^{n}\int_{\mathbb{R}^{N}}p(\boldsymbol{x}_{i}|\boldsymbol{\theta}_{i})p(\boldsymbol{\theta}_{i}|\mu,K)d\boldsymbol{\theta}_{i}$
$\displaystyle=\prod_{i=1}^{n}\int_{\mathbb{R}^{N}}\left(\prod_{j=1}^{m_{i}}p(x_{ij}|\boldsymbol{\theta}_{i})\right)p(\boldsymbol{\theta}_{i}|\mu,K)d\boldsymbol{\theta}_{i}$
$\displaystyle=\prod_{i=1}^{n}\int_{\mathbb{R}^{N}}\left(\prod_{j=1}^{m_{i}}\frac{\exp(\sum_{k=1}^{N}\theta_{ik}e_{k}(x_{ij}))}{\int_{I}\exp(\sum_{k=1}^{N}\theta_{ik}e_{k}(x))dx}\right)p(\boldsymbol{\theta}_{i}|\mu,K)d\boldsymbol{\theta}_{i}$
$\displaystyle=\prod_{i=1}^{n}\int_{\mathbb{R}^{N}}\frac{\exp(\sum_{j=1}^{m_{i}}\sum_{k=1}^{N}\theta_{ik}e_{k}(x_{ij}))p(\boldsymbol{\theta}_{i}|\mu,K)}{\left(\int_{I}\exp(\sum_{k=1}^{N}\theta_{ik}e_{k}(x))dx\right)^{m_{i}}}d\boldsymbol{\theta}_{i}$

∎

###  A.4 Proof of Lemma 2.5

###### Proof.

To show this statement we consider the three non constant additive parts of
the logarithm of the conditional distribution
$\log(p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}))$.

  * •

The first part
$\sum_{j=1}^{m_{i}}\left(\mu^{(h)}(x_{ij})+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x_{ij})\right)$
is linear in $\boldsymbol{z}_{i}$, which means there is a constant
$M_{1}\in\mathbb{R}$ such that it is $M_{1}\boldsymbol{z}_{i}$.

  * •

The second part is always negative, as

$\displaystyle-
m_{i}\log\left(\int_{I}\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)dx\right)\leq-
m_{i}\int_{I}\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)dx=0$

where the inequality is due to Jensen’s inequality and the integral is equal
to zero as all functions are in $\mathbb{L}^{2}_{0}$.

  * •

The third part consists of normal densities and is therefore quadratic in
$\boldsymbol{z}_{i}$. More precisely, we have

$\displaystyle\log(\prod_{k=1}^{N}p(z_{ik}|{\sigma_{k}^{2}}^{(h)}))=\sum_{k=1}^{N}\frac{-z_{ik}^{2}}{2{\sigma_{k}^{2}}^{(h)}}+\text{const.}\leq-\frac{1}{{2\sigma_{1}^{2}}^{(h)}}\|\boldsymbol{z}_{i}\|^{2}+\text{const.}$

Taking these three parts together this shows that there are constants
$M_{1},M_{2}\in\mathbb{R}$ and $M_{3}>0$ such that
$\log(p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}))\leq-
M_{3}\|\boldsymbol{z}_{i}\|^{2}+M_{1}\boldsymbol{z}_{i}+M_{2}$, which implies
that
$p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$
is a proper density. ∎

###  A.5 Example: The posterior mode will not necessarily be attained if the
prior is improper

###### Proof.

Consider the $1$-dimensional case $N=1$ with densities defined on the unit
interval $I=[0,1]$ and only one basis function for the clr transformed
densities given as $e_{1}=\mathds{1}_{[0,0.5]}-\mathds{1}_{]0.5,1]}$. This
means also the parameter space $\mathbb{R}$ in this case and assuming for this
parameter $z_{i}\in\mathbb{R}$ an improper $1$-dimensional normal distribution
is equivalent to assuming a flat prior.

If we further assume there is only one observation $x_{i1}=0.2$ and for the
prior mean holds $\boldsymbol{\nu}^{(h)}=0$ we compute

$\displaystyle
p(z_{i}|x_{i1},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$
$\displaystyle\propto\frac{\exp(z_{i}e_{1}(x_{i1}))}{\int_{[0,1]}\exp(z_{i}e_{1}(x))dx}=\frac{\exp(z_{i})}{\int_{[0,0.5]}\exp(z_{i})dx+\int_{[0.5,1]}\exp(-z_{i})dx}$
$\displaystyle=\frac{2\exp(z_{i})}{\exp(z_{i})+\exp(-z_{i})}=\frac{2}{1+\exp(-2z_{i})},$

which is monotonously increasing in $z_{i}\in\mathbb{R}$. Hence it does not
attain its maximum and therefore
$p(z_{i}|x_{i1},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)})$ does not
define a proper distribution in this case. ∎

###  A.6 Proof of Lemma 2.6

Like in the proof of Lemma 2.5 (Appendix A.4) we consider the three non
constant additive parts of the logarithm of the conditional distribution
$\log(p(\mathbf{z}_{i}|\boldsymbol{x}_{i},\boldsymbol{\nu}^{(h)},\boldsymbol{\Sigma}^{(h)}))$.

  * •

For the linear part
$\sum_{j=1}^{m_{i}}\left(\mu^{(h)}(x_{ij})+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x_{ij})\right)$
the gradient with respect to $\boldsymbol{z}_{i}$ is given as

$\displaystyle\frac{\partial}{\partial\boldsymbol{z}_{i}}\sum_{j=1}^{m_{i}}\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x_{ij})=\sum_{j=1}^{m_{i}}\sum_{k=1}^{N}\boldsymbol{v}^{(h)}_{k}e_{k}(x_{ij})=\sum_{k=1}^{N}\boldsymbol{v}^{(h)}_{k}\sum_{j=1}^{m_{i}}e_{k}(x_{ij})$

  * •

For the second part we compute the gradient as

$\displaystyle\phantom{=}\frac{\partial}{\partial\boldsymbol{z}_{i}}-m_{i}\log\left(\int_{I}\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)dx\right)$
$\displaystyle=\frac{-m_{i}\int_{I}\frac{\partial}{\partial\boldsymbol{z}_{i}}\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)dx}{\int_{I}\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)dx}$
$\displaystyle=\frac{-m_{i}\int_{I}\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)\sum_{k=1}^{N}\boldsymbol{v}^{(h)}_{k}e_{k}(x)dx}{\int_{I}\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)dx}$
$\displaystyle=-m_{i}\int_{I}\frac{\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)}{\int_{I}\exp\left(\mu^{(h)}(x)+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}(x)\right)dx}\sum_{k=1}^{N}\boldsymbol{v}^{(h)}_{k}e_{k}(x)dx$
$\displaystyle=-m_{i}\int_{I}\operatorname{clr}^{-1}\left(\mu^{(h)}+\sum_{k=1}^{N}\boldsymbol{z}_{i}^{T}\boldsymbol{v}^{(h)}_{k}e_{k}\right)(x)\sum_{k=1}^{N}\boldsymbol{v}^{(h)}_{k}e_{k}(x)dx$
$\displaystyle=-m_{i}\sum_{k=1}^{N}\boldsymbol{v}^{(h)}_{k}\int_{I}f_{\mathbf{z}_{i}}(x)e_{k}(x)dx$

where, in the first equation, we can interchange differentiation and
integration applying the Leibniz rule, since we have assumed that all clr
transformed densities are bounded.

  * •

The third part is the sum of logarithms of normal densities. Therefore, for
all $l=1,\dots,N$ we compute the partial derivative with respect to $z_{il}$
as

$\displaystyle\frac{\partial}{\partial
z_{il}}\log(\prod_{k=1}^{N}p(z_{ik}|{\sigma_{k}^{2}}^{(h)}))=\frac{\partial}{\partial
z_{il}}\frac{-z_{il}^{2}}{2{\sigma_{l}^{2}}^{(h)}}=\frac{-z_{il}}{{\sigma_{l}^{2}}^{(h)}}.$

Adding these three parts together gives the gradient of the logarithm of the
conditional density of the scores.

##  Appendix B Additional plots for the applications in Section 3

###  B.1 Temperature data

For the temperature data, the latent density model is estimated using kernel
density estimates with a Gaussian kernel and bandwidth $=1.5$ as initial
estimates. The number of Monte Carlo samples in the E step 2.3.1 is chosen to
be $r=50h$, where $h$ is the iteration index, i.e., the number of samples
increases over the iterations. The tuning parameter $\lambda$ for the proposal
density is set to the default value $1$. In each iteration, the dimension
reduction is at most $0.001$, i.e. we keep as many principal components as
necessary to explain at least 99.99 % of the variance.

Figure 6: Kernel density estimates (Gaussian kernel, bandwidth $=1.5$) and
predicted latent densities estimated from the latent density model (5) for the
daily maximum summer temperatures in Berlin Tempelhof from 1951 to 2022. The
kernel densities are also used as initial estimates in the latent density
model.

In Figure 6 we show the kernel density estimates and the predicted latent
densities obtained from the latent density model. For both, the trend towards
higher temperatures is evident. Note that the variance of the predicted latent
densities is smaller than the variance of the kernel density estimates, since
the latent density model effectively splits the total variance into the
variance due to the underlying stochastic process for the latent densities and
the variance due to sampling from them. Correspondingly, percentages variance
explained for the PCs depicted in Figures 1 and 7, respectively, are relative
to different total variances.

Figure 7: PCA based on the clr transformations of the kernel density estimates
(Gaussian kernel, bandwidth $=1.5$) for daily maximum temperature (°C). Top:
Effect of adding/subtracting $\hat{\sigma}_{k}\hat{\varphi}_{k}$ to the clr
transformed mean density $\mu$, where $\hat{\varphi}_{k}$ is the $k$th
principal component, with corresponding eigenvalue $\hat{\sigma}_{k}^{2}$,
$k=1,2,3,4$. Middle: Effect on the density level, i.e.
$\operatorname{clr}^{-1}$ transformations of the functions in the top row.
Bottom: Temporal trend of the corresponding predicted scores per year, with
scatterplot smoother and pointwise confidence bands overlaid.

###  B.2 Rental prices

For the rent index data, the latent density model is estimated using kernel
density estimates with a Gaussian kernel and bandwidth $=2$ as initial
estimates. The number of Monte Carlo samples in the E step 2.3.1 is chosen to
be $r=100h$, where $h$ is the iteration index, i.e., the number of samples
increases over the iterations. The parameter $\lambda$ for the proposal
density is set to $2$. In each iteration, the dimension reduction is at most
$0.0005$, i.e. we keep as many principal components as necessary to explain at
least 99.995 % of the variance.

Figure 8: Histograms of rent per square meter for each district, with overlaid
kernel density estimates (Gaussian kernel, bandwidth $=2$) in blue and
predicted latent densities based on our proposed approach in red.