# Local moment matching with Erlang mixtures under automatic roughness
penalization

Oskar Laverny E-mail address<EMAIL_ADDRESS>Aix Marseille Univ,
Inserm, IRD, SESSTIM, Sciences Economiques & Sociales de la Santé & Traitement
de l’Information Médicale, ISSPAM, Marseille, France.  Philippe Lambert
Institut de Statistique, Biostatistique et Sciences Actuarielles (ISBA),
Université catholique de Louvain, Voie du Roman Pays 20, B-1348 Louvain-la-
Neuve, Belgium  Institut de Mathématique, Université de Liège, Allée de la
Découverte 12, B-4000 Liège, Belgium.

###### Abstract

We consider the class of Erlang mixtures for the task of density estimation on
the positive real line when the only available information is given as local
moments, a histogram with potentially higher order moments in some bins. By
construction, the obtained moment problem is ill-posed and requires
regularization. Several penalties can be used for such a task, such as a lasso
penalty for sparsity of the representation, but we focus here on a simplified
roughness penalty from the P-splines literature. We show that the
corresponding hyperparameter can be selected without cross-validation through
the computation of the so-called effective dimension of the estimator, which
makes the estimator practical and adapted to these summarized information
settings. The flexibility of the local moments representations allows
interesting additions such as the enforcement of Value-at-Risk and Tail Value-
at-Risk constraints on the resulting estimator, making the procedure suitable
for the estimation of heavy-tailed densities.

Keywords: Nonparametric density estimation, grouped data, sample moments,
erlang mixtures

##  1 Introduction

Let $X_{1},...,X_{n}$ be a $n$-sample of a positive real random variable $X$,
and let $\mathcal{B}=\left(B_{1},...,B_{J}\right)$ be a finite partition of
$\mathbb{R}_{+}$, where $B_{j}=[b_{j-1},b_{j}[$ are called _bins_. Let
furthermore $\bm{k}\in\mathbb{N}^{J}$ be a vector of integers. We consider
that the only information we have about the sample is the set of empirical
moments given by:

$\displaystyle\hat{\bm{\pi}}$
$\displaystyle=\left\\{\hat{\pi}_{j}=\frac{1}{N}\sum_{i=1}^{N}\mathbb{1}_{X_{i}\in
B_{j}}:j\in 1,...,J\right\\}\text{ and}$ (1) $\displaystyle\hat{\bm{\mu}}$
$\displaystyle=\left\\{\hat{\mu}_{j,k}=\frac{1}{N}\sum_{i=1}^{N}X_{i}^{k}\mathbb{1}_{X_{i}\in
B_{j}}:j\in 1,...J,k\in 1,...,k_{j}\right\\}.$

The zeroth moments are denoted separately strictly for ease of notations. In
an insurance or reinsurance context, such binned information may arise from
attritionnal versus large cutoff in a loss random variable, and the lack of
detailed information (only a few moments are available) may typically result
from confidentiality concerns.

Indeed, the General Data Protection Regulation (GDPR) has raised important
confidentiality concerns in the insurance industry. The regulation imposes
strict rules on the collection, use, and storage of personal data, including
data related to health, which is often used in underwriting and claims
management. The GDPR’s requirements for confidentiality in insurance pose
significant challenges for insurers, who must balance the need to collect and
use personal data with the need to protect individuals’ privacy. Several
studies have examined the impact of this regulation. While providing an
implementation guide, [23] notes that the process could be used by the
industry as an opportunity for fine-tuning the data management. On the other
hand, [7] investigates the way insurer respond to the the right to explanation
and the need for transparency in premium computations.

Because of these concerns, it is not uncommon for losses data to be
transmitted in a summarized state alike Equation 1, with rather small integers
$J,k_{1},...k_{J}$, to increase the protection of personal data by not
transmitting the full losses and/or premium information. A representative
example of such a problem drawn from a LogNormal$(0,0.5)$111We denote all
probability distributions by $\texttt{Name}(\text{Parameters})$. Their
definitions are all detailed in Appendix A. is given in Table 1.

Table 1: Example of a local moment problem drawn from a LogNormal$(0,0.5)$ sample. In this example, the three first lines gives estimated boxed moments. The fourth line is a little more involved: it prescribes a Value-At-Risk $b_{J-1}$ and a Tail-Value-at-Risk $\hat{\mu}_{J-1,1}$ at the quantile level $1-\hat{\pi}_{4}$. This data structure is classical when modeling insurance losses, which are usually divided into _attritional_ and _large_ losses. $j$ | $[b_{j-1},b_{j}[$ | $N_{j}$ | $k_{j}$ | $\pi_{j}$ |  $\hat{\mu}_{j,1}$ |  $\hat{\mu}_{j,2}$ |  $\hat{\mu}_{j,3}$ |  $\hat{\mu}_{j,4}$
---|---|---|---|---|---|---|---|---
1 | $[0.000,0.948[$ | 375 | 4 | 0.500 |  0.332 |  0.235 |  0.175 |  0.136
2 | $[0.948,1.885[$ | 300 | 4 | 0.400 |  0.526 |  0.719 |  1.017 |  1.488
3 | $[1.885,3.332[$ | 67 | 4 | 0.089 |  0.206 |  0.485 |  1.167 |  2.874
4 | $[3.332,\infty[$ | 8 | 1 | 0.011 |  0.048 |  |  |

The unavailability of higher order moments on the last line of Table 1
represents the fact that these moments might not be finite and are therefore
usually not available. This flexibility of the representation allows handling
heavy-tailed random variables through VaR/TVaR insights (or even only VaR). In
what follows, we only assume that the zeroth moments in each box are entirely
part of the available information (that is $k_{j}\geq 0$ for all $j\in J$).

From this restricted information, our goal is to estimate the density for the
positive (and potentially heavy-tailed) random variable $X$. If the previous
literature on such binned moment problems (see, e.g., [17] or the related [4])
dealt mainly with finitely supported random variables, we are here
particularly interested in the estimation of the tail of distribution. Our
proposition is to tackle this problem using the class of Erlang mixtures.

Erlang Mixtures have been widely used these past 30 years in various fields,
including telecommunications, manufacturing, reliability analysis and
insurance, in which they are used to model the distribution of losses, claim
amounts, and other risk-related variables. Their flexibility allows them to
capture a wide range of distributions, including skewed and heavy-tailed
distributions that are commonly found in insurance data. This has led to their
adoption in various areas of insurance, such as pricing, reserving, and risk
management. See, e.g., [22, 21] for pioneering work on their estimation
through Expectation-Maximization (EM) algorithms, in univariate or
multivariate settings respectively. The multivariate case is commonly used to
model the joint distribution of multiple risk-related variables, including the
dependence structure, such as claim amounts, frequency, and severity among
several lines of business. See [12] for applications to ruin problems (and
references therein), or [26] for application to crop insurance pricing. More
recently, [15, 14] uses multivariate Erlang mixtures to fit censored and
truncated multivariate data with a roughness penalty, again through a
(modified) EM algorithm. This literature shows that Erlang mixtures are able
to capture complex tail behaviors and that the model can provide valuable
insights on the underlying risk factors.

Overall, univariate and multivariate Erlang mixtures offer a flexible and
powerful framework for analyzing insurance data, which makes them a valuable
tool for insurance practitioners and researchers alike. The denseness results
from [29] and subsequent multivariate extensions are probably at the root of
the popularity of this model. We are particularly interested in the work in
[5], which solve standard moment problems on the real line with Erlang
mixtures, again with risk management applications, giving an alternative to
the classical estimator from [16]. More precisely, [16] and [6] solve our
problem for a unique, fully covering bin (that is, $J=1$ and
$B_{1}=\mathbb{R}_{+}$), corresponding to a standard moment problem on the
positive real line.

If our estimation problem is a little more involved than a standard moment
problem, it is easy to see that such an inverse problem will require
regularization. Therefore, we propose to leverage a penalization of the
roughness of the density, as was already discussed in [14]. However, two
problems arise:

  * •

The direct penalization of the continuous roughness of the density proposed in
[14] comes with a large computational cost.

  * •

In our low information settings, we do not have a full dataset and the
calibration of the penalty hyperparameter is not possible through cross-
validation (like in classical EM schemes) since holding out testing data would
not be acceptable.

On the other hand, the penalization of roughness of a density is largely dealt
with in the literature around approximation by splines, notably yielding the
concept of P-splines.

P-splines, or penalized B-splines, are a popular and flexible tool for fitting
smooth curves and surfaces to data. They were introduced by Eilers and Marx
[8] as a way to overcome some limitations of traditional B-spline methods,
such as the difficulty of selecting knots and the optimal degree of smoothing.
P-splines use a penalty on the $r^{\text{th}}$-order differences of the
B-spline coefficients to achieve smoothness, and this penalty can be easily
controlled through a tuning parameter. See [9] for extensive literature review
on P-splines.

Through a Bayesian perspective, we leverage the P-spline literature to provide
in this paper a solution to the two problems encountered with the roughness
penalization of Erlang mixtures. The discussion will be organized as follows:
Section 2 derives and motivates the proposed approximation scheme, Section 3
gives a few numerical illustrations of the method, including convergence study
on simulated datasets, and Section 4 concludes.

##  2 Proposed estimation scheme

Recall that we only observe $\hat{\bm{\pi}}$ and $\hat{\bm{\mu}}$ from
Equation 1. We also consider that we know $N$, the number of observations that
were used to compute the histogram $\hat{\bm{\pi}}$ and the boxed moments
$\hat{\bm{\mu}}$. Let’s assume without loss of generality that $X$ has
theoretical counterpart moments:

$\displaystyle\bm{\pi}$ $\displaystyle=\left\\{\pi_{j}=\mathbb{P}\left(X\in
B_{j}\right):j\in 1,...,J\right\\}\text{ and}$ (2) $\displaystyle\bm{\mu}$
$\displaystyle=\left\\{\mu_{j,k}=\mathbb{E}\left(X^{k}\mathbb{1}_{X\in
B_{j}}\right):j\in 1,...J,k\in 1,...,k_{j}\right\\}.$

We consider the underlying sample $X_{1},...X_{n}$ to be independently and
identically distributed. Under this assumption, we have that
$N\hat{\bm{\pi}}\sim\texttt{Multinomial}(\bm{\pi},N)$ (see Appendix A for
distribution definitions). We then approximate, conditionally on
$\hat{\bm{\pi}}$, the distribution of the vector of moments $\hat{\bm{\mu}}$
through the law of large numbers as
$\hat{\bm{\mu}}\sim\texttt{Normal}(\bm{\mu},\bm{\Sigma}/N)$, where the
variance-covariance matrix $\bm{\Sigma}$ is defined by:

$\Sigma_{(j,k),(i,m)}=\mu_{j,k+m}\mathbb{1}_{j=i}-\mu_{j,k}\mu_{i,m}.$

We generally denote by $\ell$ log-density, neglecting additive constants, of
random objects. In particular, the loglikelihood of our model can be expressed
(neglecting additive constants) as

$\ell(\bm{\pi},\bm{\mu},\bm{\Sigma})=\left\\{\hat{\bm{\pi}}^{\prime}\log(\bm{\pi})\right\\}+\left\\{-\frac{1}{2}\log|\bm{\Sigma}|-\frac{1}{2}\lVert\bm{\mu}-\hat{\bm{\mu}}\rVert_{\Sigma}^{2}\right\\},$
(3)

where brackets delimit the respective contribution of the multinomial
distribution and the conditional multivariate normal. Here, $\log(\bm{\pi})$
is intended componentwise, $\lvert\cdot\rvert$ denotes the determinant of a
matrix and
$\lVert\bm{x}\rVert_{\Sigma}^{2}=\bm{x}^{\prime}\bm{\Sigma}^{-1}\bm{x}$
denotes the squared Mahalanobis distance. It is clear from this expression
that, under this Gaussian approximation for $\hat{\bm{\mu}}$, the triplet
$(\bm{\pi},\bm{\mu},\bm{\Sigma})$ is sufficient to evaluate the loglikelihood
of our data, whatever the underlying parametric or non-parametric model we
choose.

Due to the positive and potentially heavy-tailed nature of the random
variable, we propose to model it through the class of Erlang Mixtures. Recall
that the $\texttt{MixedGamma}(\nu,\theta)$ distribution, where
$\theta\in\mathbb{R}_{+}$ is the scale and $\nu\in\mathcal{P}(\mathbb{R}_{+})$
is the prior distribution of shapes, has density

$f_{\nu,s}(x)=\int_{\mathbb{R}_{+}}\frac{x^{\alpha-1}}{\theta^{\alpha}\Gamma(\alpha)}e^{-\frac{x}{\theta}}\nu(d\alpha).$

The restriction of the support of $\nu$ to $\mathbb{N}$, setting e.g.,
$\bm{\omega}=\left\\{\omega_{i}=\nu(i),i\in\mathbb{N}\right\\}$, give rise to
the $\texttt{MixedErlang}(\bm{\omega},\theta)$ distribution. The number of
non-zero $\omega_{i}$’s is usually finite, and in practice we denote
$n=\max_{\omega_{i}>0}i$. One of the reasons that drove the literature
interests upon this class is the denseness results by Tijms [29, p. 163].
Tijms show that this restricted class can approximate arbitrary distributions
on the positive real line. More precisely, if we let $F$ denote a distribution
function on the real line, setting weights as
$\omega_{i}=F(i\theta)-F((i-1)\theta)$, we have that the corresponding
distribution function converges to $F$ pointwise as $\theta\to 0$, that is

$\forall
x\in\mathbb{R}_{+},\;F_{\bm{\omega},\theta}(x)\xrightarrow[]{\theta\to
0}F(x).$

A multivariate extension of this result can be found in [21, Theorem 2.1]. The
result can be easily proved using the convergence of moment generating
functions. The result is constructive, and gives, for any scale $\theta$, a
set of weights $\omega_{i}=F(i\theta)-F((i-1)\theta)$ that, if we could
estimate them, would lead to a nice starting point for an iterative algorithm.
Note that using them directly as an estimator is not recommended, as, as was
already noted by [22], this clearly overfits a given dataset. Moreover, in our
case, we do not have enough data do estimate this discretization of $F$.

Note that, for $X\sim\texttt{MixedGamma}(\nu,\theta)$, $\bm{\pi},\bm{\mu}$ and
$\bm{\Sigma}$ are easily computed. Indeed, it can be checked that binned
moments are given by

$\mathbb{E}\left(X^{k}\mathbb{1}_{X\in[a,b[}\right)=\theta^{k}\int\frac{\Gamma(\alpha+k)}{\Gamma(\alpha)}\left(\gamma\left(\frac{b}{\theta},\alpha+k\right)-\gamma\left(\frac{a}{\theta},\alpha+k\right)\right)\nu(d\alpha),$
(4)

where $\gamma(x,\alpha)=\Gamma(\alpha)^{-1}\int_{0}^{x}t^{\alpha-1}e^{-t}dt$
is the lower incomplete gamma function. The expression in Equation 4 can be
efficiently computed recursively on $k$, and, for
$X\sim\texttt{MixedErlang}(\bm{\omega},\theta)$, recursively on atoms of
$\nu$. From Equation 2 we see that these integrals are enough to compute
$\bm{\pi},\bm{\mu}$ and $\bm{\Sigma}$, and from Equation 3 our loglikelihood.
With these implicit definitions, the loglikelihood of a
$\texttt{MixedErlang}(\bm{\omega},\theta)$ for our sample is denoted
$\ell(\bm{\omega},\theta)=\ell(\bm{\pi},\bm{\mu},\bm{\Sigma})$.

Due to overfitting concerns on one side, and due to spiky behavior of Erlang
mixtures, we need to include a regularization mechanism in the fitting
process. Smoothness in a density estimate is often obtained (see e.g., [14] in
the context of Erlang mixtures) by penalizing large norms of the density
derivatives through the addition of a penalty

$\widetilde{\mathrm{Pen}_{r}}(\lambda,f)=\frac{\lambda}{2}\int
f^{(r)}(x)^{2}dx,$

for a given positive integer $r$.

For $f\sim\texttt{MixedErlang}(\bm{\omega},\theta)$,
$\widetilde{\mathrm{Pen}_{r}}(\lambda,f)=\frac{\lambda}{2}\theta^{-(2r+1)}\bm{\omega}^{\prime}\tilde{\bm{P}_{r}}\bm{\omega}$,
where $\tilde{\bm{P}_{r}}$ is a fixed positive semidefinite dense matrix. The
elements of the matrix $\tilde{\bm{P}_{r}}$ have closed form expressions given
by

$\tilde{P}_{r,i,j}=\sum_{k=0}^{r}\sum_{\ell=0}^{r}c_{r,k}c_{r,\ell}\frac{\Gamma(i+j-k-\ell-1)}{\Gamma(i-k)\Gamma(j-\ell)}2^{-(i+j-k-\ell-1)}\mathbb{1}_{i-k>0,j-l>0},$

where $c_{r,k}$ are finite difference coefficients of order $r$, defined
recursively by $c_{0,0}=c_{1,0}=1$ and
$c_{r,k}=\left(c_{r-1,k}-c_{r-1,k-1}\right)\mathbb{1}_{k\in\left\\{0,...r\right\\}}$.
While looking for an optimal hyperparameter $\lambda$, the term
$\theta^{-(2r-1)}$ could be simply incorporated into $\lambda$ and removed
from the equation. However, the matrix $\tilde{\bm{P}_{r}}$ is dense, which
can be computationally problematic, and in our low information settings we do
not have the possibility to run cross-validation strategies to find the right
$\lambda$.

We solve these issues by considering instead another penalty, drawn from the
P-splines literature. Instead of penalizing the continuous roughness of $f$,
we penalize the discrete roughness of the sequence of modes of the underlying
Erlang densities. Neglecting a multiplicative $\theta^{-1}$ factor, the mode
of the $(i+1)^{\text{th}}$ Erlang density in the mixture is positioned at

$(x_{i},y_{i})=\left(\theta i,\frac{i^{i}e^{-i}}{i!}\right).$ (5)

Note that $\bm{x}=\\{x_{1},x_{2},...\\}$ are regularly distributed on the
positive real line. Therefore, we monitor the regularity of the density
estimate by the regularity of the sequence of weighted modes
$\omega_{1}y_{1},\omega_{2}y_{2},...$ Let us denote by $\bm{D}_{r}$ the
$r^{\text{th}}$ difference matrix, of size $(n-r,n)$, defined by

$D_{r,k,l}=c_{r,l-k}y_{l}\mathbb{1}_{l-k\leq r},$

such that $\bm{D}_{r}\omega$ are the finite differences of the sequence of
modes. The matrix $\bm{D}_{r}^{\prime}\bm{D}_{r}$ is sparse (only its $(2r-1)$
central diagonals are non-zero) which will be computationally efficient. The
corresponding penalty writes

$\frac{\lambda}{2}\lVert\bm{D}_{r}\bm{\omega}\rVert_{2}^{2}.$

We can interpret this penalization as a prior on $r^{\text{th}}$ order
differences of coefficients
$\bm{D}_{r}\bm{\omega}|\lambda\sim\texttt{Normal}(\bm{0},\lambda^{-1}\bm{I})$,
where $\bm{I}$ represent the identity matrix. Assuming furthermore that we
assign a $\texttt{Gamma}(a_{\lambda},b_{\lambda})$ prior on $\lambda$, with
hyperparameters $a_{\lambda},b_{\lambda}$ such that the prior is uninformative
(high variance), the final penalization is

$\mathrm{Pen}_{r}(\lambda,\bm{\omega})=\frac{1}{2}\left\\{-(n-r)\log(\lambda)+\lambda\lVert\bm{D}_{r}\bm{\omega}\rVert_{2}^{2}\right\\}+\left\\{\lambda
b_{\lambda}^{-1}-(a_{\lambda}-1)\log(\lambda)\right\\}.$

Thus, we set, still neglecting additive constants,

$\displaystyle\ell(\bm{\omega},\theta,\lambda)=$
$\displaystyle\ell(\bm{\omega},\theta)-\mathrm{Pen}_{r}(\lambda,\bm{\omega})$
$\displaystyle=$
$\displaystyle\left\\{\hat{\bm{\pi}}^{\prime}\log(\bm{\pi})\right\\}-\frac{1}{2}\left\\{\log|\bm{\Sigma}|+\lVert\bm{\mu}-\hat{\bm{\mu}}\rVert_{\Sigma}^{2}\right\\}$
$\displaystyle-\frac{1}{2}\left\\{-(n-r)\log(\lambda)+\lambda\lVert\bm{D}_{r}\bm{\omega}\rVert_{2}^{2}\right\\}$
$\displaystyle-\left\\{\lambda
b_{\lambda}^{-1}-(a_{\lambda}-1)\log(\lambda)\right\\},$

where the link between the triplet
$\left(\bm{\pi},\bm{\mu},\bm{\Sigma}\right)$ and parameters
$\left(\bm{\omega},\theta\right)$ is here implicitly assumed as given by
Equations 2 and 4. Note that the added components do not depend on
$(\omega,\theta)$ and are thus compatible with other approaches (like cross-
validation) that uses the direct penalty
$\frac{\lambda}{2}\lVert\bm{D}_{r}\bm{\omega}\rVert_{2}^{2}.$ On the other
hand, they complete the loglikelihood which allows a Bayesian interpretation.

The marginal loglikelihood for $\lambda$ can be obtained using the following
identity:

$\ell(\lambda)=\ell(\bm{\omega},\theta,\lambda)-\ell(\bm{\omega},\theta|\lambda),$
(6)

where $\ell(\bm{\omega},\theta|\lambda)$ denotes the conditional log posterior
of the Erlang mixtures parameters for a given penalty parameter $\lambda$.

Consider the following notations for the Hessian matrices involved in our
problem, which are linked together by a linearity in $\lambda$:

$\displaystyle\bm{H}(\bm{\omega},\theta)$
$\displaystyle=-\frac{\partial^{2}}{\partial^{2}(\bm{\omega},\theta)}\ell(\bm{\omega},\theta)$
$\displaystyle\bm{P}_{r}$
$\displaystyle=\frac{\partial}{\partial\lambda}\frac{\partial^{2}}{\partial^{2}(\bm{\omega},\theta)}\mathrm{Pen}_{r}(\lambda,\bm{\omega})=\begin{pmatrix}\bm{D}_{r}^{\prime}\bm{D}_{r}&0\\\
0&0\end{pmatrix}$ $\displaystyle\bm{H}(\bm{\omega},\theta,\lambda)$
$\displaystyle=-\frac{\partial^{2}}{\partial^{2}(\bm{\omega},\theta)}\ell(\bm{\omega},\theta,\lambda)=\bm{H}(\bm{\omega},\theta)+\lambda\bm{P}_{r}.$

The conditional log posterior $\ell(\bm{\omega},\theta|\lambda)$ can then be
approximated using a Laplace approximation. The Laplace approximation (see
e.g., [28]) is a computationally efficient method to approximate the posterior
distribution of penalized parameters. It has been applied to a wide range of
Bayesian models, including generalized linear models [28], mixed effects
models [3], and latent variable models [27]. Several studies have shown that
Laplace approximation can provide accurate and efficient estimates for
P-spline models, even for large datasets and complex models with high-
dimensional spline bases [10, 30, 18, 13]. Overall, Laplace approximation
remains a popular method for Bayesian inference in the P-splines literature.

Here, we approximate the posterior distribution around the posterior mode
$(\hat{\bm{\omega}},\hat{\theta})$, that is (neglecting again additive
constants):

$\ell(\bm{\omega},\theta|\lambda)\approx\frac{1}{2}\log\lvert\bm{H}(\hat{\bm{\omega}},\hat{\theta},\lambda)\rvert-\frac{1}{2}\lVert(\bm{\omega},\theta)-(\hat{\bm{\omega}},\hat{\theta})\rVert_{\bm{H}(\hat{\bm{\omega}},\hat{\theta},\lambda)^{-1}}^{2}.$
(7)

This approximation become particularly useful if we note that Equation 6 is
true for any $\bm{\omega},\theta$. Thus, by evaluating Equation 7 near the
posterior mode $(\hat{\bm{\omega}},\hat{\theta})$, with substitution in
Equation 6, we obtain:

$\ell(\bm{\omega},\theta|\lambda)\approx\frac{1}{2}\log\lvert\bm{H}(\bm{\omega},\theta,\lambda)\rvert.$

Taking the derivative with respect to $\lambda$, using the fact that for a
positive-definite matrix $\bm{A}(x)$, $\frac{\partial}{\partial
x}\log\lvert\bm{A}(x)\rvert=\mathrm{Tr}\left\\{A(x)^{-1}\frac{\partial}{\partial
x}\bm{A}(x)\right\\}$, one has:

$\displaystyle\frac{\partial}{\partial\lambda}\ell(\bm{\omega},\theta|\lambda)$
$\displaystyle\approx\frac{\partial}{\partial\lambda}\frac{1}{2}\log\lvert\bm{H}(\bm{\omega},\theta,\lambda)\rvert$
$\displaystyle=\frac{1}{2}\mathrm{Tr}\left\\{\bm{H}(\bm{\omega},\theta,\lambda)^{-1}\frac{\partial}{\partial\lambda}\bm{H}(\bm{\omega},\theta,\lambda)\right\\}$
$\displaystyle=\frac{1}{2}\mathrm{Tr}\left\\{\bm{H}(\bm{\omega},\theta,\lambda)^{-1}\bm{P}_{r}\right\\}$
$\displaystyle=\frac{1}{2}\mathrm{Tr}\left\\{\left(\bm{H}(\bm{\omega},\theta)+\lambda\bm{P}_{r}\right)^{-1}\bm{P}_{r}\right\\},$

By a recursive application of the Woodbury matrix identity, for two matrices
$\bm{A}$ and $\bm{B}$ we have that
$(\bm{A}-\bm{B})^{-1}=\sum_{k=0}^{\infty}(\bm{A}^{-1}\bm{B})^{k}\bm{A}^{-1},$
assuming the sum converges, that is assuming the spectral radius of
$\bm{A}^{-1}\bm{B}$ is less than one. In our contexts, this yields :

$\left(\bm{H}(\bm{\omega},\theta)-(-\lambda\bm{P}_{r})\right)^{-1}\bm{P}_{r}=-\sum_{k=0}^{\infty}\lambda^{k}\left(-\bm{H}(\bm{\omega},\theta)^{-1}\bm{P}_{r}\right)^{k+1}.$

Assuming that $\eta_{1},...\eta_{n+1}$ are eigenvalues of
$-\bm{H}(\bm{\omega},\theta)^{-1}\bm{P}_{r}$, we therefore have:

$\displaystyle\mathrm{Tr}\left\\{\left(\bm{H}(\bm{\omega},\theta)+\lambda\bm{P}_{r}\right)^{-1}\bm{P}_{r}\right\\}$
$\displaystyle=-\sum_{k=0}^{\infty}\lambda^{k}\mathrm{Tr}\left\\{\left(-\bm{H}(\bm{\omega},\theta)^{-1}\bm{P}_{r}\right)^{k+1}\right\\}$
$\displaystyle=-\sum_{k=0}^{\infty}\lambda^{k}\sum_{i=1}^{n+1}\eta_{i}^{k+1}$
$\displaystyle=-\sum_{i=1}^{n+1}\frac{\eta_{i}}{1-\lambda\eta_{i}}.$

Finally,

$\displaystyle\frac{\partial}{\partial\lambda}\ell(\lambda)$
$\displaystyle=\frac{\partial}{\partial\lambda}\ell(\bm{\omega},\theta,\lambda)-\frac{\partial}{\partial\lambda}\ell(\bm{\omega},\theta|\lambda)$
$\displaystyle=-\frac{1}{2}\left\\{-\frac{n-r}{\lambda}+\lVert\bm{D}_{r}^{\prime}\bm{\omega}\rVert_{2}^{2}\right\\}+\left\\{\frac{a_{\lambda}-1}{\lambda}-b_{\lambda}^{-1}\right\\}-\frac{1}{2}\mathrm{Tr}\left\\{\left(\bm{H}(\bm{\omega},\theta)+\lambda\bm{P}_{r}\right)^{-1}\bm{P}_{r}\right\\}$
$\displaystyle=\frac{-1}{2}\left\\{\lVert\bm{D}_{r}^{\prime}\bm{\omega}\rVert_{2}^{2}-\frac{(n-r)+2a_{\lambda}-2}{\lambda}+2b_{\lambda}^{-1}-\sum_{i=1}^{n+1}\frac{\eta_{i}}{1-\lambda\eta_{i}}\right\\},$
(8)

which, since $\bm{\eta}$ and
$\lVert\bm{D}_{r}^{\prime}\bm{\omega}\rVert_{2}^{2}$ are constants w.r.t.
$\lambda$, gives us a simple non-linear (polynomial of degree $n+1$) equation
to find critical points of $\ell(\lambda)$, which can be solved with, e.g., a
fixed-point method of any gradient ascent method, reaching the closest local
maxima.

Thus, Algorithm 1 uses an iterative algorithm, by recursively:

  1. 1.

Update $(\bm{\omega},\theta)=\mathrm{arg}\min\ell(\bm{\omega},\theta,\lambda)$
at the current value of $\lambda$,

  2. 2.

Update $\lambda=\mathrm{arg}\min\ell(\lambda)$ at the current value of
$(\bm{\omega},\theta)$ by setting Section 2 to $0$,

until certain convergence criterions are met. The iterative procedure would
provide, under the Bayesian specification we described, the optimal parameters
for our problem.

Assume that the optimization problem for $(\bm{\omega},\theta)$ stopped at a
mode of the penalized loglikelihood for the selected parameter $\lambda$. The
final value of $\bm{H}(\bm{\omega},\theta,\lambda)$ could be used to construct
confidence intervals through, again, Laplace approximations. For a real-valued
function $(\bm{\omega},\bm{\theta})\mapsto
f(\bm{\omega},\bm{\theta})\in\mathbb{R}$, a confidence interval for
$f(\bm{\omega},\bm{\theta})$ can be obtained by a simple Gaussian
approximation and the Delta method [25]. The Delta method gives an
approximation for the variance of $f(\bm{\omega},\bm{\theta})$ as
$N^{-1}\lVert\frac{\partial}{\partial(\bm{\omega},\theta)}f(\bm{\omega},\theta)\rVert_{\bm{H}(\bm{\omega},\theta,\lambda)}$.
This method provides fast and reliable confidence intervals for any quantity,
e.g., for the density itself (pointwise), the quantile function (again,
pointwise), or the VaR / TVaR approximations at higher levels.

##  3 Numerical examples

In this section, we provide a few numerical examples of our approach. Recall
that our final algorithm as the following form :

Input: Information about the partition and the moments
$\mathcal{B},\bm{k},\hat{\bm{\pi}},\hat{\bm{\mu}}$.

Output: Estimated parameters $\hat{\theta},\hat{\bm{\omega}}$.

Initialize $\lambda=0$ and $\bm{\omega},\theta$ with random (positive) noise.

while  _not converged_ do

$\bm{\omega},\theta\leftarrow\arg\min\ell(\bm{\omega},\theta,\lambda)$

/* using e.g., BFGS */

$\bm{H}\leftarrow-\frac{\partial^{2}}{\partial^{2}\bm{\omega}}\ell(\bm{\omega},\theta,0)$

$\eta=\mathrm{eigenvalues}(-\bm{H}^{-1}\bm{P}_{r})$

$\lambda=\arg\max\ell(\lambda)$

/* Using Section 2 */

end while

return _$\bm{\omega},\theta$_

Algorithm 1 Constrained density estimation using Erlang mixtures.

Note that the computation of the Hessian $\bm{H}$ can be extracted from the
output of the solver used at the previous step, if a second order solver is
used. If a first order solver is used, like BFGS, an approximation of the
inverse Hessian is usually available at the end of the routine. In practice,
our implementation using the Julia language uses convergence criterions on
$\bm{\omega},\theta,\lambda$. Moreover, a small $\varepsilon$ can be added to
the diagonal of $\bm{H}$ to improve the condition number of the Hessian matrix
and facilitate the numerical inversion. Technical implementation relaxes the
positivity constraint on the parameters by a square transformation.

###  3.1 Lognormal example

To study the performance and behavior of the proposed scheme, we first revisit
a simple example from [5], the log-Normal case. Although the log-normal
density is unimodal, and our approximation class is clearly targeted at
multimodal behaviors, the log-normal density is also heavy tailed and thus a
pertinent benchmark for the practitioner. Of course, when the number of atoms
of $\nu$ goes to infinity the log-normal is inside our class, and we thus
expect good results. A proper analytic expression for $\nu$ can be obtained
using the Thorin measure of the log-normal distribution from [2, 11, 20, 19]
and then the mixture representation from [1, Theorem 4.1.1] (or the one from
[24]). For practical reasons however, finite mixtures are preferable as they
allow computations to be carried out.

######  Setting 3.1 (Lognormal).

We sample $N=750$ independent observations from a random variable
$X\sim\texttt{LogNormal}(0,0.5)$. The partition is set arbitrarily at the
quantiles of the sample at levels $\alpha\in\\{0,0.5,0.9,0.99,1\\}$, so that
it has $J=4$ bins, and we set $\bm{k}=(4,4,4,1)$ for the number of observed
moments within each bin. The corresponding dataset was given in Table 1.

The choice of partition in 3.1 separates the tail from the bulk of the
distribution. We run our algorithm with $n=50$ atoms for $\nu$ (that is,
$\nu$’s support is contained in $\\{1,...,n\\}$). The obtained results are
pictured in Figure 1.

From the set of moments depicted in Table 1, we obtain the optimal Erlang
mixture depicted in Figure 1 (a, b, c). Figure 1 (d) shows the convergence of
the procedure w.r.t. the number of empirical moments in the first three bins
(1,2,3 or 4) used to estimate the density, as assessed using $S=50$ resamples
from the lognormal model. It shows that the distance between the density
estimate and the true density (as measured by $\ell_{1}$ and $\ell_{2}$
distances between quantile functions and distribution functions, as well as
the Kullback-Leibler divergence) decreases significantly with the number of
available empirical local moments. The definition of these performance metrics
is recalled in Appendix B.

Figure 1:  General results on 3.1. The histogram in (a) represents
$\hat{\bm{\pi}}$, while the red curve represents the density estimate and the
red region the (pointwise) confidence interval around the estimated density.
Couples $(x_{i},\omega_{i}y_{i})$ from Equation (5) are represented in panel
(b). Panel (c) represents quantile-quantile plots of the estimated density
against the underlying raw data and the true density. Panels (d) represents
the median of the distance statistics across resamples for an increasing
number of observed local moments, renormalized to be $1$ for the least
informative setting.

We see on Figure 1 (a) that the kernel density estimate (black curve) obtained
on the raw (non-grouped) data and our estimate based on Erlang mixtures (red
curve) from grouped data in Table 1 both have troubles reproducing correctly
the value at the mode of the true density (blue curve). Our Laplace-
approximated confidence point-wise confidence interval confirms the local
difficulty. From the quantile-quantile plots on Figure 1 (c) we see that the
fit is satisfactory, except in the far tail (last $5$ to $10$ points over
$N=750$). A Kolmogorov-Smirnov test between the true underlying density and
the estimated density has a (two-sided) p-value of $0.631$, failing to reject
the null hypothesis that the distributions are the same at that sample size.

To observe the convergence of our estimator with respect to the number of
observations, we used $S=50$ resamples from 3.1 with different number of
observations $N$. We report in Table 2 the mean value, bias, standard
deviation and RMSE of quantiles (at several levels $\alpha$) obtained from our
resamples.

Table 2: Mean, bias, standard deviation and Root Mean Square Error of obtained quantiles at several levels $\alpha$ on resamples from 3.1, fitted on several resamples of varying sizes $N$. $\alpha$ |  0.5 |  0.9 |  0.95 |  0.99 |  0.995
---|---|---|---|---|---
True values |  $1.000$ |  $1.898$ |  $2.276$ |  $3.200$ |  $3.625$
Mean |  |  |  |  |
$N=250$ |  $1.021$ |  $1.877$ |  $2.307$ |  $3.169$ |  $3.434$
$N=500$ |  $1.006$ |  $1.909$ |  $2.402$ |  $3.253$ |  $3.497$
$N=750$ |  $1.000$ |  $1.890$ |  $2.398$ |  $3.276$ |  $3.519$
$N=1000$ |  $0.997$ |  $1.879$ |  $2.372$ |  $3.252$ |  $3.516$
$N=2000$ |  $1.002$ |  $1.881$ |  $2.368$ |  $3.247$ |  $3.515$
Std. dev. |  |  |  |  |
$N=250$ |  $0.033$ |  $0.112$ |  $0.272$ |  $0.442$ |  $0.469$
$N=500$ |  $0.021$ |  $0.097$ |  $0.183$ |  $0.298$ |  $0.337$
$N=750$ |  $0.021$ |  $0.081$ |  $0.185$ |  $0.206$ |  $0.234$
$N=1000$ |  $0.020$ |  $0.085$ |  $0.193$ |  $0.178$ |  $0.182$
$N=2000$ |  $0.012$ |  $0.059$ |  $0.179$ |  $0.130$ |  $0.161$
RMSE |  |  |  |  |
$N=250$ |  $0.040$ |  $0.114$ |  $0.274$ |  $0.443$ |  $0.507$
$N=500$ |  $0.022$ |  $0.098$ |  $0.223$ |  $0.303$ |  $0.361$
$N=750$ |  $0.021$ |  $0.081$ |  $0.222$ |  $0.219$ |  $0.257$
$N=1000$ |  $0.020$ |  $0.087$ |  $0.215$ |  $0.186$ |  $0.212$
$N=2000$ |  $0.012$ |  $0.061$ |  $0.201$ |  $0.138$ |  $0.195$

We note on Table 2, especially on the root-mean-square errors, that all
quantiles seem convergent with the number of observations. In particular,
quantiles higher than the last bin edge are also more and more precisely
estimated with the increasing number of observations.

In Figure 1 (d), a few key statistics (detailed in Appendix B) are averaged on
the resamples for several values of $\bm{k}$ (but fixed number of observations
$N$). This plot shows that increasing the number of moments in each bin does
have an impact on the quality of the estimator, at least on this example.
However, the limit is not exactly the lognormal distribution, as the log-
normal distributions cannot be characterized by its moments.

Box plots of these key statistics, for a given value of $\bm{k}$ but several
values of $N$ across the resamples are available in Figure 2(b), with clear
indications of convergence.

(a)

(b)

Figure 2: Comparison of the bootstrap statistics obtained on 3.1. Panel (a) :
for an increasing number of moments with
$\bm{k}\in\\{(1,1,1,1),(2,2,2,1),(3,3,3,1),(4,4,4,1)\\}$ and a fixed number of
observations $N$. Panel (b): for an increasing number of observations
$N\in\\{250,500,750,1000,2000\\}$ and fixed $\bm{k}=(4,4,4,1)$.

While this example is overall conclusive on the quality of the estimator, we
propose a few more examples with additional specific features (such as
multimodality) for the target density.

###  3.2 A mixture of gamma distributions

######  Setting 3.2 (Mixture of gammas).

We sample $N=2000$ independent observations from a random variable $X$ that
follows a $\frac{1}{3}$ — $\frac{2}{3}$ mixture of $\texttt{Gamma}(30,1)$ and
$\texttt{Gamma}(7,1)$. The partition is set at the quantiles of the sample at
levels $\alpha\in\\{0,0.5,0.9,0.99,1\\}$, so that it has $J=4$ bins, and we
use $\bm{k}=(4,4,4,1)$ for the number of observed moments within each bin.

Our estimator should be able to recover the features of the underlying
distribution fairly correctly since the true distribution actually belongs to
the class of Erlang mixtures. However, identifying the two atoms of the
underlying measure is not possible in our $L_{2}$-penalized settings222Finding
a few atoms would require instead some kind of $L_{1}$ penalization, see e.g.,
[31, 32] and their iSCAD penalty, as confirmed by the results in Figure 3.

Figure 3:  General results on 3.2. The histogram in (a) represents
$\hat{\bm{\pi}}$, while the red curve represents the density estimate and the
red region the (pointwise) confidence interval around the estimated density.
Couples $(x_{i},\omega_{i}y_{i})$ from Equation (5) are represented in panel
(b). Panel (c) represents quantile-quantile plots of the estimated density
against the underlying raw data and the true density.

We observe a particular difficulty in the estimation of the first mode of the
density on Figure 3, but our estimator performs better than the standard
Gaussian kernel — that has all the observation available and not only the
boxed moments. On the other hand, Kolmogorov-Smirnov tests are conclusive
(p-value of $0.157$) on this example.

###  3.3 A more challenging example

To confront our estimator to a dataset that is obviously outside its range, we
construct 3.3.

######  Setting 3.3 (Mixture of a Gaussian and a (reversed) Gamma [17]).

We sample $N=500$ independent observations from a random variable $X$ that
follows a 20/80 mixture of $\texttt{Normal}(1,1/3)$ and a reversed
$\texttt{Gamma}(11,1/6)$. The partition is set arbitrarily at the quantiles of
the sample at levels $\alpha\in\\{0,0.5,0.9,0.99,1\\}$, so that it has $J=4$
bins, and we use $\bm{k}=(4,4,4,1)$ for the number of observed moments in each
bin.

Figure 4:  General results on 3.3. The histogram in (a) represents
$\hat{\bm{\pi}}$, while the red curve represents the density estimate and the
red region the (pointwise) confidence interval around the estimated density.
Couples $(x_{i},\omega_{i}y_{i})$ from Equation (5) are represented in panel
(b). Panel (c) represents quantile-quantile plots of the estimated density
against the underlying raw data and the true density. Panels (d) represents
the median of the distance statistics across resamples for an increasing
number of observed local moments, renormalized to be $1$ for the least
informative setting.

(a)

(b)

Figure 5: Comparison of the bootstrap statistics obtained on 3.3. Panel (a) : for several parameters $\bm{k}\in\\{(1,1,1,1),(2,2,2,1),(3,3,3,1),(4,4,4,1)\\}$ and fixed number of observations $N$. Panel (b): for several number of observations $N\in\\{250,500,1000,2000\\}$ and fixed $\bm{k}=(4,4,4,1)$. Table 3: Mean, bias, standard deviation and Root Mean Square Error of obtained quantiles at several levels $\alpha$ on resamples from 3.1, fitted on several resamples of varying sizes $N$ on 3.3. $\alpha$ |  0.5 |  0.9 |  0.95 |  0.99 |  0.995
---|---|---|---|---|---
True values |  $3.643$ |  $4.376$ |  $4.530$ |  $4.778$ |  $4.857$
Mean |  |  |  |  |
$N=250$ |  $3.589$ |  $4.469$ |  $4.715$ |  $5.189$ |  $5.366$
$N=500$ |  $3.606$ |  $4.483$ |  $4.729$ |  $5.202$ |  $5.378$
$N=1000$ |  $3.612$ |  $4.490$ |  $4.737$ |  $5.211$ |  $5.390$
$N=2000$ |  $3.610$ |  $4.480$ |  $4.726$ |  $5.199$ |  $5.377$
Std. dev. |  |  |  |  |
$N=250$ |  $0.054$ |  $0.048$ |  $0.050$ |  $0.056$ |  $0.059$
$N=500$ |  $0.039$ |  $0.042$ |  $0.045$ |  $0.054$ |  $0.058$
$N=1000$ |  $0.028$ |  $0.032$ |  $0.036$ |  $0.044$ |  $0.045$
$N=2000$ |  $0.014$ |  $0.023$ |  $0.027$ |  $0.034$ |  $0.038$
RMSE |  |  |  |  |
$N=250$ |  $0.077$ |  $0.105$ |  $0.191$ |  $0.414$ |  $0.513$
$N=500$ |  $0.054$ |  $0.115$ |  $0.204$ |  $0.427$ |  $0.525$
$N=1000$ |  $0.042$ |  $0.119$ |  $0.210$ |  $0.435$ |  $0.535$
$N=2000$ |  $0.037$ |  $0.107$ |  $0.198$ |  $0.422$ |  $0.522$

We observe on Table 3 that the bootstrap standard error of the quantiles’
estimation is decreasing with the number of observations, but that the root-
mean-square error for extreme quantiles of the target distribution, which is
based on a comparison of the estimated quantity with its true value, does not.
Hence, there is here an error that does not decrease with the number of
observations for the estimation of extreme quantiles, see on Figure 4, panel
(a). Still, generic distances from Figure 5(b) are nevertheless correctly
decreasing. The Kolmogorov-Smirnov p-value of our estimator on 3.3 with
$N=500$ is at $0.048$.

###  3.4 A mixture of a (truncated) normal and a beta random variables.

######  Setting 3.4 (Mixture of Gaussian and Beta).

We sample $N=750$ independent observations from a random variable $X$ that
follows a 20/80 mixture of $\texttt{Normal}(1,1/3)$ and $\texttt{Beta}(4,8)$.
The partition is set arbitrarily as the quantiles of the sample at levels
$\alpha\in\\{0,0.5,0.9,0.99,1\\}$, so that it has $J=4$ bins, and we use
$\bm{k}=(4,4,4,1)$ for the number of observed moments within each bin.

Figure 6:  General results on 3.4. The histogram in (a) represents
$\hat{\bm{\pi}}$, while the red curve represents the density estimate and the
red region the (pointwise) confidence interval around the estimated density.
Couples $(x_{i},\omega_{i}y_{i})$ from Equation (5) are represented in panel
(b). Panel (c) represents quantile-quantile plots of the estimated density
against the underlying raw data and the true density. Panels (d) represents
the median of the distance statistics across resamples for an increasing
number of observed local moments, renormalized to be $1$ for the least
informative setting.

3.4, and in particular panel (a) of Figure 6 highlights the smoothing effect
of the $L_{2}$ penalization. Indeed, by comparing the black curve (the
Gaussian kernel that has access to (non-grouped) raw data) with the red curve
(our estimator), we clearly see the impact of smoothing. The statistics from
panel (d) are not obvious to interpret, but the p-value of the Kolmogorov-
Smirnov test of $0.306$ confirms the quality of the estimated density.

##  4 Discussion

Boxed moments problems on the form of Table 1 may arise for different reasons.
Confidentiality issues, for example, is particularly present in the insurance
and reinsurance sector. The boxed moments structure is quite permissive as it,
e.g., allows to provide a few moments for the _attritional_ part, the Value-
at-Risk, and potentially the Tail-Value-at-Risk, which all can generally be
obtained directly for already-present databases. The boxing also improves the
quality of the transmitted information without increasing the moment’s order,
which could lead to numerical overflows and/or underflows (see e.g., [5]).

If [17] dealt with the bounded case, we here provide a new estimator that
allows for estimation on the (positive) real line. The class of Erlang
mixtures was a natural pick due to their denseness into the class of
distributions supported on $\mathbb{R}_{+}$, and their previous use is
insurance losses modeling. These distributions provide closed-form boxed
moments, a characteristic that contributes to the performance of our Julian
implementation.

The very nature of our estimation problem requires regularization. In a
setting where cross-validation is not possible, we adapted Bayesian techniques
from P-splines discrete penalization. We derived an explicit algorithm that
leverages this Bayesian nature to automatically tune the penalty, leaving less
hyperparameters to the user. Of course, the number of Erlang distributions in
the mixtures (our parameter $n$) is another important hyperparameter, but due
to the denseness of the class the simple rule of _bigger $n$ is better_ can be
used.

Our approach could be easily generalized to multivariate settings for
estimation on $\mathbb{R}_{+}^{d}$. For that, we could leverage the approach
from [31, 32], granted in quite different settings ($L_{1}$-penalized settings
on full datasets instead of $L_{2}$-penalized settings on restricted
information). Moreover, we think that concurrent $L_{1}$ and $L_{2}$
penalization could be beneficial to the sparsity _and_ the smoothness of the
obtained estimator, if both criterions are of interests to the user.

##  5 Acknowledgements

Authors acknowledge the support of the ARC project IMAL (grant 20/25-107)
financed by the Wallonia-Brussels Federation and granted by the Académie
Universitaire Louvain.

##  Appendix A Commonly used probability distributions

We define here the probability distributions used in this paper and there
parametrization by providing a characterizing function for each of them:

  * •

$\texttt{Normal}(\bm{\mu},\bm{\Sigma})$ has density
$f(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^{2}}$,
for all $\mu\in\mathbb{R},\sigma\in\mathbb{R}_{+}$.

  * •

$\texttt{LogNormal}(\mu,\sigma)$ has density
$f(x)=\frac{1}{x\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(\log(x)-\mu)^{2}}{2\sigma^{2}}\right),\quad
x>0,$ for all $\mu\in\mathbb{R},\sigma\in\mathbb{R}_{+}$.

  * •

$\texttt{Gamma}(\alpha,s)$ has density
$f(x)=\frac{1}{\Gamma(\alpha)s^{\alpha}}x^{\alpha-1}e^{-\frac{x}{s}}$, for all
$\alpha,s\in\mathbb{R}_{+}$ and $\Gamma$ is the gamma function.

  * •

$\texttt{Beta}(a,b)$ has density $f(x)=\frac{x^{a-1}(1-x)^{b-1}}{B(a,b)}$,
where $a,b\in\mathbb{R}_{+}$ and $B$ is the Beta function.

  * •

$\texttt{Multinomial}(\bm{\pi},N)$ has probability mass function
$p(\bm{n})=\frac{n!}{n_{1}!,...,n_{m}!}\bm{\pi}^{\bm{n}}$, where $\bm{\pi}$ is
a probability vector and $\bm{n}$ is a vector of integers summing to $N$.

  * •

$\texttt{MixedGamma}(\nu,s)$ has density $f(x)=\int
f_{\texttt{Gamma}(\alpha,s)}(x)\nu(d\alpha)$ for
$\nu\in\mathcal{M}_{+}(\mathbb{R}_{+})$ and $s\in\mathbb{R}_{+}$.

  * •

$\texttt{MixedErlang}(\bm{\omega},s)=\texttt{MixedGamma}(\sum_{k=1}^{n}\omega_{k}\delta_{k},s)$,
where $n$ is the length of the probability vector $\omega$.

Loglikelihood of these distributions is derived simply by letting
$\ell=\log\circ f$. More details on these distributions can be found in the
literature.

##  Appendix B Key evaluation statistics

The paper leveraged a few key statistics to visualize the quality of the
resampled models:

  * •

The $\ell_{2}$ squared distance between quantile function $q_{1}$ and $q_{2}$:

$\ell_{2}(q_{1},q_{2})=\int_{0}^{1}(q_{1}(x)-q_{2}(x))^{2}dx.$

  * •

The $\ell_{2}$ squared distance between distribution function $F_{1}$ and
$F_{2}$:

$\ell_{2}(F_{1},F_{2})=\int(F_{1}(x)-F_{2}(x))^{2}dx.$

  * •

The $\ell_{1}$ distance between quantile function $q_{1}$ and $q_{2}$ (it is
technically equal to the $\ell_{1}$ distance between distribution functions
$F_{1}$ and $F_{2}$.):

$\ell_{1}(q_{1},q_{2})=\int_{0}^{1}\mathrm{abs}(q_{1}(x)-q_{2}(x))dx.$

  * •

The Kullback-Leibler divergence between density functions $f_{1},f_{2}$:

$KL(f_{1},f_{2})=\int f_{1}(x)\ln\left(\frac{f_{1}(x)}{f_{2}(x)}\right)dx.$

## References

  * [1] Lennart Bondesson  “Generalized Gamma Convolutions and Related Classes of Distributions and Densities” 76, Lecture Notes in Statistics  New York, NY: Springer New York, 1992  DOI: 10.1007/978-1-4612-2948-3
  * [2] Lennart Bondesson  “On the Lévy Measure of the Lognormal and the Logcauchy Distributions”  In _Methodology and Computing in Applied Probability_ 4.3 Springer, 2002, pp. 243–256  DOI: 10.1023/A:1022533817579
  * [3] N.. Breslow and D.. Clayton  “Approximate Inference in Generalized Linear Mixed Models”  In _Journal of the American Statistical Association_ 88.421, 1993, pp. 9–25  DOI: 10.1080/01621459.1993.10594284
  * [4] P Brockett et al.  “Actuarial Usage of Grouped Data: An Approach to Incorporating Secondary Data.”  In _Insurance Mathematics and Economics_ 2.19, 1997, pp. 147
  * [5] Hélène Cossette, David Landriault, Etienne Marceau and Khouzeima Moutanabbir  “Moment-Based Approximation with Mixed Erlang Distributions”  In _Variance_ 10.1, 2016, pp. 161–182
  * [6] Hélène Cossette, Etienne Marceau, Itre Mtalai and Déry Veilleux  “Dependent Risk Models with Archimedean Copulas: A Computational Strategy Based on Common Mixtures and Applications”  In _Insurance: Mathematics and Economics_ 78, 2018, pp. 53–71  DOI: 10.1016/j.insmatheco.2017.11.002
  * [7] Jacob Dexe, Jonas Ledendal and Ulrik Franke  “An Empirical Investigation of the Right to Explanation under GDPR in Insurance”  In _Trust, Privacy and Security in Digital Business: 17th International Conference, TrustBus 2020, Bratislava, Slovakia, September 14–17, 2020, Proceedings 17_ Springer, 2020, pp. 125–139
  * [8] Paul H.. Eilers and Brian D. Marx  “Flexible Smoothing with B-splines and Penalties”  In _Statistical Science_ 11.2, 1996  DOI: 10.1214/ss/1038425655
  * [9] Paul HC Eilers and Brian D Marx  “Practical Smoothing: The Joys of P-splines”  Cambridge University Press, 2021
  * [10] L. Fahrmeir and S. Lang  “Bayesian Smoothing and Regression for Longitudinal, Spatial and Event History Data”  Springer Science & Business Media, 2001
  * [11] Edward Furman, Daniel Hackmann and Alexey Kuznetsov  “On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination”  In _SSRN Electronic Journal_ , 2017  DOI: 10.2139/ssrn.3034540
  * [12] Lan Gong  “Erlang-Based Methods in Modeling Losses in Insurance and Applications”, 2014
  * [13] Oswaldo Gressani and Philippe Lambert  “Laplace Approximations for Fast Bayesian Inference in Generalized Additive Models Based on P-splines”  In _Computational Statistics & Data Analysis_ 154 Elsevier, 2021, pp. 107088  DOI: 10.1016/j.csda.2020.107088
  * [14] Wenyong Gui, Rongtan Huang and X. Lin  “Fitting Multivariate Erlang Mixtures to Data: A Roughness Penalty Approach”  In _Journal of Computational and Applied Mathematics_ 386, 2021, pp. 113216  DOI: 10.1016/j.cam.2020.113216
  * [15] Wenyong Gui, Rongtan Huang and X. Lin  “Fitting the Erlang Mixture Model to Data via a GEM-CMM Algorithm”  In _Journal of Computational and Applied Mathematics_ 343, 2018, pp. 189–205  DOI: 10.1016/j.cam.2018.04.032
  * [16] Mary A Johnson and Michael R Taaffe  “Matching Moments to Phase Distributions: Mixtures of Erlang Distributions of Common Order”  In _Stochastic Models_ 5.4 Taylor & Francis, 1989, pp. 711–743  DOI: 10.1080/15326348908807131
  * [17] Philippe Lambert  “Nonparametric Density Estimation and Risk Quantification from Tabulated Sample Moments”  In _Insurance: Mathematics and Economics_ 108 Elsevier, 2023, pp. 177–189  DOI: 10.1016/j.insmatheco.2022.12.004
  * [18] Philippe Lambert and Oswaldo Gressani  “Penalty Parameter Selection and Asymmetry Corrections to Laplace Approximations in Bayesian P-splines Models”  In _Statistical Modelling_ 23.5-6 SAGE Publications Sage India: New Delhi, India, 2023, pp. 409–423  DOI: 10.1177/1471082X231181173
  * [19] Oskar Laverny  “Estimation of High Dimensional Gamma Convolutions through Random Projections”  arXiv, 2022  DOI: 10.48550/ARXIV.2203.13741
  * [20] Oskar Laverny, Esterina Masiello, Véronique Maume-Deschamps and Didier Rullière  “Estimation of Multivariate Generalized Gamma Convolutions through Laguerre Expansions.”  In _Electronic Journal of Statistics_ 15.2, 2021  DOI: 10.1214/21-EJS1918
  * [21] Simon Lee and X Sheldon Lin  “Modeling Dependent Risks with Multivariate Erlang Mixtures”  In _Astin Bulletin_ 42.1, 2012, pp. 153–180
  * [22] Simon C.. Lee and X. Lin  “Modeling and Evaluating Insurance Losses Via Mixtures of Erlang Distributions”  In _North American Actuarial Journal_ 14.1, 2010, pp. 107–130  DOI: 10.1080/10920277.2010.10597580
  * [23] Xenofon Liapakis  “A GDPR Implementation Guide for the Insurance Industry”  In _International Journal of Reliable and Quality E-Healthcare (IJRQEH)_ 7.4 IGI Global, 2018, pp. 34–44  DOI: 10.4018/IJRQEH.2018100103
  * [24] P.. Moschopoulos  “The Distribution of the Sum of Independent Gamma Random Variables”  In _Annals of the Institute of Statistical Mathematics_ 37.3, 1985, pp. 541–544  DOI: 10.1007/BF02481123
  * [25] Gary W Oehlert  “A Note on the Delta Method”  In _The American Statistician_ 46.1 Taylor & Francis, 1992, pp. 27–29
  * [26] Lysa Porth, Wenjun Zhu and Ken Seng Tan  “A Credibility-Based Erlang Mixture Model for Pricing Crop Reinsurance”  In _Agricultural Finance Review_ 74.2 Emerald Group Publishing Limited, 2014, pp. 162–187  DOI: 10.1108/AFR-04-2014-0006
  * [27] H. Rue and L. Held  “Gaussian Markov Random Fields: Theory and Applications”  CRC Press, 2005
  * [28] L. Tierney and J.. Kadane  “Accurate Approximations for Posterior Moments and Marginal Densities”  In _Journal of the American Statistical Association_ 81.393 Taylor & Francis, 1986, pp. 82–86  DOI: 10.1080/01621459.1986.10478240
  * [29] Henk C Tijms  “Stochastic Models: An Algorithmic Approach”  John Wiley & Sons Incorporated, 1994
  * [30] Simon N Wood  “Generalized Additive Models: An Introduction with R”  CRC press, 2017
  * [31] Cuihong Yin and X. Lin  “Efficient Estimation of Erlang Mixtures Using iSCAD Penalty with Insurance Application”  In _ASTIN Bulletin_ 46.3, 2016, pp. 779–799  DOI: 10.1017/asb.2016.14
  * [32] Cuihong Yin, X. Sheldon Lin, Rongtan Huang and Haili Yuan  “On the Consistency of Penalized MLEs for Erlang Mixtures”  In _Statistics & Probability Letters_ 145, 2019, pp. 12–20  DOI: 10.1016/j.spl.2018.08.004


* [127] Xinmu Wang, Tatini Mal-Sarkar, Aswin Krishna, Seetharam Narasimhan, and Swarup Bhunia.  Software exploitable hardware trojans in embedded processor.  In 2012 IEEE International Symposium on Defect and Fault Tolerance in VLSI and Nanotechnology Systems (DFT), pages 55–58. IEEE, 2012\.
  * [128] Rahul Sasi.  Maldrone the first backdoor for drones.  Fb1h2s aka Rahul Sasi’s Blog, 2015.
  * [129] Víctor Mayoral Vilches, Laura Alzola Kirschgens, Asier Bilbao Calvo, Alejandro Hernández Cordero, Rodrigo Izquierdo Pisón, David Mayoral Vilches, Aday Muñiz Rosas, Gorka Olalde Mendia, Lander Usategi San Juan, Irati Zamalloa Ugarte, et al.  Introducing the robot security framework (rsf), a standardized methodology to perform security assessments in robotics.  arXiv preprint arXiv:1806.04042, 2018.
  * [130] Laura Alzola Kirschgens, Irati Zamalloa Ugarte, Endika Gil Uriarte, Aday Muniz Rosas, and Víctor Mayoral Vilches.  Robot hazards: from safety to security.  arXiv preprint arXiv:1806.06681, 2018.
  * [131] Marcello Pogliani, Davide Quarta, Mario Polino, Martino Vittone, Federico Maggi, and Stefano Zanero.  Security of controlled manufacturing systems in the connected factory: The case of industrial robots.  Journal of Computer Virology and Hacking Techniques, 15:161–175, 2019.
  * [132] John R Douceur.  The sybil attack.  In International workshop on peer-to-peer systems, pages 251–260. Springer, 2002.
  * [133] Giovanni Lacava, Angelica Marotta, Fabio Martinelli, Andrea Saracino, Antonio La Marra, Endika Gil-Uriarte, and Victor Mayoral Vilches.  Cybsersecurity issues in robotics.  J. Wirel. Mob. Networks Ubiquitous Comput. Dependable Appl., 12(3):1–28, 2021.
  * [134] Davide Quarta, Marcello Pogliani, Mario Polino, Federico Maggi, Andrea Maria Zanchettin, and Stefano Zanero.  An experimental security analysis of an industrial robot controller.  In 2017 IEEE Symposium on Security and Privacy (SP), pages 268–286, 2017.
  * [135] Tamara Bonaci, J. Herron, Tariq Yusuf, Junjie Yan, T. Kohno, and H. Chizeck.  To make a robot secure: An experimental analysis of cyber security threats against teleoperated surgical robots.  ArXiv, abs/1504.04339, 2015.
  * [136] H. Alemzadeh, Daniel Chen, Xiao Li, T. Kesavadas, Z. Kalbarczyk, and R. Iyer.  Targeted attacks on teleoperated surgical robots: Dynamic model-based detection and mitigation.  2016 46th Annual IEEE/IFIP International Conference on Dependable Systems and Networks (DSN), pages 395–406, 2016.
  * [137] Pinyao Guo, Hunmin Kim, Nurali Virani, J. Xu, Minghui Zhu, and Peng Liu.  Exploiting physical dynamics to detect actuator and sensor attacks in mobile robots.  ArXiv, abs/1708.01834, 2017.
  * [138] Y. Gao, Guanghui Sun, Jianxing Liu, Yang Shi, and Ligang Wu.  State estimation and self-triggered control of cpss against joint sensor and actuator attacks.  Autom., 113:108687, 2020.
  * [139] J. Giraldo, S. Kafash, Justin Ruths, and A. Cárdenas.  Daria: Designing actuators to resist arbitrary attacks against cyber-physical systems.  2020 IEEE European Symposium on Security and Privacy (EuroS&P), pages 339–353, 2020.
  * [140] Xin Huang and Jiuxiang Dong.  Reliable control policy of cyber-physical systems against a class of frequency-constrained sensor and actuator attacks.  IEEE Transactions on Cybernetics, 48:3432–3439, 2018.
  * [141] Hangfeng He, Wenhai Qi, and Y. Kao.  Hmm-based adaptive attack-resilient control for markov jump system and application to an aircraft model.  Appl. Math. Comput., 392:125668, 2021.
  * [142] Cédric Escudero, Carlos Murguia, Paolo Massioni, and Eric Zamaï.  Enforcing safety under actuator injection attacks through input filtering.  In 2022 European Control Conference (ECC), pages 1521–1528. IEEE, 2022.
  * [143] Daehyung Park, Zackory Erickson, Tapomayukh Bhattacharjee, and Charles C Kemp.  Multimodal execution monitoring for anomaly detection during robot manipulation.  In 2016 IEEE International Conference on Robotics and Automation (ICRA), pages 407–414. IEEE, 2016.
  * [144] Pinyao Guo, Hunmin Kim, Nurali Virani, Jun Xu, Minghui Zhu, and Peng Liu.  Roboads: Anomaly detection against sensor and actuator misbehaviors in mobile robots.  In 2018 48th Annual IEEE/IFIP international conference on dependable systems and networks (DSN), pages 574–585. IEEE, 2018.
  * [145] Björn Browatzki, Vadim Tikhanoff, Giorgio Metta, Heinrich H. Bülthoff, and Christian Wallraven.  Active object recognition on a humanoid robot.  In 2012 IEEE International Conference on Robotics and Automation, pages 2021–2028, 2012.
  * [146] F. Alonso-Martín and Miguel A. Salichs.  Integration of a voice recognition system in a social robot.  Cybernetics and Systems, 42(4):215–245, 2011.
  * [147] Claudia Cruz, L. Enrique Sucar, and Eduardo F. Morales.  Real-time face recognition for human-robot interaction.  In 2008 8th IEEE International Conference on Automatic Face Gesture Recognition, pages 1–6, 2008.
  * [148] Marco Barreno, Blaine Nelson, Anthony D Joseph, and J D Tygar.  The security of machine learning.  Machine Learning, 81(2):121–148, 2010.
  * [149] Davide Castelvecchi.  Can we open the black box of ai?, 2016.
  * [150] In 2016, microsoft’s racist chatbot revealed the dangers of online conversation.  https://spectrum.ieee.org/in-2016-microsofts-racist-chatbot-revealed-the-dangers-of-online-conversation.  Accessed: 2021-08-23.
  * [151] How data poisoning attacks can corrupt machine learning models.  https://www.linkedin.com/pulse/how-data-poisoning-attacks-can-corrupt-machine-learning-misra.  Accessed: 2021-08-23.
  * [152] K. Liu, Brendan Dolan-Gavitt, and S. Garg.  Fine-pruning: Defending against backdooring attacks on deep neural networks.  In RAID, 2018.
  * [153] Micah Goldblum, D. Tsipras, Chulin Xie, Xinyun Chen, Avi Schwarzschild, D. Song, A. Madry, Bo Li, and T. Goldstein.  Dataset security for machine learning: Data poisoning, backdoor attacks, and defenses.  ArXiv, abs/2012.10544, 2020.
  * [154] Moses Charikar, Jacob Steinhardt, and Gregory Valiant.  Learning from untrusted data.  In Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2017, page 47–60, New York, NY, USA, 2017. Association for Computing Machinery.
  * [155] Giorgio Severi, Jim Meyer, Scott Coull, and Alina Oprea.  Explanation-Guided backdoor poisoning attacks against malware classifiers.  In 30th USENIX Security Symposium (USENIX Security 21), pages 1487–1504. USENIX Association, August 2021.
  * [156] Y. Liu, Y. Xie, and A. Srivastava.  Neural trojans.  In 2017 IEEE 35th International Conference on Computer Design (ICCD), pages 45–48, Los Alamitos, CA, USA, nov 2017. IEEE Computer Society.
  * [157] Matthew Jagielski, Alina Oprea, B. Biggio, Chang Liu, C. Nita-Rotaru, and Bo Li.  Manipulating machine learning: Poisoning attacks and countermeasures for regression learning.  2018 IEEE Symposium on Security and Privacy (SP), pages 19–35, 2018\.
  * [158] Liuwan Zhu, R. Ning, Cong Wang, Chunsheng Xin, and Hongyi Wu.  Gangsweep: Sweep out neural backdoors by gan.  Proceedings of the 28th ACM International Conference on Multimedia, 2020.
  * [159] Reza Shokri, Marco Stronati, Congzheng Song, and Vitaly Shmatikov.  Membership inference attacks against machine learning models.  In 2017 IEEE symposium on security and privacy (SP), pages 3–18. IEEE, 2017.
  * [160] Matt Fredrikson, Somesh Jha, and Thomas Ristenpart.  Model inversion attacks that exploit confidence information and basic countermeasures.  In Proceedings of the 22nd ACM SIGSAC conference on computer and communications security, pages 1322–1333, 2015.
  * [161] Yigitcan Kaya and T. Dumitras.  When does data augmentation help with membership inference attacks?  In ICML, 2021.
  * [162] Gabriel Resende Machado, Eugênio Silva, and Ronaldo Ribeiro Goldschmidt.  Adversarial machine learning in image classification: A survey toward the defender’s perspective.  ACM Comput. Surv., 55(1), nov 2021.
  * [163] Xiaoyi Chen, A. Salem, M. Backes, Shiqing Ma, and Y. Zhang.  Badnl: Backdoor attacks against nlp models.  ArXiv, abs/2006.01043, 2020.
  * [164] Aarti Kashyap, Syed Mubashir Iqbal, K. Pattabiraman, and Margo I. Seltzer.  Relusyn: Synthesizing stealthy attacks for deep neural network based cyber-physical systems.  ArXiv, abs/2105.10393, 2021.
  * [165] Jiawang Bai, Baoyuan Wu, Yong Zhang, Y. Li, Zhifeng Li, and Shutao Xia.  Targeted attack against deep neural networks via flipping limited weight bits.  ArXiv, abs/2102.10496, 2021.
  * [166] Xiangyu Qi, Jifeng Zhu, Chulin Xie, and Yong Yang.  Subnet replacement: Deployment-stage backdoor attack against deep neural networks in gray-box setting.  ArXiv, abs/2107.07240, 2021.
  * [167] Jonathan Pan.  Blackbox trojanising of deep learning models: Using non-intrusive network structure and binary alterations.  2020 IEEE REGION 10 CONFERENCE (TENCON), pages 891–896, 2020.
  * [168] Yunqian Wen, Li Song, Bo Liu, Ming Ding, and Rong Xie.  Identitydp: Differential private identification protection for face images.  ArXiv, abs/2103.01745, 2021.
  * [169] Maxim Maximov, Ismail Elezi, and L. Leal-Taix’e.  Ciagan: Conditional identity anonymization generative adversarial networks.  2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 5446–5455, 2020.
  * [170] Santiago Zanella Béguelin, L. Wutschitz, Shruti Tople, Victor Rühle, Andrew J. Paverd, O. Ohrimenko, Boris Köpf, and Marc Brockschmidt.  Analyzing information leakage of updates to natural language models.  Proceedings of the 2020 ACM SIGSAC Conference on Computer and Communications Security, 2020.
  * [171] Takashi Shibata, Go Irie, Daiki Ikami, and Yu Mitsuzumi.  Learning with selective forgetting.  In IJCAI, 2021.
  * [172] Shiqi Wang, Kexin Pei, J. Whitehouse, Junfeng Yang, and S. Jana.  Formal security analysis of neural networks using symbolic intervals.  In USENIX Security Symposium, 2018.
  * [173] Tribhuvanesh Orekondy, Bernt Schiele, and Mario Fritz.  Prediction poisoning: Towards defenses against dnn model stealing attacks.  arXiv preprint arXiv:1906.10908, 2019.
  * [174] Adi Ashkenazy and Shahar Zini.  Cylance, i kill you, 2019.
  * [175] Keywhan Chung, Xiao Li, Peicheng Tang, Zeran Zhu, Zbigniew T Kalbarczyk, Ravishankar K Iyer, and Thenkurussi Kesavadas.  Smart malware that uses leaked control data of robotic applications: The case of $\\{$Raven-II$\\}$ surgical robots.  In 22nd International symposium on research in attacks, intrusions and defenses (RAID 2019), pages 337–351, 2019.
  * [176] Víctor Mayoral-Vilches, Unai Ayucar Carbajo, and Endika Gil-Uriarte.  Industrial robot ransomware: Akerbeltz.  In 2020 fourth IEEE international conference on robotic computing (IRC), pages 432–435. IEEE, 2020.
  * [177] CG Leela Krishna and Robin R Murphy.  A review on cybersecurity vulnerabilities for unmanned aerial vehicles.  In 2017 IEEE International Symposium on Safety, Security and Rescue Robotics (SSRR), pages 194–199. IEEE, 2017.
  * [178] Gianmarco Veruggio.  The birth of roboethics.  2005\.
  * [179] Robot ethics.  https://www.ieee-ras.org/robot-ethics, 2023.  Accessed: 2023-08-20.
  * [180] Spyros G Tzafestas.  Roboethics.  A navigating overview. Heilberg: Springer, 2016.
  * [181] Accountability Act.  Health insurance portability and accountability act of 1996.  Public law, 104:191, 1996.
  * [182] Protection Regulation.  Regulation (eu) 2016/679 of the european parliament and of the council.  Regulation (eu), 679:2016, 2016.
  * [183] Hong Wang, Amir Khajepour, Dongpu Cao, and Teng Liu.  Ethical decision making in autonomous vehicles: Challenges and research progress.  IEEE Intelligent Transportation Systems Magazine, 14(1):6–17, 2020\.
  * [184] Noah J Goodall.  Can you program ethics into a self-driving car?  Ieee Spectrum, 53(6):28–58, 2016.
  * [185] Bryce Goodman and Seth Flaxman.  European union regulations on algorithmic decision-making and a “right to explanation”.  AI magazine, 38(3):50–57, 2017.
  * [186] Patrick Lin.  The ethics of saving lives with autonomous cars are far murkier than you think.  Wired Opinion, 30, 2013.
  * [187] Erin A McDaniel.  Robot wars: Legal and ethical dilemmas of using unmanned robotic systems in 21st century warfare and beyond.  PhD thesis, Fort Leavenworth, KS: US Army Command and General Staff College, 2008.
  * [188] Isaac Asimov.  Three laws of robotics.  Asimov, I. Runaround, 2, 1941.
  * [189] Peter M Asaro.  The liability problem for autonomous artificial agents.  In 2016 AAAI Spring symposium series, 2016.
  * [190] Shane O’Sullivan, Nathalie Nevejans, Colin Allen, Andrew Blyth, Simon Leonard, Ugo Pagallo, Katharina Holzinger, Andreas Holzinger, Mohammed Imran Sajid, and Hutan Ashrafian.  Legal, regulatory, and ethical frameworks for development of standards in artificial intelligence (ai) and autonomous robotic surgery.  The international journal of medical robotics and computer assisted surgery, 15(1):e1968, 2019.
  * [191] Preparing for the future of artificial intelligence.  https://obamawhitehouse.archives.gov/sites/default/files/whitehouse_files/microsites/ostp/NSTC/preparing_for_the_future_of_ai.pdf, 2023\.  Accessed: 2023-09-30.
  * [192] Kaori Ishii.  Comparative legal study on privacy and personal data protection for robots equipped with artificial intelligence: looking at functional and technological aspects.  AI & society, 34:509–533, 2019.
  * [193] Richmond Y Wong, Andrew Chong, and R Cooper Aspegren.  Privacy legislation as business risks: How gdpr and ccpa are represented in technology companies’ investment risk disclosures.  Proceedings of the ACM on Human-Computer Interaction, 7(CSCW1):1–26, 2023.
  * [194] Laura Graves, Vineel Nagisetty, and Vijay Ganesh.  Amnesiac machine learning.  In AAAI, 2021.
  * [195] E. F. Villaronga, Peter Kieseberg, and Tiffany C. Li.  Humans forget, machines remember: Artificial intelligence and the right to be forgotten.  Comput. Law Secur. Rev., 34:304–313, 2018.
  * [196] Aditya Golatkar, A. Achille, and Stefano Soatto.  Eternal sunshine of the spotless net: Selective forgetting in deep networks.  2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 9301–9309, 2020.
  * [197] Kelly Caine, Selma Šabanovic, and Mary Carter.  The effect of monitoring by cameras and robots on the privacy enhancing behaviors of older adults.  In Proceedings of the Seventh Annual ACM/IEEE International Conference on Human-Robot Interaction, HRI ’12, page 343–350, New York, NY, USA, 2012. Association for Computing Machinery.
  * [198] Delphine Reinhardt, Monisha Khurana, and Luca Hernández Acosta.  “i still need my privacy”: Exploring the level of comfort and privacy preferences of german-speaking older adults in the case of mobile assistant robots.  Pervasive and Mobile Computing, 74:101397, 2021.
  * [199] Jakub Złotowski, Kumar Yogeeswaran, and Christoph Bartneck.  Can we control it? autonomous robots threaten human identity, uniqueness, safety, and resources.  International Journal of Human-Computer Studies, 100:48–54, 2017\.
  * [200] P. A. Hancock, D. R. Billings, and K. E. Schaefer.  Can you trust your robot?  Ergonomics in Design, 19(3):24–29, 2011.
  * [201] Alan R Wagner and Ronald C Arkin.  Acting deceptively: Providing robots with the capacity for deception.  International Journal of Social Robotics, 3(1):5–26, 2011.
  * [202] Dimosthenis Kontogiorgos, Sanne van Waveren, Olle Wallberg, Andre Pereira, Iolanda Leite, and Joakim Gustafson.  Embodiment effects in interactions with failing robots.  In Proceedings of the 2020 CHI conference on human factors in computing systems, pages 1–14, 2020.
  * [203] Nicole Mirnig, Gerald Stollnberger, Markus Miksch, Susanne Stadler, Manuel Giuliani, and Manfred Tscheligi.  To err is robot: how humans assess and act toward an erroneous social robot.  Frontiers, 4(21):1, 2017.
  * [204] Emily McReynolds, Sarah Hubbard, Timothy Lau, Aditya Saraf, Maya Cakmak, and Franziska Roesner.  Toys That Listen: A Study of Parents, Children, and Internet-Connected Toys, page 5197–5207.  Association for Computing Machinery, New York, NY, USA, 2017.
  * [205] Ha Manh Do, Minh Pham, Weihua Sheng, Dan Yang, and Meiqin Liu.  Rish: A robot-integrated smart home for elderly care.  Robotics and Autonomous Systems, 101:74–92, 2018.
  * [206] Eric Zeng, Shrirang Mare, and Franziska Roesner.  End user security and privacy concerns with smart homes.  In Thirteenth Symposium on Usable Privacy and Security (SOUPS 2017), pages 65–80, Santa Clara, CA, July 2017. USENIX Association.
  * [207] Hasmik Atoyan, Jean-Rémi Duquet, and Jean-Marc Robert.  Trust in new decision aid systems.  In Proceedings of the 18th Conference on l’Interaction Homme-Machine, pages 115–122, 2006.
  * [208] John D Lee and Katrina A See.  Trust in automation: Designing for appropriate reliance.  Human factors, 46(1):50–80, 2004.
  * [209] Peter A Hancock, Deborah R Billings, and Kristen E Schaefer.  Can you trust your robot?  Ergonomics in Design, 19(3):24–29, 2011.
  * [210] Zahra Rezaei Khavas, S Reza Ahmadzadeh, and Paul Robinette.  Modeling trust in human-robot interaction: A survey.  In International conference on social robotics, pages 529–541. Springer, 2020.
  * [211] Maha Salem, Gabriella Lakatos, Farshid Amirabdollahian, and Kerstin Dautenhahn.  Would you trust a (faulty) robot? effects of error, task type and personality on human-robot cooperation and trust.  In Proceedings of the tenth annual ACM/IEEE international conference on human-robot interaction, pages 141–148, 2015.
  * [212] Michael W Boyce, Jessie YC Chen, Anthony R Selkowitz, and Shan G Lakhmani.  Effects of agent transparency on operator trust.  In Proceedings of the Tenth Annual ACM/IEEE International Conference on Human-Robot Interaction Extended Abstracts, pages 179–180, 2015\.
  * [213] Kornelia Lazanyi and Greta Maraczi.  Dispositional trust—do we trust autonomous cars?  In 2017 IEEE 15th International Symposium on Intelligent Systems and Informatics (SISY), pages 000135–000140. IEEE, 2017.
  * [214] Munjal Desai, Mikhail Medvedev, Marynel Vázquez, Sean McSheehy, Sofia Gadea-Omelchenko, Christian Bruggeman, Aaron Steinfeld, and Holly Yanco.  Effects of changing reliability on trust of robot systems.  In Proceedings of the seventh annual ACM/IEEE international conference on Human-Robot Interaction, pages 73–80, 2012.
  * [215] Paul Robinette, Alan R Wagner, and Ayanna M Howard.  Building and maintaining trust between humans and guidance robots in an emergency.  2013\.
  * [216] Benedict Tay, Younbo Jung, and Taezoon Park.  When stereotypes meet robots: the double-edge sword of robot gender and personality in human–robot interaction.  Computers in Human Behavior, 38:75–84, 2014.
  * [217] Paul Robinette, Wenchen Li, Robert Allen, Ayanna M Howard, and Alan R Wagner.  Overtrust of robots in emergency evacuation scenarios.  In 2016 11th ACM/IEEE International Conference on Human-Robot Interaction (HRI), pages 101–108. IEEE, 2016.
  * [218] Mohammad Obaid, Eduardo B Sandoval, Jakub Złotowski, Elena Moltchanova, Christina A Basedow, and Christoph Bartneck.  Stop! that is close enough. how body postures influence human-robot proximity.  In 2016 25th IEEE international symposium on robot and human interactive communication (RO-MAN), pages 354–361. IEEE, 2016.
  * [219] X Jessie Yang, Vaibhav V Unhelkar, Kevin Li, and Julie A Shah.  Evaluating effects of user experience and system transparency on trust in automation.  In Proceedings of the 2017 ACM/IEEE international conference on human-robot interaction, pages 408–416, 2017.
  * [220] Michael L Walters, Mohammedreza A Oskoei, Dag Sverre Syrdal, and Kerstin Dautenhahn.  A long-term human-robot proxemic study.  In 2011 RO-MAN, pages 137–142. IEEE, 2011.
  * [221] Scott Ososky, David Schuster, Elizabeth Phillips, and Florian G Jentsch.  Building appropriate trust in human-robot teams.  In 2013 AAAI spring symposium series, 2013.
  * [222] Simone Fischer-Hübner.  IT-security and privacy: design and use of privacy-enhancing security mechanisms.  Springer, 2001.
  * [223] Andreas Pfitzmann and Marit Hansen.  Anonymity, unlinkability, unobservability, pseudonymity, and identity management-a consolidated proposal for terminology.  2005\.
  * [224] Privacy enhancing technologies.  https://www.frbsf.org/banking/wp-content/uploads/sites/5/Privacy-Enhancing-Technologies_FINAL_V2_TOC-Update.pdf, 2023\.  Accessed: 2023-09-30.
  * [225] Nesrine Kaaniche, Maryline Laurent, and Sana Belguith.  Privacy enhancing technologies for solving the privacy-personalization paradox: Taxonomy and survey.  Journal of Network and Computer Applications, 171:102807, 2020.
  * [226] Stephen Eick and Annie I Antón.  Enhancing privacy in robotics via judicious sensor selection.  In 2020 IEEE International Conference on Robotics and Automation (ICRA), pages 7156–7165. IEEE, 2020.
  * [227] Björn Horstmann, Niels Diekmann, Hendrik Buschmeier, and Teena Hassan.  Towards designing privacy-compliant social robots for use in private households: A use case based identification of privacy implications and potential technical measures for mitigation.  In 2020 29th IEEE International Conference on Robot and Human Interactive Communication (RO-MAN), pages 869–876. IEEE, 2020.
  * [228] Daniel J Butler, Justin Huang, Franziska Roesner, and Maya Cakmak.  The privacy-utility tradeoff for remotely teleoperated robots.  In Proceedings of the tenth annual ACM/IEEE international conference on human-robot interaction, pages 27–34, 2015.
  * [229] Marco Dorigo, Guy Theraulaz, and Vito Trianni.  Reflections on the future of swarm robotics.  Science Robotics, 5(49):eabe4385, 2020.
  * [230] Sanjay Sarma Oruganti Venkata, Ramviyas Parasuraman, and Ramana Pidaparti.  Kt-bt: A framework for knowledge transfer through behavior trees in multirobot systems.  IEEE Transactions on Robotics, 2023.
  * [231] Chuanqiang Shan, Huiyun Jiao, and Jie Fu.  Towards representation identical privacy-preserving graph neural network via split learning.  arXiv preprint arXiv:2107.05917, 2021.
  * [232] Shaswata Mitra, Aritran Piplai, Sudip Mittal, and Anupam Joshi.  Combating fake cyber threat intelligence using provenance in cybersecurity knowledge graphs.  In 2021 IEEE International Conference on Big Data (Big Data), pages 3316–3323. IEEE, 2021.
  * [233] Tathagata Chakraborti, Yu Zhang, David E Smith, and Subbarao Kambhampati.  Planning with resource conflicts in human-robot cohabitation.  In Proceedings of the 2016 international conference on autonomous agents & multiagent systems, pages 1069–1077, 2016.
  * [234] Bing Cai Kok and Harold Soh.  Trust in robots: Challenges and opportunities.  Current Robotics Reports, 1:297–309, 2020.
  * [235] Aidan Fuller, Zhong Fan, Charles Day, and Chris Barlow.  Digital twin: Enabling technologies, challenges and open research.  IEEE access, 8:108952–108971, 2020.
  * [236] Arne Bilberg and Ali Ahmad Malik.  Digital twin driven human–robot collaborative assembly.  CIRP annals, 68(1):499–502, 2019.
  * [237] Baicun Wang, Huiying Zhou, Geng Yang, Xingyu Li, and Huayong Yang.  Human digital twin (hdt) driven human-cyber-physical systems: Key technologies and applications.  Chinese Journal of Mechanical Engineering, 35(1):11, 2022.
  * [238] Subash Neupane, Ivan A Fernandez, Wilson Patterson, Sudip Mittal, Milan Parmar, and Shahram Rahimi.  Twinexplainer: Explaining predictions of an automotive digital twin.  arXiv preprint arXiv:2302.00152, 2023.
  * [239] Antonio Roque and Suresh K Damodaran.  Explainable ai for security of human-interactive robots.  International Journal of Human–Computer Interaction, 38(18-20):1789–1807, 2022.
  * [240] Subash Neupane, Jesse Ables, William Anderson, Sudip Mittal, Shahram Rahimi, Ioana Banicescu, and Maria Seale.  Explainable intrusion detection systems (x-ids): A survey of current methods, challenges, and opportunities.  IEEE Access, 10:112392–112415, 2022.
  * [241] Devleena Das, Siddhartha Banerjee, and Sonia Chernova.  Explainable ai for robot failures: Generating explanations that improve user assistance in fault recovery.  In Proceedings of the 2021 ACM/IEEE International Conference on Human-Robot Interaction, HRI ’21, page 351–360, New York, NY, USA, 2021. Association for Computing Machinery.
  * [242] Kexin Chen, Tiffany Hwu, Hirak J Kashyap, Jeffrey L Krichmar, Kenneth Stewart, Jinwei Xing, and Xinyun Zou.  Neurorobots as a means toward neuroethology and explainable ai.  Frontiers in Neurorobotics, 14:570308, 2020.
  * [243] Nathan Fulton and André Platzer.  Verifiably safe off-model reinforcement learning.  In Tomáš Vojnar and Lijun Zhang, editors, Tools and Algorithms for the Construction and Analysis of Systems, pages 413–430, Cham, 2019. Springer International Publishing.
  * [244] Lukas Brunke, Melissa Greeff, Adam W Hall, Zhaocong Yuan, Siqi Zhou, Jacopo Panerati, and Angela P Schoellig.  Safe learning in robotics: From learning-based control to safe reinforcement learning.  Annual Review of Control, Robotics, and Autonomous Systems, 5:411–444, 2022.
  * [245] Felix Berkenkamp, Matteo Turchetta, Angela Schoellig, and Andreas Krause.  Safe model-based reinforcement learning with stability guarantees.  Advances in neural information processing systems, 30, 2017.
  * [246] Scott Thiebes, Sebastian Lins, and Ali Sunyaev.  Trustworthy artificial intelligence.  Electronic Markets, 31:447–464, 2021.
  * [247] Daniel Votipka, Eric Zhang, and Michelle L Mazurek.  Hacked: A pedagogical analysis of online vulnerability discovery exercises.  In 2021 IEEE Symposium on Security and Privacy (SP), pages 1268–1285. IEEE, 2021.
  * [248] Affan Yasin, Lin Liu, Tong Li, Jianmin Wang, and Didar Zowghi.  Design and preliminary evaluation of a cyber security requirements education game (sreg).  Information and Software Technology, 95:179–200, 2018.
  * [249] Sylvain Frey, Awais Rashid, Pauline Anthonysamy, Maria Pinto-Albuquerque, and Syed Asad Naqvi.  The good, the bad and the ugly: A study of security decisions in a cyber-physical systems game.  IEEE Transactions on Software Engineering, 45(5):521–536, 2019\.
  * [250] Association for advancing automation.  https://www.automate.org/a3-content/who-we-are, 2023.  Accessed: 2023-10-04.