#  Empirical Density Estimation based on Spline Quasi-Interpolation with
applications to Copulas clustering modeling

Cristiano Tamborrino
Department of Computer Science
University of Bari
Italy
<EMAIL_ADDRESS>
&Antonella Falini
Department of Computer Science
University of Bari
Italy
<EMAIL_ADDRESS>
&Francesca Mazzia
Department of Computer Science
University of Bari
Italy
<EMAIL_ADDRESS>

###### Abstract

Density estimation is a fundamental technique employed in various fields to
model and to understand the underlying distribution of data. The primary
objective of density estimation is to estimate the probability density
function of a random variable. This process is particularly valuable when
dealing with univariate or multivariate data and is essential for tasks such
as clustering, anomaly detection, and generative modeling. In this paper we
propose the mono-variate approximation of the density using spline quasi
interpolation and we applied it in the context of clustering modeling. The
clustering technique used is based on the construction of suitable
multivariate distributions which rely on the estimation of the monovariate
empirical densities (marginals). Such an approximation is achieved by using
the proposed spline quasi-interpolation, while the joint distributions to
model the sought clustering partition is constructed with the use of copulas
functions. In particular, since copulas can capture the dependence between the
features of the data independently from the marginal distributions, a finite
mixture copula model is proposed. The presented algorithm is validated on
artificial and real datasets.

_K_ eywords Spline Quasi-Interpolation, Copulas, Empirical Density Estimation,
Clustering

##  1 Introduction

Density estimation is a fundamental technique employed in various fields to
model and understand the underlying distribution of data. It plays a pivotal
role in capturing the inherent patterns and structures within a dataset,
making it a crucial component in statistical modeling, machine learning, and
data analysis. The primary objective of density estimation is to estimate the
probability density function (PDF) of a random variable. This process is
particularly valuable when dealing with univariate or multivariate data and is
essential for tasks such as clustering, anomaly detection, and generative
modeling. Several methods have been developed to address the challenge of
estimating the underlying density function. Classical approaches include
histogram-based methods or kernel density estimation (KDE). These methods vary
in complexity, assumptions, and computational efficiency, offering a range of
options to suit different types of data and analytical requirements [25, 40].
Other approaches in literature utilize splines [14, 38, 21]. Understanding the
strengths and limitations of each method is crucial for selecting an
appropriate density estimation technique based on the characteristics of the
data at hand. In this work, we aim to investigate a novel method for
estimating the probability density through the utilization of a technique
known as B-spline Hermite quasi-interpolation [30]. We will propose its
application for the development of a new Copula-Based clustering algorithm,
where density estimation plays a crucial role. Clustering is considered an
effective method for organizing data into groups based on the similarities in
features and characteristics among data points [17]. In recent years, various
algorithms have been developed for this purpose [1, 50, 51, 6, 10, 17, 22,
20]. For a comprehensive overview of the diverse approaches in the literature,
e.g., [48, 12, 13, 42]. One of the most well-known and widely used approaches
involves finite mixture modeling, which serves as a flexible and robust
probabilistic tool for both univariate and multivariate data. However, it is
important to note that Gaussian distributions, commonly used in finite
mixtures, may not always accurately represent real-world data aggregation. To
overcome this drawback, the use of alternative distributions, based on
copulas, have gained significant interest in recent times for their potential
to enhance the accuracy and applicability of clustering algorithms in a wider
range of scenarios.

Copulas offer an alternative and increasingly popular approach to modeling
data dependencies and aggregations [19, 35, 11]. They provide a more versatile
framework that does not rely on specific distributional assumptions like
Gaussian distributions. Copulas can capture complex and non-linear
relationships between variables, making them particularly useful for
situations where data aggregation behaviors deviate from traditional
distributions. The copulas have already been used as a possible solution to
the clustering problem, for example, Di Lascio et all. [10], address the
problem of the clusterization by using different family of copulas. In [44],
[36], the clusterization process is done by the assumption that the
multivariate dependencies is modelled by a new family of copula called
Gaussian Mixture Copula, in [22] the authors analyze models of finite mixtures
with different families of copulas, assuming that the marginals follow known
distributions. The widely used method to fit the finite mixture model to
observed data is the expectation maximization (EM) algorithm [9, 31], to
estimate the maximum likelihood. Variations and/or adaptations to special
situations of the EM algorithm exist, as the stochastic EM (SEM) algorithm
(e.g., [16, 28, 33]), the classification EM (CEM) algorithm (e.g.[6]), the
Monte CarloEM (MCEM) algorithm (e.g., [43, 46]) and those developed by [39,
32]. An important point to take into account is that the choice of the copula
model-based clustering imposes distributional assumptions on the marginals,
along each dimension, and these marginal distributions are assumed or forced
to be identical (e.g. a multivariate normal imposes univariate normal
distribution on each marginal); such assumptions restrict the modeling
flexibility. These restrictive assumptions could lead to erroneous modeling.
To address this issue, a semiparametric approach is employed, in which
marginal distributions are empirically estimated using kernel density
estimation (KDE) [15, 18]. In this work, we propose a new strategy that
enables the estimation of marginal densities through the use of the introduced
B-Spline Quasi-Interpolants (BSHQI) and we describe a mixture model for
density estimation based on Copulas that allows us to automatically choose a
different Copula for each cluster. The paper is organized as follows, in
Section 2, the BSHQI density estimation is described together with its
theoretical consistency properties and some statistical tests are conducted to
validate the proved theoretical results. In Section 3 we revise some
preliminary concepts related to copulas. In Section 4 the used EM algorithm
for Copula Mixture models is detailed. In Section 5 artificial and real
datasets are analysed and finally some conclusive remarks are presented in
Section 6.

##  2 BSHQI density estimation

Density estimation is a fundamental task in statistical analysis, involving
the determination of the underlying probability distribution for a set of
observed data.

Let $X_{1},\ldots,X_{n}$ be independent and identically distributed (i.i.d)
random variables with an unknown Cumulative Distribution Function (CDF)
$F(x)$. A non parametric estimator for $F(x)$ is provided by the Empirical
Cumulative Distribution Function (ECDF) $F_{n}(x)=\sum_{i=1}^{n}I(X_{i}\leq
x)/n$, where $I(X_{i}\leq x)$ is the indicator function, equal to 1 if
$X_{i}\leq x$, and equal to 0 otherwise. However, the ECDF is discontinuous as
it jumps with size $1/n$ when $x=X_{i},i=1\ldots,n$ and this is inconvenient
since, often, the CDF itself is a continuous function. The information given
by the ECDF is the starting point to estimate the PDF. Widely used methods to
obtain an estimator of the PDF are based on Kernel Density Estimation (KDE)
[34, 49, 4]. KDE employs a kernel function, which serves as a “base shape", to
estimate the density of the data distribution. There are various methods and
kernels available for KDE, each with its own characteristics. The choice of
kernel and method can significantly influence the accuracy of the density
estimation. Some common kernels include the Gaussian kernel, the rectangular
kernel, the Epanechnikov kernel, and others [40]. The most used kernel
function is the uniform kernel:

$K(x)=\begin{cases}1&\text{if }x\in[-1/2,1/2],\\\
0&\text{otherwise.}\end{cases}$

and the resulting estimation of the density is called naive kernel
$\hat{f}_{K}(x)$, [25]:

$\hat{f}_{K}(x)=\frac{1}{nh}\sum_{i=1}^{n}K\left(\frac{X_{i}-x}{h}\right),$
(1)

where $h$ is the bandwidth parameter and its value affects the width of the
KDE curves, and consequently the accuracy of the estimation.

We propose to estimate the CDF by applying the B-spline Hermite quasi-
interpolant [30] to compute $\hat{f}(x)$, approximation of the probability
density $f(x)=F^{\prime}(x)$, so that the final approximation $\hat{F}(x)$ of
the CDF will be given by integrating $\hat{f}(x)$.

Quasi-interpolation is a technique that allows to construct a local
approximant by keeping low the computational cost, see e.g, [7, 26] and
references therein. Generally, the common way to express a univariate spline
quasi-interpolant (QI) of $d$-degree reads as,

$\displaystyle Q_{d}\;f(\cdot)=\sum_{j=-d}^{N-1}\lambda_{j}(f)B_{j,d}(\cdot),$
(2)

where $B_{j,d}$ are $d$-degree B-splines assumed to be defined on an extended
knot vector
$\tau:=\\{\tau_{-d},\ldots,\tau_{N+d}\\},\,\tau_{j}\leq\tau_{j+1},$ and
spanning the space, $\mathbb{S}^{\pi}_{d}:=\langle
B_{-d,d},\ldots,B_{N-1,d}\rangle.$ The local linear functionals $\lambda_{j}$
in (2) can be computed by using several methodologies, such as differential,
integral methods, and discrete approaches see,e.g., [8, 24, 26, 41]. The main
advantage of QI is that they have a direct construction without solving big
linear systems. Moreover, it is local, in the sense that the value of Q$f(x)$
depends only on values of $f$ in a neighbourhood of $x$.

Given an interval $[a,b]$ such that $X_{i}\in[a,b]\;\mbox{ for }i=1,\ldots,n$
and a uniform mesh $\pi=\\{a=x_{0},x_{1},\dots,x_{N}=b\\}$ defined by a
constant stepsize $h=(b-a)/N$. Note that the choice of $h$ is important and
depends on $n$ as it plays the bandwidth role for the kernel density. For the
estimation approximation of the PDF we use the B-spline Hermite quasi-
interpolant BSHQI defined in [30]. BSHQI computes the $\lambda_{j}(f)$ as
linear combination of the function $f$ and its derivatives evaluated at the
mesh points.

We define the BSHQI with uniform knot vector $\pi$, coincident auxiliary knots
and $d=2$. Hence, in the following to ease the notation, $B_{j,d}=B_{j}$. A
discrete approximation of the sought CDF is expressed as,

$F_{{h}}(x_{j}):=\frac{1}{n}\sum_{i=1}^{n}I(X_{i}\leq x_{j}),\qquad
j=0,\ldots,N$ (3)

Starting from $F_{h}(x)$ it is possible to approximate $f(x)$ at the mesh
points by computing the first derivative using finite differences:

$\begin{array}[]{l}f(x_{j})=F^{\prime}(x_{j})\approx
F^{\prime}_{h,j}=\displaystyle{\frac{F_{h}(x_{j+1})-F_{h}(x_{j-1})}{2h}},\qquad
j=1,\ldots,N-1\\\\[8.5359pt] f(x_{0})=F^{\prime}(x_{0})\approx
F^{\prime}_{h,0}=\displaystyle{\frac{F_{h}(x_{1})-F_{h}(x_{0})}{h}},\qquad
f(x_{N})=F^{\prime}(x_{N})\approx
F^{\prime}_{h,N}=\displaystyle{\frac{F_{h}(x_{N})-F_{h}(x_{N-1})}{h}.}\end{array}$
(4)

Since the used quasi-interpolant is of Hermite type we also need an
approximation of the first derivative at the same mesh points:

$\begin{array}[]{l}f^{\prime}(x_{j})=F^{\prime\prime}(x_{j})\approx
F^{\prime\prime}_{h,j}=\displaystyle{\frac{F_{h}(x_{j+1})-2F_{h}(x_{j})+F_{h}(x_{j-1})}{h^{2}}},\qquad
j=1,\ldots,N-1\\\\[8.5359pt] f^{\prime}(x_{0})=F^{\prime\prime}(x_{0})\approx
F^{\prime\prime}_{h,0}=0,\qquad
f^{\prime}(x_{N})=F^{\prime\prime}(x_{N})\approx
F^{\prime\prime}_{h,N}=0.\end{array}$ (5)

Note that the values attained by the density at $x_{0}$ and $x_{N}$ has been
set to be zero, as this is the expected value of the sought CDF.

From the definition of the coefficients of BSHQI for $d=2$, we get

$\begin{array}[]{l}\lambda_{j}=\frac{1}{2}\left(F^{\prime}_{h,(j+1)}+F^{\prime}_{h,(j+2)}\right)-\frac{1}{4}h\left(-F^{\prime\prime}_{h,(j+1)}+F^{\prime\prime}_{h,(j+2)}\right),\qquad
j=-1,\ldots,N-2\\\\[8.5359pt]
\lambda_{-2}=F^{\prime}_{h,0},\quad\lambda_{N-1}=F^{\prime}_{h,N}.\\\
\end{array}$ (6)

The following theorem proves that the constructed $\hat{f}$ is indeed a
continuous density function having first derivative continuous as well.

######  Theorem 2.1.

The function $\hat{f}$, BSHQI estimation of $f$ in a given interval $[a,b]$,

$\hat{f}(\cdot)=\sum_{j=-d}^{N-1}\lambda_{j}(f)B_{j}(\cdot),$ (7)

with $\lambda_{j}$ as defined in (6), is a density function. In particular:

  * a)

$\hat{f}(\cdot)\geq 0$,

  * b)

$\int_{-\infty}^{+\infty}\hat{f}(x)\,dx=1$,

  * c)

$\hat{f}\in C^{1}[a,b]$.

Proof.

To prove a) it is sufficient to show that $\lambda_{j}\geq 0$, for
$j=-d,\ldots,N-1$. In particular, since $h$ is constant, setting

$F_{j}:=F_{h}(x_{j}),$

it can be shown that

$\lambda_{j}=\frac{F_{(j+2)}-F_{(j+1)}}{h},\qquad j=-1,\ldots,N-2$ (8)
$\lambda_{-2}=\lambda_{-1}\qquad\lambda_{N-1}=\lambda_{N-2}.$

by substituting (4) and (5) into the equation (6). Therefore, it is
straightforward to see that they are always positive.

To prove b), we set $\hat{f}(x)=0$ outside the interval $[a,b]$. Then, knowing
that the integral of a B-spline is given by

$\int_{a}^{b}B_{i,d}(x)\,dx=\int_{\tau_{i}}^{\tau_{i+d+1}}B_{i,d}(x)\,dx=\frac{\tau_{i+d+1}-\tau_{i}}{d+1},$

recalling that $d=2$, we have

$\tau_{i+d+1}-\tau_{i}=\left\\{\begin{array}[]{lll}3h&\mbox{for
}&i=0,\ldots,N-3,\\\ 2h&\mbox{for }&i=-2,N-2,\\\ h&\mbox{for }&i=-1,N-1.\\\
\end{array}\right.$

Therefore,

$\int_{-\infty}^{+\infty}\hat{f}(x)\,dx=\int_{a}^{b}\sum_{j=-d}^{N-1}\lambda_{j}(f)B_{j}(x)\,dx=\sum_{j=-d}^{N-1}\lambda_{j}(f)\int_{\tau_{j}}^{\tau_{j+d+1}}B_{j}(x)\,dx=$
$\begin{split}&=\frac{F_{1}-F_{0}}{h}\frac{h}{3}+\frac{F_{1}-F_{0}}{h}\frac{2h}{3}+\frac{F_{2}-F_{1}}{h}h+\cdots+\frac{F_{N-3}-F_{N-2}}{h}h+\frac{F_{N-1}-F_{N}}{h}\frac{2h}{3}+\frac{F_{N-1}-F_{N}}{h}\frac{h}{3}\\\
&=-F_{0}+F_{N}=0+1=1\end{split}$

The point c) descends from the properties of the B-spline functions of degree
$2$.

∎

In the following we investigate the consistency of the derived density
function following the analysis proposed in [14]. In particular, the next
Lemma will be useful as a preliminary result.

######  Lemma 2.2.

The coefficients in (8) can be obtained by evaluating $\hat{f}_{K}(x)$ at
$c_{j}:=x_{j+1}+h/2,j=-1,\ldots,N-2$, i.e.:

$\lambda_{j}=\hat{f}_{K}(c_{j})=\frac{1}{nh}\sum_{i=1}^{n}K\left(\frac{X_{i}-c_{j}}{h}\right),\qquad
j=-1,\ldots,N-2$

and

$\lambda_{-2}=\lambda_{-1},\qquad\lambda_{N-1}=\lambda_{N-2},$

######  Theorem 2.3.

Let $X_{1},\ldots,X_{n}$ denote i.i.d. observations having a PDF $f(x)\in
C^{1}[a,b]$, $f$ and $f^{\prime}$ bounded, and let $B_{j}(x)$ denote the
$j-th$, $2$-nd degree B-spline basis. Let $x\in[a,b]$ and let $\hat{f_{n}}(x)$
be as in (7), with uniform mesh defined in $[a,b]$ choosing a constant $h$
such that as $n\to\infty$, as $nh\to\infty$ and as $h\to 0$, then
$\hat{f}_{n}$ is a uniformly consistent estimator of $f$.

###### Proof.

To prove the pointwise consistency of $\hat{f}_{n}$ it is necessary to show
that the $MSE(\hat{f}_{n}(x))\rightarrow 0$ as $n\to\infty$, as $nh\to\infty$
and as $h\to 0$. In the following we use the well-known "Bias-Variance"
formulation:

$MSE(\hat{f}_{n}(x))\equiv
E\left[\left|\hat{f}_{n}(x)-f(x)\right|^{2}\right]=\text{Var}(\hat{f}_{n}(x))+\left[E(\hat{f}_{n}(x))-f(x)\right]^{2}\equiv\text{Var}(\hat{f}_{n}(x))+\text{Bias}^{2}(\hat{f}_{n}(x)).$

We start considering the absolute value of the Bias:

$\begin{split}\left|\text{Bias}(\hat{f}_{n}(x))\right|&=\left|E\left(\frac{1}{nh}\sum_{j}\lambda_{j}B_{j}(x)\right)-f(x)\right|\\\
&=\left|E\left(\frac{1}{nh}\sum_{j}\sum_{i=1}^{n}K\left(\frac{X_{i}-c_{j}}{h}\right)B_{j}(x)\right)-f(x)\right|\\\
&=\left|\frac{1}{h}\sum_{j}B_{j}(x)E\left(K\left(\frac{X-c_{j}}{h}\right)\right)-f(x)\right|\\\
&=\left|\frac{1}{h}\sum_{j}B_{j}(x)\left(\int
f(X)K\left(\frac{X-c_{j}}{h}\right)dX\right)-f(x)\left(\int\frac{1}{h}\sum_{j}K\left(\frac{X-c_{j}}{h}\right)B_{j}(x)dX\right)\right|\end{split}$

since
$\displaystyle{\int\frac{1}{h}\sum_{j}K\left(\frac{X-c_{j}}{h}\right)B_{j}(x)dX=1}$
and denoting with  $I_{x,h}=[x_{j+1},x_{j+2}]=supp\\{B_{j}(x)\\}\cap
supp\\{K\left(\frac{X-c_{j}}{h}\right)\\}$ we have,

$\begin{split}&=\left|\frac{1}{h}\sum_{j}B_{j}(x)\left(\int
f(X)K\left(\frac{X-c_{j}}{h}\right)dX\right)-f(x)\frac{1}{h}\sum_{j}B_{j}(x)\int
K\left(\frac{X-c_{j}}{h}\right)dX\right|\\\
&=\left|\frac{1}{h}\sum_{j}B_{j}(x)\left(\int(f(X)-f(x))K\left(\frac{X-c_{j}}{h}\right)dX\right)\right|\\\
&\leq\sup_{X\in I_{x,h}}|(f(X)-f(x)|\frac{1}{h}\sum_{j}B_{j}(x)\int
K\left(\frac{X-c_{j}}{h}\right)dX=\sup_{X\in I_{x,h}}|f(X)-f(x)|\\\
&\leq\sup_{\xi\in I_{x,h}}|f^{\prime}(\xi)|h\\\ \end{split}$

therefore

$\text{Bias}^{2}\left(\hat{f}_{n}(x)\right)\leq h^{2}\left(\sup_{\xi\in
I_{x,h}}|f^{\prime}(\xi)|\right)^{2}.$

Considering now the variance we deduce that:

$\begin{split}\text{Var}(\hat{f}_{n}(x))&=\text{Var}\left(\sum_{j}\lambda_{j}B_{j}(x)\right)=\\\
&=\text{Var}\left(\frac{1}{nh}\sum_{j}\sum_{i=1}^{n}K\left(\frac{X_{i}-c_{j}}{h}\right)B_{j}(x)\right)=\\\
&=\mathbf{E}\left[\left(\frac{1}{nh}\sum_{j}\sum_{i=1}^{n}K\left(\frac{X_{i}-c_{j}}{h}\right)B_{j}(x)\right)^{2}\right]-\left(\mathbf{E}\left[\frac{1}{nh}\sum_{j}\sum_{i=1}^{n}K\left(\frac{X_{i}-c_{j}}{h}\right)B_{j}(x)\right]\right)^{2}\\\
&=\frac{1}{nh^{2}}\sum_{j}\sum_{z}B_{j}(x)B_{z}(x)\int
K\left(\frac{X-c_{j}}{h}\right)K\left(\frac{X-c_{z}}{h}\right)f(X)dX-\\\
&\frac{1}{nh^{2}}\sum_{j}\sum_{z}B_{j}(x)B_{z}(x)\int
K\left(\frac{X-c_{j}}{h}\right)f(X)dX\int
K\left(\frac{X-c_{z}}{h}\right)f(X)dX\\\
&=\frac{1}{nh^{2}}\sum_{j}\sum_{z}B_{j}(x)B_{z}(x)\\\ &\left(\int
K\left(\frac{X-c_{j}}{h}\right)K\left(\frac{X-c_{z}}{h}\right)f(X)dX-\int
K\left(\frac{X-c_{j}}{h}\right)f(X)dX\int
K\left(\frac{X-c_{z}}{h}\right)f(X)dX\right)\\\ \end{split}$

since the kernel $K$ is evaluated at the mesh points it is easy to prove that
this quantity is bounded and

$h=\int K\left(\frac{X-c_{j}}{h}\right)K\left(\frac{X-c_{z}}{h}\right)dX,$

hence,

$\begin{split}\text{Var}(\hat{f}_{n}(x))=&\leq\frac{1}{nh^{2}}\sum_{j}\sum_{z}B_{j}(x)B_{z}(x)\left(\sup_{X\in\mathbb{R}}|f(X)|h+\sup_{X\in\mathbb{R}}(f(X))^{2}\left(\int
K\left(\frac{X-c_{j}}{h}\right)dX\right)\left(\int
K\left(\frac{X-c_{z}}{h}\right)dX\right)\right)\leq\\\
&\leq\frac{1}{nh^{2}}\sum_{j}\sum_{z}B_{j}(x)B_{z}(x)\left(\sup_{X\in\mathbb{R}}|f(X)|h+\sup_{X\in\mathbb{R}}(f(X))^{2}h^{2}\right)\\\
&=\frac{1}{nh}\left(\sup_{X\in\mathbb{R}}|f(X)|+\sup_{X\in\mathbb{R}}(f(X))^{2}h\right)\left(\sum_{j}B_{j}(x)\right)\left(\sum_{z}B_{z}(x)\right)\\\
&=\frac{1}{nh}\sup_{X\in\mathbb{R}}|f(X)|+\frac{1}{n}\sup_{X\in\mathbb{R}}(f(X))^{2}.\end{split}$

∎

Then we have

$\text{MSE}(\hat{f}_{n}(x))=\text{Bias}^{2}(\hat{f}_{n}(x))+\text{Var}(\hat{f}_{n}(x))\leq
h^{2}\left(\sup_{\xi\in
I_{x,h}}|f^{\prime}(\xi)|\right)^{2}+\frac{1}{nh}\sup_{X\in\mathbb{R}}|f(X)|+\frac{1}{n}\sup_{X\in\mathbb{R}}(f(X))^{2}.$
(9)

This quantity will tend to zero as as $n\rightarrow\infty$, as
$nh\rightarrow\infty$, and as $h\rightarrow 0$.

Since $\displaystyle{\sup_{\xi\in I_{x,h}}|f^{\prime}(\xi)|\leq
S_{1}}:=\displaystyle{\sup_{\xi\in[a,b]}|f^{\prime}(\xi)|}$, denoting by
$\displaystyle{S_{0}:=\sup_{X\in[a,b]}|f(X)|}$, then, the upper bound for the
MSE can be written as:

$MSE(\hat{f}_{n}(x))\leq\text{Bias}^{2}(\hat{f}_{n}(x))+\text{Var}(\hat{f}_{n}(x))\leq
h^{2}(S_{1})^{2}+\frac{1}{nh}S_{0}+\frac{1}{n}S_{0}^{2}.$

The above upper bound does not depend on $x$ as so the uniformly consistency
is proved. ∎

In the BSHQI density estimation, the bandwidth $h$ can be freely chosen as
long as the assumption of Theorem 2.3 are satisfied. Thus, the optimal
bandwidth $h$ can be chosen by minimizing the MSE neglecting the smallest term
$\frac{1}{n}S_{0}^{2}$:

$h_{\text{opt}}(x):=\left(\frac{1}{2n}\frac{S_{0}}{(S_{1})^{2}}\right)^{\frac{1}{3}}\sim
n^{-1/3}.$ (10)

This choice for the smoothing bandwidth leads to a MSE at the rate

$\text{MSE}_{\text{opt}}(\hat{f}_{n}(x))=O(n^{-\frac{2}{3}}).$ (11)

In the previous analysis, our focus was solely on a single point, $x$.
However, in a broader context, our goal is to manage the overall MSE for every
point. In such cases, a straightforward extension is the mean integrated
square error (MISE) of $\hat{f}_{n}(x)$. We have the following corollary:

######  Corollary 2.4.

Let $f$ be a probability density function on $\mathbb{R}$, and let
$A\subseteq\mathbb{R}$ be an open region with $A=\\{x\mid f(x)\neq 0\\}$,
where $f\in C^{1}(\mathbb{R})$, and both $f$ and $f^{\prime}$ are bounded. If
$A$ is contained in a closed and bounded region, then, with $\hat{f}_{n}$ as
defined above, if $h\to 0$ and $nh\to\infty$ as $n\to\infty$, then

$\text{MISE}(\hat{f}_{n})=\int\text{MSE}(\hat{f}_{n}(x))\,dx=\int\mathbf{E}\left(\hat{f}_{n}(x)-f(x)\right)^{2}dx\to
0,$

as $n\to\infty$, i.e., $\hat{f}_{n}$ is a consistent estimator of $f$ in mean
integrated squared error (MISE).

The above evidence concludes the investigation of the consistency results of
the B-spline estimator with the proposed approach. However, in this context,
how to choose $h$ is an unsolved problem in statistics known as bandwidth
selection. Most bandwidth selection approaches either suggest an estimate of
AMISE and then aim to minimize the estimated AMISE. For more details, we refer
to [37].

######  Observation 2.5.

Given the weighted CDF

$F_{{h,w}}(x_{j})={\frac{1}{\sum_{i=1}^{n}w_{i}}}\sum_{i=1}^{n}w_{i}I(X_{i}\leq
x_{j}),\qquad j=0,\ldots,N,$ (12)

with $w_{i}$ a weight associated to the corresponding observation $X_{i}$,
then the results of Theorem 2.1 and Theorem 2.3 continue to hold.

######  Observation 2.6.

Note that the continuous CDF $\hat{F}$ is then computed by integrating the
density $\hat{f}$ in Equation (7).

###  2.1 Statistical Tests for marginals fitting with BSHQI spline

In this subsection we compare the BSHQI density estimation with the classical
empirical approach constructed by using the Gaussian Kernel Density function.
The careful selection of the number of bins is of crucial importance in
approximating density, as it directly influences the accuracy and the visual
representation of the data distribution pattern. For both procedures
considered in this work, it is possible to select different criteria for
choosing the number of bins, in particular we consider, the so called Rice’s
Rule [45], where the number of bins is equal to $2\times\lceil n^{1/3}\rceil$.

There is no single optimal criterion for selecting the most suitable bins. For
the experiments conducted in this work, unless otherwise indicated, we will
use the Rice rule. To assess the goodness of the produced model, we conduct
two statistical tests: the Kolmogorov-Smirnov (KS) Test [29] and the Cramér-
von Mises (CvM) Test [2]. Additionally, we show the error in terms of Average
Mean Integrated Squared Error (AMISE) and Root Mean Square Error (RMSE) for
the computed probability density functions.

We perform the tests on three different distributions: a normal distribution
$X\sim\mathcal{N}(\mu,\sigma^{2})$ with $\mu=5$ and $\sigma^{2}=0.3$, an
exponential distribution $X\sim\text{Exp}(\lambda)$ with $\lambda=1$ and a
third distribution consisting of a mixture of Gaussians with different means
and variances.

(a) Pdf

(b) Cdf

Figure 1: Comparison of samples generated from $X\sim\mathcal{N}(5,0.3)$ with
the KDEpy and BSHQI method for probability density (a) and for the cumulative
distribution (b).

AMISE RMSE KS-Test Cramér–von Mises statistic p-value statistic p-value BSHQI
3.43e-06 1.13e-04 7.75e-03 2.77e-01 2.05e-01 2.58e-01 KDEpy 1.08e-05 3.53e-04
1.21e-02 1.58e-02 9.08e-01 4.05e-03

Table 1: Statistics ran on the results for Normal Density Estimation— Rice’s
Rule for bins.

All the numerical experiments are performed using Python 3.10 on a computing
system equipped with Windows 11 operating system, 16 GB of RAM, and powered by
an Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz with a base clock speed of 2.59
GHz. The python package
KDEpy111https://kdepy.readthedocs.io/en/latest/introduction.html has been
chosen for the comparison since, in our opinion, it is the most efficient
among the available Python routines for this task; the Kernel density
estimation is constructed following the theory in [23].

(a) Pdf

(b) Cdf

Figure 2: Comparison of samples generated with the KDE and BSHQI method for
probability density (a) and for the cumulative distribution (b)

AMISE RMSE KS-Test Cramér–von Mises statistic p-value statistic p-value
EMP_BSHQI 2.27e-06 2.96e-05 7.78e-03 2.73e-01 1.18e-01 5.02e-01 EMP_KDEpy
1.68e-04 2.19e-03 6.96e-02 0 5.79e+01 2.01e-08

Table 2: Statistics ran on the results for Exponential Density Estimation—
Rice’s Rule for the bins.

For all the considered distributions, we generated a group of $n=2^{15}$
samplings for $20$ iterations. This iterative process allows us to calculate
the AMISE, the RMSE, the values of both statistics and relative p-values, and
we derive their average values as comprehensive measures of performance.
Furthermore, we evaluated efficiency in terms of computational time by
calculating the mean and standard deviation after the set number of
iterations. Regarding the statistical tests, the null hypothesis states that
the true underlying distribution and the empirical one are identical; the
alternative hypothesis suggests that they are not. The used statistics is the
maximum absolute difference between the exact values and the ones computed by
the empirical distribution functions at the same samples. If the KS or CvM
statistics are large, then the p-value will be small, and this may be taken as
evidence against the null hypothesis in favor of the alternative.

By observing the Tables 1, 2, 3 it can be seen that the density approximation
with BSHQI is preferable to the classical empirical evaluation of the
distribution taken into account.

(a) Pdf

(b) Cdf

Figure 3: Comparison of samples generated with the KDE and BSHQI method for
probability density (a) and for the cumulative distribution (b).

AMISE RMSE KS-Test Cramér–von Mises statistic p-value statistic p-value
EMP_BSHQI 1.84e-06 2.16e-05 4.91e-03 8.22e-01 8.09e-02 6.87e-01 EMP_KDEpy
9.35e-06 1.10e-04 1.06e-02 5.04e-02 5.18e-01 3.58e-02

Table 3: Statistics ran on the results for Mixture Gaussian Density
Estimation— Rice’s Rule for the bins.

Indeed, comparing the p-values of both tests in each of the Tables 1, 2, 3,
allows us to accept the hypothesis when using the BSHQI, contrary to what can
be concluded when referring to KDEpy.

Moreover, for all three experiments, the AMISE and RMSE obtained with the
proposed approach are lower compared to the AMISE and RMSE obtained using
KDEpy. This can be observed in Figures 1, 3, and 2, where the estimates of the
PDF and the CDF for both methods are compared to the real distribution. Note
that different results could be obtained by choosing a different rule for the
number of bins. In particular, for the normal distribution in Figure 1 we can
observe a good agreement for both the methodologies; for the exponential and
the mixture ones in Figures 3, and 2, the result produced by BSHQI seems to
better alienate with the original distribution while KDE seems to be less
accurate as it always produce an under-estimate of the original distribution.
The visual output is also supported by the results shown in Tables 1,2,3: if
on the one hand the reached accuracy seems similar in terms of AMISE and RMSE,
on the other hand, the statistics conducted to test the validity of the null
hypothesis clears out any doubt in assessing a better goodness of fit, under
statistical point of view, of the BSHQI method.

Algorithm | Mean Time (ms) | Standard Deviation (ms)
---|---|---
KDEpy | $0.0241$ | $\pm 0.0011$
BSHQI | $0.0124$ | $\pm 0.0009$
Table 4: Comparison of Execution Times for KDE and BSQH Algorithms

In addition examining the results in Table 4, it can be seen an advantage in
terms of computational time when utilizing the BSHQI, compared to the
alternative approach. While the difference is not substantial, it is
sufficient to highlight the efficiency of our approach as the time is almost
halved when using BSHQI.

##  3 Copulas Mixture Model

In this section, we leverage the density estimation strategy introduced in
Section 2 in the context of copulas. Our objective is to introduce the concept
of copula and emphasize its profound connection to marginal distributions,
with a specific focus on its application in formulating a novel clustering
algorithm based on copula mixture. Copulas functions are a useful tool
employed to easily express multivariate distributions by specifying the
marginals. In this Section the principal concepts are revised following the
setting in [11].

######  Definition 3.1.

A D-dimensional copula is a CDF with uniform marginals:

$C:[0,1]^{D}\rightarrow[0,1]\;\mbox{such that }C(u)=C(u_{1},\ldots,u_{D}).$

######  Theorem 3.2.

(Sklar’s Theorem) Consider a D-dimensional CDF, $G$, with marginals
$F_{1},\ldots,F_{D}$. Then there exist a copula, $C$, such that

$G(x_{1},\ldots,x_{D})=C(F_{1}(x_{1}),\ldots,F_{D}(x_{D}))$ (13)

for all $x_{i}\in[-\infty,\infty]$ and $i=1,\ldots,D$. If $F_{i}$ is
continuous for all $i=1,\ldots,D$ then $C$ is unique. In the opposite
direction, given a copula $C$ and univariate CDFs, $F_{1},\ldots,F_{D}$, then
$G$ as in (13) is a multivariate CDF with marginals $F_{1},\ldots,F_{D}.$

Given a CDF $G$ with PDF $g$, and a copula $C$ defined as in Definition 3.1,
the density copula function $c$ can be computed as,

$c(u_{1},\ldots,u_{D})=\frac{g(F_{1}^{-1}(u_{1}),\ldots,F_{D}^{-1}(u_{D}))}{f_{1}(F_{1}^{-1}(u_{1}))\cdots
f_{D}(F_{D}^{-1}(u_{D}))},$ (14)

where $f_{1},\ldots,f_{D}$ are the PDFs of the marginals.

A novel algorithm, which takes advantage of the proposed construction for
empirical cumulative distribution based on BSHQI, for estimating the marginals
of a chosen copula is presented. As main application, we show the performance
of such algorithm in the clustering context. Therefore, the following
formulation will be framed within the clustering setting.

Our goal is to implement a model-based algorithm capable of correctly
identifying how the instances of the dataset can be grouped into different
clusters. In this sense, using copulas provides a way to fit data that have
different probability distributions, thus having an advantage in better
discriminating the possible clusters of a dataset. We can assume an a priori
model made with $K$ clusters and the data in each single cluster are
distributed as a multidimensional Copula belonging to the Elliptical family,
in particular Gaussian Copula, and Archimedean family, in which we consider
the Clayton, Gumbel and Frank copulas [19]. We therefore stress the fact that
our mixture distribution is composed of a linear combination of Copulas. This
linear combination is called “Copula Mixture Model". In the following we
describe in details the derived formulation.

######  Definition 3.3 (Semiparametric approach).

A Copula Mixture is a function consisting of several Copula density functions
$c_{k}$, with $k\in\\{1,\ldots,K\\}$ and $K$ denoting the number of clusters
of the considered dataset. Each Copula $c_{k}$ in the mixture is characterized
by a vector $\boldsymbol{\omega}$ that defines the parameters of the specific
copula chosen for the mixture, and by the methods chosen for the approximation
of the marginals. Moreover, for each Copula density function $c_{k}$ is
defined a mixing probability $\pi_{k}$, referred to as mixing coefficient,
such that:

$\sum_{k=1}^{K}\pi_{k}=1.$

Let us assume a dataset $\mathbf{X}=(X_{1},X_{2},\ldots,X_{D})$ where each
$X_{i}$ consists of $n$ i.i.d. observations, and where the $i$-th observation
is $\mathbf{x}_{i}=(x_{i,1},x_{i,2},\ldots,x_{i,D})^{T}$ for $i=1\ldots,n$.

We express the probabilistic model of mixture copulas for all observations in
the following form:

$p(\mathbf{X}|\boldsymbol{\theta})=\prod_{i=1}^{n}p(\mathbf{x}_{i})=\prod_{i=1}^{n}\sum_{k=1}^{K}\pi_{k}g_{k}(\mathbf{x}_{i}|\boldsymbol{\omega}_{k})$
(15)

where $p$ is the probability,
$\boldsymbol{\theta}:=\\{\pi_{k},\boldsymbol{\omega}_{k}\\}$ indicates the
parameters of the model: $\pi_{k}$ is the mixing probability and
$\boldsymbol{\omega_{k}}$ represents the vector of parameters with respect to
the chosen copula, while $g_{k}({\bf x}_{i}|\boldsymbol{\omega_{k}})$
represents the multivariate distribution constructed through the copula
density $c_{k}$, explicitly:

$\begin{split}g_{k}({\bf
x}_{i}|\boldsymbol{\omega}_{k})&=c_{k}(F_{1}(x_{i,1}),\ldots,F_{D}(x_{i,D})|\boldsymbol{\omega}_{k})\bigg{(}f_{1}(x_{i,1})\times\cdots\times
f_{D}(x_{i,D})\bigg{)}\\\
&=c_{k}(F_{1}(x_{i,1}),\ldots,F_{D}(x_{i,D})|\boldsymbol{\omega}_{k})\prod_{j=1}^{D}f_{j}(x_{i,j}).\end{split}$

Since there are different families of copulas, the parameter
$\boldsymbol{\omega}_{k}$ is the one related to the specific copula that is
chosen to model the cluster $k$. The goal of the mixture model is to find the
optimal parameters in $\boldsymbol{\theta}$, in the equation (15), that
maximize the log-likelihood $\mathcal{L}(\mathbf{X}|\boldsymbol{\theta})$:

$\underset{\boldsymbol{\theta}}{\operatorname{arg\,max}}\mathcal{L}(\mathbf{X}|\boldsymbol{\theta})=\underset{\boldsymbol{\theta}}{\operatorname{arg\,max}}\log
p(\mathbf{X}|\boldsymbol{\theta})=\underset{\pi_{k},\boldsymbol{\omega}_{k}}{\operatorname{arg\,max}}\sum_{i=1}^{n}\log\sum_{k=1}^{K}\pi_{k}g_{k}({\bf
x}_{i}|\boldsymbol{\omega}_{k}).$ (16)

Usually, in order to solve the optimization problem in (16), the Expectation-
Maximization (EM) [9] algorithm is employed together with mixture models. To
derive the probability that an observation $\mathbf{x}_{i}$, is drawn from
$g_{k}$, we introduce the latent variable
$\mathbf{z}_{i}=(z_{i1},\cdots,z_{iK})^{T}$ such that $z_{ik}\in\\{0,1\\}$
with $k\in\\{1,\cdots,K\\}$. In this context the introduction of latent
variables is a common practice [5] that enhances both the theoretical
framework and the rationale behind employing the expectation-maximization
algorithm. We define the joint distribution $p({\bf x}_{i},{\bf z}_{i})$ in
terms of a marginal distribution $p({\bf z}_{i})$ and a conditional
distribution $p({\bf x}_{i}|{\bf z}_{i})$,

$p({\bf x}_{i},\mathbf{z}_{i}):=p({\bf
x}_{i}|\mathbf{z}_{i})p(\mathbf{z}_{i}).$ (17)

The marginal distribution for ${\bf z}_{i}$ is characterized by the mixing
coefficients, with the specification that $p(z_{ik}=1)=\pi_{k}$.

We know beforehand that each $z_{ik}$ occurs independently and that it can
only take the value one when $k$ is equal to the cluster from which the
observation comes from, then the overall probability is:

$p(\mathbf{z}_{i})=p(z_{i1}=1)^{z_{i1}}p(z_{i2}=1)^{z_{i2}}\cdots
p(z_{iK}=1)^{z_{iK}}=\prod_{k=1}^{K}\pi_{k}^{z_{ik}},$

while the conditional distribution $p({\bf x}_{i}|{\bf z}_{i})$ can be written
as,

$p(\mathbf{x}_{i}|\mathbf{z}_{i})=\prod_{k=1}^{K}g_{k}(\mathbf{x}_{i}|\boldsymbol{\bf\omega}_{k})^{z_{ik}}.$
(18)

Let us introduce ${\bf Z}$ as the matrix whose $i^{th}$ row is the vector of
the latent variables $\mathbf{z}_{i}$, then we have the overall joint
distribution,

$p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})=\prod_{i=1}^{n}\prod_{k=1}^{K}\pi_{k}^{z_{i_{k}}}g_{k}({\bf
x}_{i}|\boldsymbol{\omega}_{k})^{z_{i_{k}}}$

and hence the log-likelihood is,

$\log(p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}))=\sum_{i=1}^{n}\sum_{k=1}^{K}z_{i_{k}}\left(\log\pi_{k}+\log
g_{k}({\bf x}_{i}|\boldsymbol{\omega}_{k})\right).$ (19)

The formulated expression for the joint distribution
$p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})$ leads to significant
simplifications in the EM algorithm. However this cannot be computed, since
${\bf z}_{i}$ is unknown, then we evaluate the conditional probability of
${\bf z}_{i}$ given ${\bf x}_{i}$, and its value can be determined using
Bayes’ theorem,

$p(z_{ik}=1|{\bf x}_{i})=\frac{p({\bf
x}_{i}|z_{ik}=1)\,p(z_{ik}=1)}{p(\bf{x}_{i})}=\frac{p({\bf
x}_{i}|z_{ik}=1)\,p(z_{ik}=1)}{\displaystyle{\sum_{k^{\prime}=1}^{K}\pi_{k^{\prime}}p({\bf
x}_{i}|z_{i{k^{\prime}}}=1)p(z_{i{k^{\prime}}}=1)}},$

where the quantity $p({\bf x}_{i})$ has been computed via marginalization.
Knowing that $p(z_{ik}=1)=\pi_{k}$ and $p({\bf x}_{i}|z_{ik}=1)=g_{k}({\bf
x}_{i}|\boldsymbol{\omega}_{k})$, then, the above equation becomes:

$p(z_{ik}=1|{\bf x}_{i})=\frac{\pi_{k}g_{k}({\bf
x}_{i}|\boldsymbol{\omega}_{k})}{\displaystyle{\sum_{k^{\prime}=1}^{K}\pi_{k^{\prime}}g_{k^{\prime}}({\bf
x}_{i}|\boldsymbol{\omega}_{k^{\prime}})}}.$ (20)

The quantity in (20) is called responsibility of $k$-th cluster to observation
$i$ and from now on will be denoted with the symbol $\gamma_{ik}$. Then for
every cluster we have an array of responsabilities. This quantity is crucial
in the Expectation step of the EM algorithm for the maximization of the
complete log-likelihood.

##  4 Expectation-Maximization for Copula Mixture Model

We formalize the expectations maximization algorithm for copula mixture model
in the general form, this can be implemented in different ways:

  * •

a completely parametric way in which the parameters of the copula and the
parameters of the marginal probability densities are estimated. In this case
there are two approaches: Inference For Marginal (IFM) [19], and
Expectation/Conditional Maximization ECM [33]. This last one, although it may
work well for small-sized data, requires high computational costs when there
is a lack of a-priori knowledge about the marginals. In such cases, one needs
to search for the distribution that fits the data well within a set of
distributions;

  * •

a semiparametric way in which, the estimation of marginals is approached
empirically. The semiparametric nature of this approach strikes a balance
between flexibility and computational efficiency.

In this work, we adopt the semiparametric approach, leveraging the
approximation properties of densities with the previously introduced BSHQI.
Below, we describe the details of the proposed algorithm.

In the Expectation step, we employ the existing parameter values
$\boldsymbol{\theta}:=\boldsymbol{\theta}^{(t)}$ to determine the posterior
distribution of the latent variables at the $t$-step of the algorithm, denoted
as $p({\bf Z}|{\bf X},\boldsymbol{\theta}^{(t)})$. Subsequently, we utilize
this posterior distribution to calculate the expectation of the complete-data
log likelihood, evaluated for a general parameter vector
$\boldsymbol{\theta}$. This expectation, is called auxiliary function
represented as $Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})$, is given by

$Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})=\mathbb{E}_{(\mathbf{Z}\mid\mathbf{X},\boldsymbol{\theta}^{(t)})}\left(\log
p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})\right).$

To simplify the notation, we shall consider
$\mathbb{E}_{(\mathbf{Z}\mid\mathbf{X},\boldsymbol{\theta}^{(t)})}\left(\log
p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})\right)=\mathbb{E}\left(\log
p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})\right)$ and so we have:

$\begin{split}Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})&=\mathbb{E}[\log
p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}^{(t)})]=\sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})\log
p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})=\sum_{i=1}^{n}\sum_{k=1}^{K}\gamma_{ik}^{(t)}\left(\log\pi_{k}+\log
g_{k}(x_{i}|\boldsymbol{\omega}_{k})\right)\\\
&=\sum_{i=1}^{n}\sum_{k=1}^{K}\gamma_{ik}^{(t)}\bigg{(}\log\pi_{k}+\log\Big{(}c_{k}(\hat{F}_{1}(x_{i,1}),\ldots,\hat{F}_{D}(x_{i,D})|\boldsymbol{\omega}_{k}\Big{)}+\sum_{j=1}^{D}\log\hat{f}_{j}(x_{ij})\bigg{)}\end{split}$
(21)

in which $\gamma_{ik}^{(t)}$ is the Responsibility introduced in equation
(20).

##### Maximization step:

In the maximization step we update the parameters in
$\boldsymbol{\theta}^{(t+1)}$ by computing:

$\boldsymbol{\theta}^{(t+1)}=\underset{\boldsymbol{\theta}}{\operatorname{arg\,max}}\
Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})$

This is the most complex step of the algorithm and in the following we address
separately the computation of the optimal parameters $\pi_{k}^{(t+1)}$ and
$\boldsymbol{\omega}_{k}^{(t+1)}$.

For $\pi_{k}^{(t+1)}$, a closed form can be derived. Note that the
maximization of the function
$Q(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(t)})$ should take into account
the restriction that $\sum_{k=1}^{K}\pi_{k}=1$. Hence, we can add a Lagrange
multiplier to (21) ,

$Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})=\sum_{i=1}^{n}\sum_{k=1}^{K}\gamma_{ik}^{(t)}\left(\log\pi_{k}+\log
g_{k}(x_{i}|\boldsymbol{\omega}_{k})\right)-\lambda\Bigg{(}\sum_{k=1}^{K}\pi_{k}-1\Bigg{)}.$

Taking the derivative of $Q$ with respect to $\pi_{k}$ and setting it equal to
zero, leads to

$\frac{\partial
Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})}{\partial\pi_{k}}=\sum_{i=1}^{n}\frac{\gamma_{ik}^{(t)}}{\pi_{k}}-\lambda=0.$
(22)

Then, by rearranging the terms and applying a summation over $k$ to both sides
of the equation, we obtain:

$\sum_{i=1}^{n}\gamma_{ik}^{(t)}=\pi_{k}\lambda\rightarrow\sum_{k=1}^{K}\sum_{i=1}^{n}\gamma_{ik}^{(t)}=\sum_{k=1}^{K}\pi_{k}\lambda$

we know that the summation of all mixing coefficients equals one. In addition,
we know that summing up the responsabilities $\gamma$ over $k$ will also give
us 1. Thus we get $\lambda=n$. Using this result, we can solve (22) for
$\pi_{k}$:

$\pi_{k}^{(t+1)}:=\frac{\sum_{i=1}^{n}\gamma_{ik}^{(t)}}{n}.$

Regarding the optimization with respect to $\boldsymbol{\omega}_{k}$,
systematic decomposition into distinct maximization steps is carried out in
order to get the optimization of the complete data log-likelihood with respect
to the parameters of the model.

One notable feature of our algorithm lies in its inherent flexibility in the
selection of copulas from the outset. When a single copula is chosen
initially, the subsequent maximization step, following the updating of
marginals, exclusively targets the parameters of the selected copula until the
log-likelihood function converges. Conversely, when the initial choice
encompasses various copulas, a comprehensive fitting process is initiated
after updating the marginals. The copula that best aligns with the updated
data is then chosen based on the maximum likelihood, enhancing the
adaptability and performance of our algorithm in diverse scenarios, formally:

  * •

Maximization first step: For each cluster, use the data ${\bf X}_{j}$,
$j\in\\{1,2,\ldots,D\\}$ to update the CDFs $\hat{F}_{j}$ according to the
equation (12) in which the weights are the _responsabilities_ and compute the
PDFs $\hat{f}_{j}$, $j\in\\{1,2,\ldots,D\\}$ with the BSHQI strategy described
in Section 2.

  * •

Maximization second step (one copula): By looking at (21), for each clusters,
we need to maximize only the following,

$\sum_{i=1}^{n}\sum_{k=1}^{K}\gamma_{ik}^{(t)}\bigg{(}\log\Big{(}c_{k}(\hat{F}_{1}(x_{i,1}),\ldots,\hat{F}_{D}(x_{i,D});\boldsymbol{\omega}_{k}\Big{)}\bigg{)},$
(23)

with respect to $\boldsymbol{\omega}_{k}$ for the optimal copula parameters.

  * •

Maximization second step (two or more copulas): If the choice of different
copulas is enabled, then, the optimization of (23) is carried out also for
each copula.

Note that, in this case for the maximization of the log-likelihood we use the
Limited-memory Broyden, Fletcher, Goldfarb, Shanno (L-BFGS-B)
method222https://docs.scipy.org/doc/scipy/reference/optimize.minimize-
lbfgsb.html, that is an optimization algorithm in the family of quasi-Newton
methods that approximates the Broyden–Fletcher–Goldfarb–Shanno algorithm
(BFGS) using a limited amount of computer memory [3, 27].

To start the algorithm we use the following initialization procedure:

##### Initialization

  * 0.

Choose a set $S$ of copulas function among Gaussian, Clayton, Gumbel, Frank.

  * 1.

Create a random clustering partition of the given data. We call $X_{j,r}$ the
subset of the $j$-th column given by the produced partition.

  * 2.

For every cluster, compute the the marginals $\hat{F}_{j}(X_{j,r})$ for
$j=1,\ldots,D$ either with the KDE methodology or with the BSHQI (the choice
is at the user discretion).

  * 3.

For every cluster, find the best copula in $S$ according to the maximum
likelihood achieved and keep their parameters $\boldsymbol{\omega}_{k}$ fitted
by the produced data
$u_{1}=\hat{F}_{1}(X_{1,r}),u_{2}=\hat{F}_{2}(X_{2,r}),\ldots,u_{D}=\hat{F}_{D}(X_{D,r})$,
and compute the associated $\pi_{k}$ for $k=1,\ldots,K$.

  * 4.

Repeat steps 1-3 for $5$ times and select the copulas partition with the
maximum likelihood so that
$\boldsymbol{\theta}^{(0)}:=\\{\pi_{k}^{(0)},\boldsymbol{\omega}_{k}^{(0)}\\}=\\{\pi_{k}^{\mbox{best}},\boldsymbol{\omega}_{k}^{\mbox{best}}\\}$
for $k=1,\ldots,K$.

##### Expectation

In this step we evaluate the responsabilities.

##### Maximization

The algorithm as just described monotonically approaches a local minimum of
the cost function. After computing the new estimates, we set
$\theta^{(t)}=(\pi_{k}^{(t)},\omega_{k}^{(t)})$ for $k=1\ldots,k$, and go to
the next Expectation step. We set the tolerance $tol\leq 10^{-4}$ than, the
best parameters are obtained when the the convergence of the log-likelihood is
reached, i.e.

$\frac{|\mathcal{L}(\mathbf{X}|\boldsymbol{\theta}^{(t+1)})-\mathcal{L}(\mathbf{X}|\boldsymbol{\theta}^{(t)})|}{1+|\mathcal{L}(\mathbf{X}|\boldsymbol{\theta}^{(t+1)})|}<tol\quad\textnormal{for}\quad
t=1,\ldots,(\emph{iter}-1).$ (24)

##  5 Experiments

In this section, we describe the experiments conducted to validate the
accuracy of our algorithm. Specifically, we commence by examining synthetic
datasets constructed with diverse families of copulas. Subsequently, we delve
into the assessment of real-world datasets. The clustering strategy using
copulas and density estimation with the proposed BSHQI technique, called
CopMixM_BSHQI, will be compared by using marginals estimation with kernel
density estimation using the Python package KDEpy, referred to as
CopMixM_KDEpy here, which is found to be the fastest among various python
packages for empirical density estimators as it utilizes the convolution Fast
Fourier Transform for calculations. Additionally, we will compare the results
obtained with clustering algorithms available from scikit-learn: Gaussian
Mixture Model (GMM), where we set the initialization to ‘random’, the
tolerance equal to $10^{-4}$ (as in our algorithm) and default parameters
setting otherwise, and K-Means where we set the initialization to random. To
measure the achieved performance we use some classical metrics, available from
scikit-learn, such as:

  * •

Silhouette Score: a value between $-1$ and $1$, with $1$ being the best;

  * •

Calinski-Harabasz Index: highest values indicate better performance.

  * •

Davies-Bouldin Score: the minimum is $0$, lower values indicate better output.

In those tests where ground truth labels are accessible, the evaluation is
done also by using the Adjusted Rand Score, Homogeneity Score, Rand Score, and
Completeness Score which are permutation invariant and the highest value $1$
indicates optimal clustering.

###  5.1 Synthetic Dataset

(a) $\mathcal{X}_{1}$

(b) $\mathcal{X}_{2}$

(c) $\mathcal{X}_{3}$

(d) $\mathcal{X}_{4}$

(e) Pairwise $\mathcal{X}_{4}$

Figure 4: Synthetic dataset: (a) Ground truth $\mathcal{X}_{1}$, (b) Ground
truth $\mathcal{X}_{2}$, (c) Ground truth $\mathcal{X}_{3}$, (d) Ground truth
$\mathcal{X}_{4}$, (e) Pairwise ground truth $\mathcal{X}_{4}$

(a) GMM

(b) CopMixM_BSHQI

(c) CopMixM_KDEpy

Figure 5: Synthetic dataset $\mathcal{X}_{1}$: (a) GMM, (b) CopMixM_BSHQI, (c)
CopMixM_KDEpy

K-Means GMM CopMixM_BSHQI CopMixM_KDEpy Silhouette Score 0.660 0.512 0.659
0.628 Calinski-Harabasz Index 5977 1570 5890 3940 Davies-Bouldin Score 0.482
0.529 0.486 0521 Adjusted Rand Score 0.974 0.437 1 0.963 Homogeneity Score
0.963 0.455 1 0.952 Rand Score 0.989 0.707 1 0.985 Completeness Score 0.964
0.797 1 0.957

Table 5: Clustering metrics for Synthetic Dataset $\mathcal{X}_{1}$. In bold
the best values.

Clusters # Points Different Copulas # Points One Copula 1 500 Gaussian 505
Gaussian 2 500 Gaussian 498 Gaussian 3 300 Gaussian 299 Gaussian 4 200 Clayton
198 Gaussian Log-likelihood -7535 -8125

Table 6: Results of CopMixM_BSHQI for the Synthetic Dataset $\mathcal{X}_{1}$
obtained with different copulas and the Gaussian copula.

(a) GMM

(b) CopMixM_BSQI

(c) CopMixM_KDEpy

Figure 6: Synthetic dataset $\mathcal{X}_{2}$: (a) GMM, (b) CopMixM_BSQI, (c)
CopMixM_KDEpy

K-Means GMM CopMixM_BSHQI CopMixM_KDEpy Silhouette Score 0.726 0.269 0.726
0.370 Calinski-Harabasz Index 11495 2640 11500 2990 Davies-Bouldin Score 0.365
1.48 0.365 0.964 Adjusted Rand Score 0.998 0.382 1 0.605 Homogeneity Score
0.994 0.488 0.999 0.738 Rand Score 0.999 0.758 1 0.845 Completeness Score
0.994 0.512 0.998 0.720

Table 7: Clustering metrics for Synthetic Dataset $\mathcal{X}_{2}$. In bold
the best values.

Clusters # Points Different Copulas # Points One_Copula 1 1013 Clayton 1152
Gaussian 2 1014 Clayton 980 Gaussian 3 985 Gumbel 858 Gaussian 4 988 Gumbel
1010 Gaussian Log-likelihood -13324 -14873

Table 8: Results of CopMixM_BSHQI for the Synthetic Dataset $\mathcal{X}_{2}$
obtained with different copulas and the Gaussian copula.

(a) GMM

(b) CopMixM_BSQI

(c) CopMixM_Emp

Figure 7: Synthetic dataset $\mathcal{X}_{3}$: (a) GMM, (b) CopMixM_BSQI, (c)
CopMixM_KDEpy

K-Means GMM CopMixM_BSHQI CopMixM_KDEpy Silhouette Score 0.638 0.253 0.638
0.384 Calinski-Harabasz Index 7874 1420 7710 702 Davies-Bouldin Score 0.474
3.8 0.476 1.47 Adjusted Rand Score 0.960 0.476 0.999 0.615 Homogeneity Score
0.932 0.518 0.997 0.598 Rand Score 0.982 0.724 1 0.804 Completeness Score
0.931 0.841 0.997 0.771

Table 9: Clustering metrics for Synthetic Dataset $\mathcal{X}_{3}$. In bold
the best values.

Clusters # Points Different Copulas # Points One Copula 1 702 Frank 698
Gaussian 2 1001 Frank 935 Gaussian 3 998 Frank 1067 Gaussian Log-likelihood
-12644 -13835

Table 10: Results of CopMixM_BSHQI for the Synthetic Dataset $\mathcal{X}_{3}$
obtained with different copulas and the Gaussian copula.

(a) GMM

(b) CopMixM_BSHQI

Figure 8: Synthetic dataset $\mathcal{X}_{4}$: (a) GMM , (b) CopMixM_BSHQI

K-Means GMM CopMixM_BSHQI CopMixM_KDEpy Silhouette Score 0.576 0.504 0.561
0.536 Calinski-Harabasz Score 2030 1420 1850 1630 Davies-Bouldin Score 0.601
0.592 0.618 0.661 Adjusted Rand Score 0.595 0.544 0.695 0.685 Homogeneity
Score 0.488 0.528 0.589 0.659 Rand Score 0.797 0.772 0.848 0.843 Completeness
Score 0.489 0.555 0.594 0.595 Miss-classification Rate: 0.12 0.11 ($\pm 0.02$)
0.08 ($\mathbf{\pm 0.01}$) 0.09 ($\pm 0.01$)

Table 11: Clustering metrics for Synthetic Dataset $\mathcal{X}_{4}$. In bold
the best values.

We examine four synthetic datasets
$\mathcal{X}_{1},\mathcal{X}_{2},\mathcal{X}_{3},\mathcal{X}_{4}$
characterized as follows:

  * •

the samples in $\mathcal{X}_{1}$ are drawn from two Claytons copula, one Frank
and one Gumbel, see Figure 4 (a);

  * •

the samples in $\mathcal{X}_{2}$ constitutes $4$ clusters and drawn from 2
Clayton copulas and 2 Gumbel copulas, see Figure 4 (b);

  * •

the samples in $\mathcal{X}_{3}$ constitutes $3$ clusters drawn from the Frank
copula, see Figure 4 (c);

  * •

the samples in $\mathcal{X}_{4}$ constitutes $2$ clusters of 3D scattered
points, see Figure 4 (d)-(e).

All the presented results are obtained with random initialization and with the
selection of different copulas. Regarding dataset $\mathcal{X}_{1}$, in Figure
5, the outcomes from GMM, CopMixM_BSHQI and clustering with copulas fitting
the marginals using the KDEpy approach are shown. Visually, the best results
are obtained with our approach, see Figure 5 (b), and KDEpy Figure 5 (c), as
they closely aligns with the ground truth in Figure 4 (a). The metrics
detailed in Table 5 show a more quantitative comparison for the four used
clustering algorithms. Notably, CopMixM_BSHQI outperforms GMM and gives better
results compared to the other methods with respect to various metrics. The
highest Silhouette Score, Calinski-Harabasz Index, and Homogeneity Score
indicate a superior cluster quality and better separation. For the Davies-
Bouldin Score, where a lower value is desirable in presence of well-defined
and compact clusters, CopMixM_BSHQI again excels. The Adjusted Rand Score,
Rand Score, and Completeness Score further support the dominance of
CopMixM_BSHQI, as it achieves perfect scores$=1.0$ in these metrics,
signifying strong agreement with the ground truth. In contrast, K-Means, GMM
and CopMixM_KDEpy exhibit lower scores, reflecting a lower level of agreement
with the true cluster assignments. Moreover, in Table 6, we present the
results of CopMixM_BSHQI on the synthetic dataset $\mathcal{X}_{1}$ under
various copula configurations. Two primary scenarios were considered: one
utilizing diverse copulas for specific clusters, and the other employing a
single Gaussian copula, which is commonly used in practice and is often
considered the default choice in mixture models. Evaluation was based on the
maximum likelihood for the choice of the cluster’s copula and the Log-
likelihood for the overall mixture model performance. The results suggest to
adopt diverse copulas tailored to specific clusters rather than using a single
one. Indeed, the higher log-likelihood highlights the superior model
adaptability, to the considered synthetic dataset $\mathcal{X}_{1}$, compared
to the simplistic use of a single Gaussian copula.

Moreover, with respect to datasets $\mathcal{X}_{2}$ and $\mathcal{X}_{3}$,
the findings depicted in Figures 6 and 7, along with the corresponding metrics
presented in Tables 8 and 9, confirm the superiority of the CopMixM_BSHQI
methodology over the other three methodologies. This validate that
CopMixM_BSHQI works better not only in effectively distinguishing distinct
clusters but also in appropriately assigning the suitable copula for each of
them. This is further illustrated in Tables 7 and 10, which show again the
better performance in the use of diverse copulas with respect to a single
copula. In these instances as well, the log-likelihood supports this
conclusion.

In the latest experiment, we show the effectiveness of the algorithm using the
synthetic dataset $\mathcal{X}_{4}$, which is in 3D and comprises of two
clusters. The results are presented in Figure 8 where only GMM and
CopMixM_BSHQI are shown as the results obtained with K-Means and CopMixM_KDEpy
are similar to the last one. Moreover in Table 11 the achieved metrics are
reported. It can be seen that the worst performance is given by GMM, while the
best sometimes is achieved by CopMixM_KDEpy, altought the discrepancy with
CopMixM_BSHQI is minimal. For this test, since there are only two classes, we
have also included accuracy in terms of miss-classification rate. As can be
observed in this scenario the CopMixM_BSHQI performance surpasses that of the
conventional K-Means, GMM and KDEpy based approach. In summary, these results
highlight the algorithm’s effectiveness in capturing complex structures within
synthetic datasets, positioning it as a robust choice for clustering
activities in similar contexts. The findings suggest that the choice of copula
significantly influences goodness of fit, with specific copula types being
favored by certain clusters. Therefore, the model has an overall robust
performance, particularly in contexts where the option to choose from multiple
copulas is available, compared to relying on a single copula.

###  5.2 Real Datasets

We conduct experiments to validate the effectiveness of the CopMixM_BSHQI
algorithm on real-world datasets. Our analysis considers three datasets: the
Australian Institute of Sport (AIS), Breast Cancer Wisconsin, and a case study
involving the clustering of textual data using the 20newsgroups
dataset333https://scikit-learn.org/0.19/datasets/twenty_newsgroups.html. For
the AIS444https://rdrr.io/cran/GLMsData/man/AIS.html and Breast Cancer
Wisconsin datasets555https://archive.ics.uci.edu/dataset/14/breast+cancer, we
compare our approach against a specific method referenced in [51], where our
initialization settings are considered the same as in [51] wherever possible.
In the case of the 20newsgroups dataset, we compare our results with the
conventional GMM.

####  5.2.1 AIS

The AIS dataset, comprising 13 measurements collected from 102 male athletes
and 100 female athletes. The primary objective of clustering this dataset is
to assess the effectiveness of our methodologies in producing results aligned
with the given ground-truth. We consider the same settings as in [51],
focusing on a subset of five variables: lean body mass (LBM), weight (Wt),
body mass index (BMI), white blood cell count (WBC), and percentage of body
fat (PBF).

This particular example exhibits a non-Gaussian distribution, as the contours
are non-elliptical and displays asymmetric dependency patterns. Therefore, the
application of copulas seems to be the optimal one for capturing this more
complex form of dependency. Our algorithm was ran with both K-Means and random
initialization. Furthermore, we explored the use of different copula choices.
We utilized the miss-classification rate metric to assess the performance of
the clustering algorithm. The results obtained with our approach, employing
both types of initialization, are presented in the Table 12 alongside with
those reported in [51].

k-means VCMM k-means CopMixM_BSHQI k-means CopMixM_BSHQI random Miss-
classification rate (Bins=$\lceil n^{1/3}\rceil$) 0.21 0.040 0.090 0.045 Miss-
classification rate (Bins=Rice) 0.21 0.040 0.040 0.035

Table 12: Results for AIS dataset. In bold the best values.

The methods called VCMM (k-means) used in [51] for this dataset, initialized
with k-means reaches the best Miss-classification value if in our model we set
the number of bins equal to $\lceil n^{1/3}\rceil$, however the CopMixM_BSHQI
with random initialization, shows a competitive result and in particular when
the Rice Rule is adopted, it achieves the best results. This discrepancy in
performance between the two initialization methods indicates that random
initialization produces more favorable results for CopMixM_BSHQI in the
context of the AIS dataset.

These results generally highlight the sensitivity of the mixture model
algorithm to the choice of initialization method. While the algorithm
CopMixM_BSHQI generally works well, it is critical to consider the impact of
initialization on its effectiveness. The higher miss-classification rate
observed with K-Means initialization should be taken into account when
implementing CopMixM_BSHQI in scenarios similar to the AIS dataset and its
further investigation will be object of future work especially as it seems to
improve by setting a specific number of bins. Indeed, by choosing the Rice
rule, the CopMixM_BSHQI results to be very competitive with respect to both
initializations.

####  5.2.2 Breast Cancer Wisconsin (Diagnostic)

The Breast Cancer Wisconsin (Diagnostic) dataset from the UCI Machine Learning
Repository [47] consists of digitized images of fine needle aspirates from
breast masses from 569 patients. For each of the considered ten features, the
mean value, extreme value (mean of the three largest values), and standard
error are computed, resulting in 30 total features. The dataset comprises of
benign (352 patients) and malignant (212 patients) diagnoses, enabling the
measurement of miss-classification rates for binary classification algorithms.
In this study, the same features as in [51] were considered, specifically:
perimeter standard error (PSE), extreme smoothness (ES), extreme concavity
(EC), and extreme concave points (ECP).

The copMxM_BSHQI algorithm was applied with both random initialization and
K-Means and by setting two different numbers for the bins. Our results are
reported in the Table 13, and a comparison is made with the outcomes presented
in the study in [51].

K-Means VCMM k-means VCMM (C-vine) Multivariate normal Multivariate skew
normal Multivariate t Multivariate skew t CopMixM_BSHQI k-means CopMixM_BSHQI
random Miss-classification rate Bins=$\lceil n^{1/3}\rceil$ 0.14 0.10 0.18
0.12 0.15 0.11 0.15 0.09 0.10 Miss-classification rate (Bins=Rice) 0.14 0.10
0.18 0.12 0.15 0.11 0.15 0.08 0.08

Table 13: Results for Breast Cancer Dataset. In bold the best values.

The results indicate that CopMixM_BSHQI with K-Means and with the two choices
of bins outperforms the other approaches, achieving a low miss-classification
rate of 0.09 and 0.084, respectively. In this case, with the random
initialization we can observe a miss-classification rate $0.082$ only by
choosing the Rice bins. This indicate a high degree of accuracy in grouping
data points, demonstrating the capability of our methodologies to discern
underlying patterns within the dataset. In contrast, the other methods from
[51], such as VCMM, VCMM (C-vine), Multivariate normal, Multivariate skew
normal, Multivariate t, and Multivariate skew t, exhibit comparatively higher
miss-classification rates ranging from 0.10 to 0.18. The results show the
effectiveness of the CopMixM_BSHQI approach, when combined with either K-Means
or random initialization.

####  5.2.3 Text Clustering

For the last experiment, we focus on the well-known 20newsgroups dataset,
accessible through Scikit-Learn’s API and designed for linguistic analysis.
This dataset contains online discussions, or newsgroups, categorized into
different thematic groups. For our clustering analysis, we specifically
selected texts related to technology, religion, and sports.

To convert the textual data into a format suitable for quantitative analysis,
we employ the TfidfVectorizer (implemented in scikit-learn library in
Python666https://scikit-
learn.org/stable/modules/generated/sklearn.feature_extraction.text.TfidfVectorizer.html).
This essential tool facilitates the transformation of text into numerical
vectors by calculating Term Frequency-Inverse Document Frequency (TF-IDF)
values for each term. TF-IDF reflects the importance of each term within
individual documents relative to the entire dataset, allowing us to capture
semantic information.

It is crucial to note that the TfidfVectorizer is also pivotal for the
clustering process. We transform textual data into numerical representations,
creating a dataset with $18,846$ rows and $24,471$ features. Considering the
big size of the dataset, a preliminary step in our clustering analysis
involves the application of a dimensionality reduction technique.
Specifically, we employ Truncated Singular Value Decomposition
(T-SVD)777https://scikit-
learn.org/stable/modules/generated/sklearn.decomposition.TruncatedSVD.html
with a predefined number of components set to 2. This preliminary step is
fundamental to mitigate the challenges posed by high-dimensional data,
allowing for a more efficient and manageable representation of the underlying
structure.

We present scatter plots for GMM, CopMixM_BSHQI and CopMixM_KDEpy in Figure 9.
Additionally, the Table 14 provides a comparative overview of clustering
performance between K-Means, GMM, CopMixM_BSHQI and CopMixM_KDEpy algorithms
applied to the 20newsgroup dataset. Again the proposed procedure achieves the
best results with respect to Silhouette score and Davies-Bouldin Score, while
is the runner up for the Calinski-Harabasz score, but we observe a small
discrepancy with the best result provided by K-Means.

(a) GMM

(b) CopMixM_BSHQI

(c) CopMixM_KDEpy

Figure 9: Text Clustering results.

K-Means GMM CopMixM_BSHQI CopMixM_KDEpy Silhouette Score 0.423 0.154 0.425
0.373 Calinski-Harabasz Score 2326 386 2200 1900 Davies-Bouldin Score 0.896
1.89 0.849 0.856

Table 14: Clustering metrics for 20newsgroup Dataset. In bold the best values.

The highest Silhouette Score emphasizes a better ability to create well-
defined clusters similarly, the highest Davis-Boudin score assesses a better
compactness and separation between clusters. The Calinski-Harabasz Score
reflects the effectiveness in achieving a favorable ratio of between-cluster
to within-cluster variance. In summary, across multiple clustering metrics,
CopMixM_BSHQI outperforms GMM, with random initialization, as highlighted by
the values in bold and provides a correct interpretation of the undrlying
statistical distribution.

##  6 Conclusions

In this paper we presented a novel algorithm for empirical density estimation
and we used it for cluster modeling based on the use of Copulas. In
particular, the multivariate copulas distribution rely on the estimation of
the marginal distributions based on the Hermite Quasi-interpolant in [30]. The
proposed construction is superior in terms of statistical significance with
respect to classical approaches based on empirical kernel density estimation
and provides consistent cumulative distribution functions as outlined in the
detailed analysis of Section 2. The novel clustering algorithm allows for the
automatic selection of the best copula among a certain set of copulas families
and provides a robust strategy with respect to a random seed as
initialization. Moreover, the obtained results show a rather good agreement
with the ground-truth (when provided), and mostly we are able to correctly
identify the underlying statistical distribution from where the given points
where drawn. The obtained clusters exhibit good shape parameters, in terms of
Silhouette Score, Calinski-Harabasz Score and Davis-Bouldin Score, and achieve
good accuracy in terms of permutation invariant metrics. Future work will be
devoted to investigate the choice of the optimal bandwidth and further to
deeply analyze how to deal with overlapping clusters, like in the text mining
example.

## Acknowledgments

Cristiano Tamborrino and Francesca Mazzia acknowledge the support of the PNRR
project FAIR - Future AI Research (PE00000013), Spoke 6 - Symbiotic AI (CUP
H97G22000210007) under the NRRP MUR program funded by the NextGenerationEU.
The research of Antonella Falini is founded by PON Ricerca e Innovazione
2014-202 FSE REACT-EU, Azione IV.4 “Dottorati e contratti di ricerca su
tematiche dell’innovazione” CUP H95F21001230006. The authors thank the GNCS
for its valuable support under the INdAM-GNCS project CUP E55F22000270001.

## References

  * [1] L. M. Abualigah, A. T. Khader, and E. S. Hanandeh, A new feature selection method to improve the document clustering using particle swarm optimization algorithm, Journal of Computational Science, 25 (2018), pp. 456–466.
  * [2] T. W. Anderson, On the Distribution of the Two-Sample Cramer-von Mises Criterion, The Annals of Mathematical Statistics, 33 (1962), pp. 1148 – 1159.
  * [3] G. Andrew and J. Gao, Scalable training of regularized log-linear models, in Proceedings of the 24th International Conference on Machine Learning, New York, NY, USA, 2007, Association for Computing Machinery, p. 33–40.
  * [4] A. AZZALINI, A note on the estimation of a distribution function and quantiles by a kernel method, Biometrika, 68 (1981).
  * [5] C. M. Bishop, Pattern Recognition and Machine Learning (Information Science and Statistics), Springer-Verlag, Berlin, Heidelberg, 2006.
  * [6] G. Celeux and G. Govaert, A classification em algorithm for clustering and two stochastic versions, Computational Statistics & Data Analysis, 14 (1992), pp. 315–332.
  * [7] C. Dagnino, P. Lamberti, and S. Remogna, On spline quasi-interpolation through dimensions, ANNALI DELL’UNIVERSITA’DI FERRARA, 68 (2022), pp. 397–415.
  * [8] C. De Boor and G. Fix, Spline approximation by quasiinterpolants, Journal of Approximation Theory, 8 (1973), pp. 19–45.
  * [9] 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 (1977), pp. 1–38.
  * [10] F. M. L. Di Lascio, F. Durante, and R. Pappadà, Copula–based clustering methods, in Copulas and Dependence Models with Applications, M. Úbeda Flores, E. de Amo Artero, F. Durante, and J. Fernández Sánchez, eds., Cham, 2017, Springer International Publishing, pp. 49–67.
  * [11] F. Durante and C. Sempi, Principles of copula theory, CRC press, 2015.
  * [12] A. Ezugwu, A. Ikotun, O. Olaide, L. Abualigah, O. Agushaka, C. Eke, and A. Akinyelu, A comprehensive survey of clustering algorithms: State-of-the-art machine learning applications, taxonomy, challenges, and future research prospects, Engineering Applications of Artificial Intelligence, 110 (2022), p. 104743.
  * [13] A. Fahad, N. Alshatri, Z. Tari, A. Alamri, I. Khalil, A. Y. Zomaya, S. Foufou, and A. Bouras, A survey of clustering algorithms for big data: Taxonomy and empirical analysis, IEEE Transactions on Emerging Topics in Computing, 2 (2014), pp. 267–279.
  * [14] K. R. Gehringer and R. A. Redner, Nonparametric probability density estimation using normalized b–splines, Communications in Statistics - Simulation and Computation, 21 (1992), pp. 849–878.
  * [15] C. GENEST, K. Ghoudi, and L.-P. Rivest, A semiparametric estimation procedure of dependence parameters in multivariate families of distributions, Biometrika, 82 (1995).
  * [16] D. C. Gilles. Celeux and J. Diebolt, Stochastic versions of the em algorithm: an experimental study in the mixture case, Journal of Statistical Computation and Simulation, 55 (1996), pp. 287–314.
  * [17] A. K. Jain, Data clustering: 50 years beyond k-means, Pattern Recognition Letters, 31 (2010), pp. 651–666.  Award winning papers from the 19th International Conference on Pattern Recognition (ICPR).
  * [18] H. Joe, Dependence modeling with copulas, CRC press, 2014.
  * [19] H. Joe and J. J. Xu, The estimation method of inference functions for margins for multivariate models, 10 1996.
  * [20] S. R. Kasa, S. Bhattacharya, and V. Rajan, Gaussian mixture copulas for high-dimensional clustering and dependency-based subtyping, Bioinformatics, 36 (2019), pp. 621–628.
  * [21] J. Kirkby, A. Leitao Rodriguez, and D. Nguyen, Spline local basis methods for nonparametric density estimation, Statistics Surveys, 17 (2023).
  * [22] I. Kosmidis and D. Karlis, Model-based clustering using copulas with applications, Statistics and Computing, 26 (2014).
  * [23] H. Läuter, Silverman, bw: Density estimation for statistics and data analysis. chapman & hall, london–new york 1986, 175 pp.,£ 12.—, Biometrical Journal, 30 (1988), pp. 876–877.
  * [24] B.-G. Lee, T. Lyche, and K. Mørken, Some examples of quasi-interpolants constructed from local spline projectors, Mathematical methods for curves and surfaces, Oslo, 2000 (2000), pp. 243–252.
  * [25] Q. Li and J. S. Racine, Nonparametric econometrics: Theory and practice, Economics & Finance, 23 (1994), pp. 2059–2078.
  * [26] T. Lyche and L. L. Schumaker, Local spline approximation methods, Journal of Approximation Theory, 15 (1975), pp. 294–325.
  * [27] R. Malouf, A comparison of algorithms for maximum entropy parameter estimation, in Proceedings of the 6th Conference on Natural Language Learning - Volume 20, COLING-02, USA, 2002, Association for Computational Linguistics, p. 1–7.
  * [28] I. Marschner, On stochastic versions of the em algorithm, Biometrika, 88 (2001).
  * [29] F. J. Massey, The kolmogorov-smirnov test for goodness of fit, Journal of the American Statistical Association, 46 (1951), pp. 68–78.
  * [30] F. Mazzia and A. Sestini, The bs class of hermite spline quasi-interpolants on nonuniform knot distributions, BIT Numerical Mathematics, 49 (2009), pp. 611–628.
  * [31] G. Mclachlan and T. Krishnan, The EM Algorithm and Extensions (Wiley Series in Probability and Statistics), 01 2007.
  * [32] G. J. McLachlan and D. Peel, Finite mixture models, vol. 299 of Probability and Statistics – Applied Probability and Statistics Section, Wiley, New York, 2000.
  * [33] X.-L. MENG and D. B. RUBIN, Maximum likelihood estimation via the ECM algorithm: A general framework, Biometrika, 80 (1993), pp. 267–278.
  * [34] E. Nadaraya, Some new estimates for distribution functions, Theory of Probability and Its Applications - THEOR PROBAB APPL-ENGL TR, 9 (1964).
  * [35] R. B. Nelsen, An Introduction to Copulas (Springer Series in Statistics), Springer-Verlag, Berlin, Heidelberg, 2006.
  * [36] V. Rajan and S. Bhattacharya, Dependency clustering of mixed data with gaussian mixture copulas, in Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, IJCAI’16, AAAI Press, 2016, p. 1967–1973.
  * [37] R. Redner, Convergence rates for uniform b-spline density estimators part i: One dimension, Siam Journal on Scientific Computing, 20 (1999).
  * [38] R. A. Redner and K. Gehringer, Function estimation using partitions of unity, Communications in Statistics - Theory and Methods, 23 (1994), pp. 2059–2078.
  * [39] R. A. Redner and H. F. Walker, Mixture densities, maximum likelihood and the em algorithm, SIAM Review, 26 (1984), pp. 195–239.
  * [40] M. Rosenblatt, Remarks on Some Nonparametric Estimates of a Density Function, The Annals of Mathematical Statistics, 27 (1956), pp. 832 – 837.
  * [41] P. Sablonniere, Recent progress on univariate and multivariate polynomial and spline quasi-interpolants, Trends and applications in constructive approximation, (2005), pp. 229–245.
  * [42] A. Saxena, M. Prasad, A. Gupta, N. Bharill, O. P. Patel, A. Tiwari, M. J. Er, W. Ding, and C.-T. Lin, A review of clustering techniques and developments, Neurocomputing, 267 (2017), pp. 664–681.
  * [43] M. A. Tanner and W. H. Wong, The calculation of posterior distributions by data augmentation, Journal of the American Statistical Association, 82 (1987), pp. 528–540.
  * [44] A. Tewari, M. Giering, and A. Raghunathan, Parametric characterization of multimodal distributions with non-gaussian modes, vol. 11, 12 2011, pp. 286–292.
  * [45] M. P. Wand and M. C. Jones, Kernel Smoothing, no. 60 in Chapman & Hall/CRC Monographs on Statistics & Applied Probability, Chapman & Hall, Boca Raton, FL, U.S., December 1994.
  * [46] G. C. G. Wei and M. 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 (1990), pp. 699–704.
  * [47] S. N. Wolberg William, Mangasarian Olvi and S. W., Breast Cancer Wisconsin (Diagnostic).  UCI Machine Learning Repository, 1995.
  * [48] R. Xu and D. Wunsch, Survey of clustering algorithms, IEEE Transactions on Neural Networks, 16 (2005), pp. 645–678.
  * [49] H. e. a. Yamato, Uniform convergence of an estimator of a distribution function, Bulletin of Mathematical Statistics, 15 (1973).
  * [50] Y. Zhou, H. Wu, Q. Luo, and M. Abdel-Baset, Automatic data clustering using nature-inspired symbiotic organism search algorithm, Knowledge-Based Systems, 163 (2019), pp. 546–557.
  * [51] Özge Sahin and C. Czado, Vine copula mixture models and clustering for non-gaussian data, Econometrics and Statistics, 22 (2022), pp. 136–158.  The 2nd Special issue on Mixture Models.