Nonparametric Multivariate Density Estimation: Case Study of Cauchy Mixture Model

Abstract

Estimation of probability density functions (pdf) is considered an essential part of statistical modelling. Heteroskedasticity and outliers are the problems that make data analysis harder. The Cauchy mixture model helps us to cover both of them. This paper studies five different significant types of non-parametric multivariate density estimation techniques algorithmically and empirically. At the same time, we do not make assumptions about the origin of data from any known parametric families of distribution. The method of the inversion formula is made when the cluster of noise is involved in the general mixture model. The effectiveness of the method is demonstrated through a simulation study. The relationship between the accuracy of evaluation and complicated multidimensional Cauchy mixture models (CMM) is analyzed using the Monte Carlo method. For larger dimensions (d ~ 5) and small samples (n ~ 50), the adaptive kernel method is recommended. If the sample is n ~ 100, it is recommended to use a modified inversion formula (MIDE). It is better for larger samples with overlapping distributions to use a semi-parametric kernel estimation and more isolated distribution-modified inversion methods. For the mean absolute percentage error, it is recommended to use a semi-parametric kernel estimation when the sample has overlapping distributions. In the smaller dimensions (d = 2) and a sample is with overlapping distributions, it is recommended to use the semi-parametric kernel method (SKDE) and for isolated distributions, it is recommended to use modified inversion formula (MIDE). The inversion formula algorithm shows that with noise cluster, the results of the inversion formula improved significantly.

1. Introduction

Estimation of probability density functions (pdf) is considered an essential part of statistical modelling. It expresses random variables as functions of other variables, making it possible to detect hidden relationships between data [1]. In a significant number of machine learning algorithms, it is essential to determine a previously unknown function of the distribution density of the data. The function of the distribution density is applied in the Bayesian classifier [2,3], in density-based clustering algorithms [4,5,6], or information-based feature selection algorithms [7,8]. Effective density estimates must be carefully created in advance to obtain unknown functions of probability density. Nowadays, there is still much focus on developing innovative density estimation procedures [9,10]. Density estimation is an open research topic in the fast-growing area of deep learning. Scientists have begun proposing robust density estimators based on neural networks such as Parzen neural networks [11], soft-constrained neural networks [12], and others [13].

Let us say that the random vector X∈R^d satisfies the distribution mixture model if its distribution density f(x) satisfies the equation $f\left(x\right)={{\displaystyle \sum }}_{k=1}^q{p}_k{f}_k\left(x\right)=f\left(x,\theta \right)$. The parameter q is the number of mixture clusters, and p_k is the a priori probability. These conditions must also be met: $p_k>0$ and ${{\displaystyle \sum }}_{k=1}^q{p}_k=1$. The function $f_k\left(x\right)$ is a function of the distribution density, and $\theta$ is a multidimensional parameter of the model. Suppose $X$ is a d-dimensional random vector with a distribution density $f\left(x\right)$, and there is a sample of independent copies of $X$, where $X=\left(X\left(1\right),\dots ,X\left(n\right)\right)$. It can be argued that the sample satisfies the mixture model if $X\left(t\right)$ satisfies $f\left(x\right)={{\displaystyle \sum }}_{k=1}^q{p}_k{f}_k\left(x\right)=f\left(x,\theta \right)$.

One of the statistical tasks is to estimate the density of the observed random variable. Suppose the available sample’s distribution type is known (Normal, Poisson, and others). In that case, the distribution density of the data can be estimated simply using mean and covariance matrix estimates, fitting them to a defined distribution [14,15,16]. Thus, the standard parametric method is applied when the assumptions about the density form are met. When estimating density in a parametric way, the value of the multidimensional distribution parameter θ needs to be found, which is not straightforward because the number of parameters increases rapidly as the dimension d increases. For example, in the case of a mixture of Gaussian distributions, $\dim \theta =\frac{1}{2qd\left(d+1\right)}+qd+q-1$, and even with a small dimension $d=q=5$, the model will consist of $\dim \left(\theta \right)=104$ parameters. When searching for parameter estimates, it may be necessary to solve the optimization problem in the 104-dimensional space. In practice, the number of clusters $q$ may also be unknown, and it needs to be estimated. The parametric method is not proper when the random size distribution is unknown. In this case, non-parametric methods are used to determine certain forms of density estimates [17,18,19].

The histogram is one of the simplest and oldest estimates of density. To the best of our knowledge, data in the form of histograms (without graphical representation) were first presented in 1661 to determine mortality probabilities in different age groups [20]. To approximate the density $f\left(x\right)$ in the area $\mathsf{\Omega }$, the number of observations $X\left(t\right)$ falling into $\mathsf{\Omega }$ is calculated and divided by n and the volume of the area $\mathsf{\Omega }$. The area of space to which all observations fall is first found. That means the fluctuation intervals of all $X$ projections on the axes $X^{\left(1\right)}$,$X^{\left(2\right)}$,…,$X^{\left(d\right)}$ are found. The fluctuation intervals of the observations are divided into l partial intervals and in the hypercubes ${\mathsf{\Omega }}_j\left(j=1,\dots ,r\right)$ bounded by them, the density estimate is calculated as

$$\widehat{f}\left(x\right)=\frac{n\left({\mathsf{\Omega }}_j\right)}{n\cdot h_1\cdot h_2\dots \cdot h_d}.$$

(1)

Here $n\left({\mathsf{\Omega }}_j\right)$ is the number of observations entering the hypercube ${\mathsf{\Omega }}_j$ and $h_j,j=1,\dots ,d$ are the edges of the hypercube [21,22]. It is recommended to select the number of hypercubes [17,23,24], and to choose $r\cong 1+3.32\log \left(n\right)$, and $l=\sqrt[d]{r}$ has to be an integer number, so $r$ is chosen that ${\lceil \sqrt[d]{1+3.32logn}}^d\rceil$.

A histogram is one of the simplest means of presenting data that is easy to understand and convenient. This estimate is a function that acquires non-negative values, and its integral is equal to one. However, it is not continuous. That poses problems when knowing the density estimate derivatives is essential, mainly when density estimation is used in intermediate steps of other methods, such as clustering using a gradient algorithm or plotting high-measurement data-level lines. Remarkably, the histogram stood as the only non-parametric density estimator until the 1950’s when substantial and simultaneous progress was made in density estimation and spectral density estimation. In 1951, in a little-known paper, Fix and Hodges [25] introduced the basic algorithm of non-parametric density estimation; an unpublished technical report was formally published as a review by Silverman and Jones in 1989 [26]. They addressed the problem of statistical discrimination when the parametric form of the sampling density was not known. During the following decade, several general algorithms and alternative theoretical modes of analysis were introduced by Rosenblatt in 1956 [27], Parzen in 1962 [28], and Cencov in 1962 [29]. Then followed the second wave of essential and primarily theoretical papers by Watson and Leadbetter in 1963 [30], Loftsgaarden and Quesenberry in 1965 [31], Schwartz in 1967 [32], Epanechnikov in 1969 [33], Tarter and Kronmal in 1970 [34], and Kimeldorf and Wahba in 1971 [35]. Next, Cacoullos introduced the natural multivariate generalization in 1966 [36]. Finally, in the 1970s, the first papers focusing on the practical application of these methods were published by Scott et al. in 1979 [24] and Silverman in 1978 [37]. These and later multivariate applications awaited the computing revolution.

Modern data analysis uses several non-parametric methods for statistically estimating the distribution density of multivariate random variables. Kernel estimates are particularly prevalent [38,39]. Quite popular and spline [40,41] and semi-parametric [42,43,44,45,46] algorithms. However, detailed comparisons of the effectiveness of existing popular estimates for multimodal density are lacking. With the most popular non-parametric estimation procedures, optimal selection of their parameters is encountered in practice. The most crucial element in the design of kernel estimates is the width of the smoothing. It is not easy to select the nodes of the spline estimates. Although several adaptive procedures for the selection of these parameters have been developed [39,47,48,49,50,51,52], however, they are not efficient enough when the sample volume is not large, especially then the observational dimension is large. In the latter case, it is appropriate to apply data design [53,54,55,56] because of the more extensive the dimension of the observed random vectors, the more complex the task of parameter selection.

The main idea of this paper is to estimate the performance of different density estimators by using density mixtures to show another type of problem, which may result from data heteroscedasticity and outliers. The relationship between the accuracy of evaluation and complicated multidimensional Cauchy mixture models (CMM) is analyzed using the Monte Carlo method. For example, Kalantan and Einbeck [57] used engineering data and, for computer vision, used CMM, comparing it with the Gaussian mixture model. Azzari and Foi [58] used harmony between Gauss and heavy-tailed Cauchy to find noise-model parameters that make outlier estimation robust when imaged dominated by texture. Finally, Teimouri [59] analyzed patients with Cushing’s syndrome and their diagnostic tests. The focus was on the tetra hydrocortisone urine release rate (mg/24 h) and evaluating parameters in the EM algorithm and Cauchy mixture model.

Scientific novelty. Evaluation accuracy comparative analysis is made by using different probability density estimation procedures. Density function estimates are chosen as popular different technique estimates, which other researchers have already analyzed. This research is essential because it focuses on Cauchy distributions.

2. The Density Estimation Algorithms

This section aims to present the density estimation algorithms used in the study theoretically. All algorithms are presented with algorithms theoretical substantiation. When making the histogram, each $X\left(t\right)$ can be imagined as a separate column with a height of $1/n$. Then it makes sense to change the centre of the column to $X\left(t\right)$ itself and get the following function:

$$\widehat{f}\left(x\right)=\frac{1}{n\cdot h_1\cdot h_2\dots \cdot h_d}{{\displaystyle \sum }}_{t=1}^n{I}_{C_h\left(X\left(t\right)\right)}\left(x\right).$$

(2)

Here $C_h$ is a hypercube with centre $X\left(t\right)$, and the lengths of the edges are $h_1,\dots ,h_d$. In summary, instead of the indicator function, a smooth “prominence”—the kernel function—can be used at each observed point. The multidimensional fixed-width bandwidth estimate with the kernel function K and the fixed (global) kernel width parameter h, which can be used to estimate the density $\widehat{f}\left(x\right)$ of the multidimensional data X∈R^d, is then defined as follows:

$$\widehat{f}\left(x\right)=\frac{1}{n{h}^d}{{\displaystyle \sum }}_{t=1}^{n}K\left(\frac{x-X\left(t\right)}{h}\right).$$

(3)

These are some of the most common non-parametric estimates of distribution density [38,39,60,61]. The kernel function is selected to meet the following condition:

$${{\displaystyle \int }}_{R^d}K\left(x\right)dx=1.$$

(4)

The standard normal distribution density function $\phi$ is often used as the kernel [62,63]:

$$\varphi \left(x\right)=(2\pi {)}^{-d/2}\exp (-\frac{1}{2}{x}^{\prime }x).$$

(5)

Often the observations are not evenly distributed in all directions. Therefore, it is desirable to scale the data by eliminating the most significant dispersion differences in the different coordinate directions. One suitable method for this [64] is data standardization. That means the sample’s effect on a linear transformation. The mean of the transformed data is zero, the covariance matrix is unitary, and (3) apply the Equation to already standardized data. For example, suppose Z is a standardized random vector,

$$Z=S^{-\frac{1}{2}}\left(X-\overline{X}\right),$$

(6)

here $\overline{X}$ is the empirical mean of the sample, and S∈R^d×d is the empirical covariance matrix. Based on the fixed kernel width density estimate (3), a more complex standardized data density estimate has been constructed:

$$\widehat{f_z}\left(z\right)=\frac{1}{n{h}^d}{{\displaystyle \sum }}_{t=1}^{n}K\left(\frac{z-Z\left(t\right)}{h}\right).$$

(7)

$$\acute{f}\left(x\right)=\frac{{\left(detS\right)}^{-1/2}}{n{h}^d}{{\displaystyle \sum }}_{t=1}^{n}K\left(S^{-1/2}\frac{x-X\left(t\right)}{h}\right).$$

(8)

The optimal kernel width $h^{*}$ for a fixed core width is determined by minimizing the average integral root mean square error (MISE) [65]. For example, when the distribution of observations is normal with a unit covariance matrix in Gaussian kernel, the expression $h^{*}$ proposed by [65] is $h^{*}=A{n}^{-\frac{1}{d+4}}$, here $A={[4/\left(2d+1\right)]}^{\frac{1}{d+4}}$. More sophisticated kernel width selection methods (such as the least-squares cross-checking method) are obtained by more complex and lengthy calculations [66,67,68,69,70].

In practical research, the kernel width is often selected experimentally. If the value of h is small, the density function estimate has more modes that correspond to the layout of the observed data. A higher value of h means more significant smoothing of the estimate.

Although fixed-core width density estimates are widely used to estimate non-parametric densities, they often have some practical drawbacks [65]. For example, fixed-core width density estimates do not ensure the distribution ends’ integrity without over-smoothing the underlying bulk density.

2.1. Adapted Kernel Density Estimate (AKDE)

A good improvement on the fixed kernel width density estimate is the adapted kernel density estimate [65]. The adapted kernel density estimate is constructed similarly to the fixed kernel width density estimate. The kernel describes the density at each observed point. In this case, the kernel width is already considered when moving from one observation to another. In areas of different smoothness, it is appropriate to take different kernel widths. This method consists of two steps: estimation of the adapted kernel width and density estimation by the kernel method, using the information obtained in the first step. The algorithm can be summarized as follows:

Step 1: The elements of sample $X=\left(X\left(1\right),\dots ,X\left(n\right)\right)$ are standardized to $Z=\left(Z\left(1\right),\dots ,Z\left(n\right)\right)$ such that $\widehat{E}\left[Z\right]=0$ and $\widehat{E}\left[Z{Z}^{\prime }\right]=I$.

Step 2: Estimates ${\tilde{f}}_Z\left(z\right)$ of the fixed kernel density estimate (3) satisfying the condition ${\tilde{f}}_Z\left(Z\left(t\right)\right)>0$, ∀t.

Step 3: The local width parameter is determined ${\lambda }_t={\left(\frac{{\tilde{f}}_Z\left(Z\left(t\right)\right)}{g}\right)}^{-\gamma }$, where g is ${\tilde{f}}_Z\left(z\right)$ the geometric mean, $\log g=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\log \widehat{f_Z}\left(Z\left(t\right)\right)$ and γ is the sensitivity parameter: $0\le \gamma \le 1$.

Step 4: An adapted kernel estimate is made with variable-width kernels: ${\widehat{f}}_Z\left(z\right)=\frac{1}{n}{{\displaystyle \sum }}_{t=1}^n{h}^{-d}{\lambda }_t^{-d}K\left(\frac{z-Z\left(t\right)}{h{\lambda }_t}\right)$.

Where h is the same global smoothness parameter as in Equation (3), the higher $\gamma$, the more sensitive the density selection. Quite often, the parameter value is selected as follows $\gamma =\frac{1}{2}$ [65,71].

2.2. Semi-Parametric Kernel Density Estimate (SKDE)

When data are scarce, parametric estimates are often applied even when the unknown density is not parameterized. Therefore, it is essential to mention the combination of parametric and non-parametric estimates. For example, one of the semi-parametric estimates of kernel density was examined by F. Hoti and L. Holmström [46]. This estimate divides the random vector into two subvectors and estimates the distribution density of one of them by the kernel method. Afterward, another relative density is approximated by the Normal distribution density [46]. For example, suppose d and s are positive integers d ≥ 2, 1 ≤ s ≤ d − 1. Using this method, the d-dimensional vector X∈R^d is decomposed into two s and (d–s) dimensioned subvectors $X=\left(\begin{array}{c}Y\\ Z\end{array}\right)$, and the sample is decomposed accordingly: $X=\left(\begin{array}{c}Y\\ Z\end{array}\right)$, where Y∈R^s, Z∈R^d-s. The evaluated density function is expressed as the product of the distribution density of the random vector Y and the conditional distribution density of the random vector Z: $f_X\left(x\right)=f_{\left(Y,Z\right)}\left(y,z\right)=f_Y\left(y\right)f_{Z|Y=y}(z|y)$, $x=\left(\begin{array}{c}y\\ z\end{array}\right)\in R^d$. Here $f_X$ and $f_Y$ are the densities of $X$ and $Y$. $f_{Z|Y=y}$ is the density of $Z$ when $Y=y$.

Suppose that the relative density $Y=y$ is multidimensional normal Gaussian, but the density $f_Y$ does not belong to any family of parametric functions. The density $f_X$ is then obtained by estimating $f_Y$ in a non-parametric manner and applying a multidimensional Normal density to each $f_{Z|Y=y}$. The density function $f_Y\left(y\right)$, as with (8), is evaluated by the kernel method [65]. Since the sample elements are not standardized, the smoothness parameter is not the same in all directions. Therefore, using the kernel method, it is replaced by the s-dimensional matrix H:

$$\widehat{f}\left(y\right)=\frac{1}{n}{{\displaystyle \sum }}_{t=1}^n\frac{1}{\det (H)}K\left(H^{-1}\left(y-Y\left(t\right)\right)\right).$$

(9)

Usually, the shape of H is chosen diagonally—$H=diag\left(h_1,\dots ,h_s\right)$, and the smoothness parameters are calculated as follow

$$h_j={\left(\frac{4}{s+2}\right)}^{1/\left(s+4\right)}{n}^{-1/\left(s+4\right)}{\sigma }_j.$$

(10)

It should be noted that this form, when s = 1, was proposed by B. W. Silverman [65]. Replacing the standard deviation σ_j of the component Yj with its estimate ${\widehat{\sigma }}_j=\sqrt{\frac{\sum {\left(X_j-\overline{X_j}\right)}^2}{n_j}}$ and by the rule of D. W. Scott [39] first multiplier is always between 0.924 and 1.059 ${\widehat{h}}_j$ can be calculated as follows

$${\widehat{h}}_j=n^{-1/\left(s+4\right)}{\widehat{\sigma }}_j.$$

(11)

This Scott’s rule is easy to summarize for the smoothness matrix H:

$$\widehat{H}=n^{-1/\left(s+4\right)}{ \mathsf{\Sigma } ^}^{1/2}.$$

(12)

Here $\widehat{\Sigma }=\mathrm{diag}\left({\widehat{\sigma }}_1^2,{\widehat{\sigma }}_2^2,\dots ,{\widehat{\sigma }}_s^2\right)$ is the diagonal matrix of Y empirical variances.

The conditional density $f_{Z|Y}(\cdot |y)$ is approximated by the Gaussian distribution N(m(y), C(y)), where m(y), C(y) denote the conditional mean of the vector Y and the conditional covariance matrix: $m\left(y\right)=E(Z|Y=y)$, y∈R^s, $C\left(y\right)=E\left[{(Z-m\left(y\right))}^{\prime }\left(Z-m\left(y\right)\right)|Y=y\right]$, y∈R^s. For the estimation of m(y) and C(y), it is proposed to apply the kernel smoothing:

$$\widehat{m}\left(y\right)=\frac{{{\displaystyle \sum }}_{t=1}^n{K}_{H_2}\left(y-Y\left(t\right)\right)Z\left(t\right)}{{{\displaystyle \sum }}_{j=1}^n{K}_{H_2}\left(y-Y\left(j\right)\right)}={{\displaystyle \sum }}_{t=1}^n{W}_{H_2}\left(y-Y\left(t\right)\right)Z\left(t\right), y\in R^s.$$

(13)

Here are the weights $W_{H_2}\left(y-Y\left(t\right)\right)=\frac{K_{H_2}\left(y-Y\left(t\right)\right)}{{{\displaystyle \sum }}_{j=1}^n{K}_{H_2}\left(y-Y\left(j\right)\right)}$.

The sum of which is equal to one. The formula (13) can be understood as a regression estimate of the conditional mean function m of Nadaraya and Watson [72,73]. The conditional covariance matrix can be evaluated similarly $\widehat{C}\left(y\right)={{\displaystyle \sum }}_{t=1}^n{W}_{H_3}\left(y-Y\left(t\right)\right){(Z\left(t\right)-\widehat{m}\left(y\right))}^{\prime }\left(Z\left(t\right)-\widehat{m}\left(y\right)\right)$, y∈R^s. The parametric estimate of the relative density f(Z|Y) = y looks akin to this ${\widehat{f}}_{Z|Y=y}\left(z\right)={[{(2\pi )}^{d-s}\widehat{\mathit{det}C}\left(y\right)]}^{-1/2}\mathit{exp}\left\{-\frac{1}{2}\left(z-\widehat{m}\left(y\right)\right)\widehat{C}{(y)}^{-1}{(z-\widehat{m}\left(y\right))}^{\prime }\right\}$, z∈R^d-s. The estimate of the distribution density f_X of X then is: ${\widehat{f}}_X\left(x\right)={\widehat{f}}_{\left(Y,Z\right)}\left(y,z\right)={\widehat{f}}_Y\left(y\right){\widehat{f}}_{Z|Y=y}\left(z\right)$, x = (y,z)∈R^d.

The procedure described above is called the semi-parametric kernel density estimate. In practice, even if the conditional assumption of the normality of several random vector components is satisfied. The decomposition dimensions also influence the accuracy of the density estimation results, and the choice of the coordinates influences the accuracy of the density estimation results. One way to select them is to use the least-squares method or the maximum likelihood cross-entropy method recommended by original method authors [46]. The authors propose the parameters H₂ and H₃ [46] to select 2H.

2.3. Logspline Estimation (LSDE)

This subsection describes the logspline estimation (LSDE) calculation. One-dimensional polynomial splines are called partial polynomials of a certain degree. Breakpoints that contain a transition from one polynomial to another are called nodes. Suppose that the vector $t=\left(t_1,\dots ,t_K\right)\in R^K$ defines a set of such K points. Splines describe smooth connections, showing how different areas are separated by nodes [74]. These constraints are precisely defined by expressing partial polynomials in the number of continuous derivatives s. These include partially linear curves. If there are no restrictions, breakpoints are allowed in the nodes of these functions. Assuming that the functions are globally continuous, it is required that the individual linear parts meet at each node. If greater smoothness is needed (for continuous first-order derivatives), then the flexibility of the nodes is lost. Moreover, the curves are considered simple linear functions. The term “linear spline” is applied to a continuous partial linear function in the literature on approximation theory. Accordingly, the term “cubic spline” is assigned to continuous cubic functions with second-order continuous derivatives and nodes that allow jumps of third-order derivatives. If the polynomial degree is b and the vector of the nodes is t, then the set of polynomial splines with s continuous derivatives forms a linear space. For example, a set of linear splines with nodes in the sequence t is defined by function

$$1,x,\left(x-t_1\right)+,\dots ,\left(x-t_K\right)+.$$

(14)

Here (٠)₊ = max(٠,0). We will rely on this set as the base of space. The base of the spline space of degree b and s smoothness consists of monomial whose form ${(x-t_k)}_{+}^{s+j}$, here 1 ≤ j ≤ b − s. Using this formula, in the case of classical cubic splines, where b = 3 and s = 2, the base consists of elements

$$1,x,x^2,x^3,{\left(x-t_1\right)}^3+,\dots ,{\left(x-t_K\right)}^3+.$$

(15)

From the model point of view, this base is convenient because the individual functions at the nodes are merged. In expressions (14) and (15), each function is precisely associated with one of the nodes, and removing this function essentially corresponds to removing the node itself. It is known that the numerical properties of functions (14) and (15) are poor. For example, the solution matrix deteriorates as rapidly as the number of nodes decreases in linear regression problems. A practical alternative is the so-called B-spline base [75,76]. These functions are designed to be supported in several contiguous intervals defined by nodes (b + 1 contiguous intervals are used for the smoothest splines). Suppose we can find the basis for splines of space B₁(x; t), …, B_J(x; t) b with smoothness s and a sequence of nodes t so that any function in space can be written as $g\left(x;\beta ,t\right)={\beta }_1{B}_1\left(x;t\right)+\dots +{\beta }_J{B}_J\left(x;t\right)$. Where the corresponding coefficient vector is β = (β₁, …, β_J)′. As seen from (14) and (15), then spline spaces of maximum smoothness are used J = K + b + 1.

According to the title of the subsection, the object of this analysis is the logarithmic density. Suppose X is a random vector that takes values from the interval (L, U). In the individual case, L and U can be ±∞. The parameters L and U are set to 2t₁ − t₂ and 2t_K − t_{K − 1}, respectively. If β₁ ≥ 0 or β_{K − 1}≥0, then the adjustment is made 2L_old − t₁ and U_new = 2U_old − t_K is performed. The method of Kooperberg and Stone [52,77,78,79], known as logspline, is implemented with cubic spline. The cubic spline is described in (15). These functions are also continuously differentiated, and the partial polynomials are defined accordingly in the sequence of nodes t = (t₁, …, t_K). In each interval [t₁, t₂], …, [t_K−1, t_K] cubic splines are also cubic polynomials, but at the edges (L, t₁] and [t_K, U) are linear functions. The minimum number of nodes is K ≥ 3 (otherwise, a linear function or constant can be obtained). The basis form is 1, B₁(x; t), …, B_J(x; t), where J = K − 1.

It is said that the vector β = (β₁, …, β_J)′∈R^J exists then $C\left(\beta ,t\right)=log\left({{\displaystyle \int }}_L^{U}exp({\beta }_1{B}_1\left(x;t\right)+\dots +{\beta }_J{B}_J\left(x;t\right))dx\right)<\infty$. Suppose B denotes a set of such possible vectors. After selecting β∈B, the family of positive density functions in the interval (L, U) is defined the form of which is

$$g\left(x;\beta ,t\right)=exp({\beta }_1{B}_1\left(x;t\right)+\dots +{\beta }_J{B}_J\left(x;t\right)-C\left(\beta ,t\right)), L<x<U.$$

(16)

Now, having a random sample n of magnitude X(1), …, X(n) from the interval (L, U) with an unknown density function f, the logical probability function corresponding to the model of logsplines (16) is

$$l\left(\beta ,t\right)={{\displaystyle \sum }}_{i}log\left(g\left(X_i;\beta ,t\right)\right)={{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{\beta }_j{B}_j\left(X_i;t\right)-nC\left(\beta ,t\right), \beta \in B.$$

(17)

where estimation of maximum likelihood $\widehat{\beta }=argma{x}_{\beta \in B}l\left(\beta ,t\right)$ and an estimate of density $\widehat{f}=g\left(x;\widehat{\beta },t\right)$, L < x < U.

Let us say that during the stepwise determination procedure, the sequence of models is denoted by v, the vth model has J_v base functions. The Generalized Akaike Information Criterion (AIC) selects the best model [80]. Suppose that ${\widehat{l}}_v$ defines the estimate of the logic-likelihood function (17) for the vth model. The Equation defines the Akaike information criterion ${\mathrm{AIC}}_{a,v}\left(t\right)=-2{\widehat{l}}_v\left(t\right)+a{J}_v$ for which the model has a loss parameter a. From many models, the one whose value of v minimizes AIC_a,v. Stone [52] recommends the use of a = log n.

2.4. PPDE Algorithm. Estimation of the Projection Density of the Target

The projection pursuit density estimator (PPDE) proposed by Friedman is based on the target projection and consistent projection Gaussianization. The essence of J. H. Friedman and coauthors [54,55,81] in estimating the target projection density is to search for “interesting”, small-measurement data projections. The distribution structures, where the projections have distributions that are very different (in the sense of some projection index) from Gaussian. Huber [82] made a heuristic proposal to consider the Gaussian distribution as the least interesting. This proposal is based on the facts that:

• The multidimensional Gaussian distribution is entirely defined by its linear structure (mean and covariance matrices). Therefore, it is desired to feel a data structure independent of the correlation and linear transformations such as the scale parameter.
• All projections of a multidimensional Gaussian distribution are also Gaussian distributions. Thus, if the projection differs insignificantly from the Gaussian distribution, it indicates that distribution is also close to the Gaussian.
• For multidimensional data with a structure in multiple projection directions, many projections will have a distribution close to normal. This statement follows from the central limit theorem.
• In the case of constant variance, the Gaussian distribution is considered to be the least informative.

Friedman developed Huber’s idea and proposed an algorithm called exploratory target projection to estimate multidimensional non-parametric density. This procedure consists of five steps:

(1)

Data standardization simplifies layout, scalability, and correlation structures;

(2)

Projection index: the degrees of ‘interest’ in various directions are determined.

(3)

Optimization strategy: search for the direction in which the projection index is the largest.

(4)

Data transformation: the one-dimensional density is calculated in the chosen direction, and the data are multiplied.

(5)

Density formation: multidimensional density is formed from the calculated one-dimensional densities. Multidimensional density is a function of one-dimensional densities.

The following projection index construction has been proposed. It is known that all projections of a multidimensional Gaussian distribution are one-dimensional Gaussian distributions. If the distribution in one direction is not Gaussian, then the multidimensional distribution is also not Gaussian. Therefore, the projection index $I\left(\tau \right)$ shows how far the one-dimensional density $f_{\tau }\left(y\right)$ is in the direction $\tau \left(Y={\tau }^{\prime }Z\right)$ from the Gaussian distribution when $Z$ is a standardized quantity [83]:

$$\tilde{I}\left(\tau \right)={{\displaystyle \int }}_{-\infty }^{\infty }{(f_{\tau }\left(y\right)-\varphi \left(y\right))}^{2}dy, \mathrm{where} \varphi \left(y\right)=\frac{1}{\sqrt{2\pi }}{e}^{\frac{-y^2}{2}}.$$

(18)

The projection direction $\mathsf{\tau }$, which maximizes the projection of a distribution $\tilde{I}\left(\tau \right)$, is chosen to highlight the multimodal or other nonlinear structure of that distribution. We transform the data $y$ by equality $R=2\Phi \left(Y\right)-1=2\Phi \left({\tau }^{\prime }Z\right)-1, R\in \left[-1, 1\right]$, where $\Phi \left(u\right)$ is a function of the standard normal distribution $\Phi \left(u\right)=\frac{1}{\sqrt{2\pi }}{{\displaystyle \int }}_{-\infty }^u{e}^{\frac{-t^2}{2}}dt$. The distribution density of the transformed quantity $R$, function $f_R\left(r\right)$ can be rewritten as

$$f_R\left(r\right)=\frac{f_{\tau }\left(y\right)}{\left|\frac{\partial r}{\partial y}\right|}=\frac{f_{\tau }\left(y\right)}{2\varphi \left(y\right)}.$$

(19)

Equation (18) can be rewritten by changing the variable $y$ to $r: \tilde{I}\left(\tau \right)={{\displaystyle \int }}_{-1}^{1}2\varphi \left(y\right){(f_R\left(r\right)-1/2)}^{2}dr={{\displaystyle \int }}_{-1}^{1}2\varphi \left({\Phi }^{-1}\left(\frac{R+1}{2}\right)\right){(f_R\left(r\right)-1/2)}^{2}dr.$ Friedman [55] proposed a slightly different form of the projection index $I\left(\tau \right)$, taking the integrated square error as a measure of $R$ inequality:

$$I\left(\tau \right)={{\displaystyle \int }}_{-1}^1{(f_R\left(r\right)-1/2)}^{2}dr={{\displaystyle \int }}_{-1}^1{f}_R^2\left(r\right)dr-1/2.$$

(20)

Note that if the distribution of $Y$ is Gaussian, then $f_R\left(r\right)=\frac{1}{2}$, ∀r, and the projection index $I\left(\tau \right)$ is zero. The more the $Y$ distribution differs from the normal, the higher the value of the index $I\left(\tau \right)$. Since $R\in \left[-1, 1\right]$, $f_R\left(r\right)$ can be decomposed by orthogonal Lagrangian polynomials, ${\left\{{\psi }_j\right\}}_{j=0}^{\infty }$, i.e., $f_R\left(r\right)={{\displaystyle \sum }}_{j=0}^{\infty }{b}_j{\psi }_j\left(r\right)$:

$$I\left(\tau \right)={{\displaystyle \int }}_{-1}^1{f}_R^2\left(r\right)dr-1/2={{\displaystyle \int }}_{-1}^1\left[{{\displaystyle \sum }}_{j=0}^{\infty }{b}_j{\psi }_j\left(r\right)\right]f_R\left(r\right)dr-1/2.$$

(21)

An iterative expression defines orthogonal Lagrangian polynomials. ${\psi }_0\left(r\right)=1$ and ${\psi }_1\left(r\right)=r. {\psi }_j\left(r\right)=\frac{\left(2j-1\right)r{\psi }_{j-1}\left(r\right)-\left(j-1\right){\psi }_{j-2}\left(r\right)}{j}$, then $j\ge 2$. It follows from the orthogonality property that the coefficients $b_j$ can be calculated as follows $b_j=\frac{2j+1}{2}{{\displaystyle \int }}_{-1}^1{\psi }_j\left(r\right)f_R\left(r\right)dr=\frac{2j+1}{2}{E}_R\left[{\psi }_j\left(r\right)\right]=\frac{2j+1}{2}\frac{1}{n}{{\displaystyle \sum }}_{t=1}^n{\psi }_j\left(2\Phi \left(Y\left(t\right)\right)-1\right),$where ${{\displaystyle \int }}_{-1}^1{\psi }_j\left(r\right)f_R\left(r\right)dr=E_R\left[{\psi }_j\left(r\right)\right]$ is the mean of the sample approximates the expression. Thus, equality can be written as

$$I\left(\tau \right)={{\displaystyle \int }}_{-1}^1{f}_R^2\left(r\right)dr-1/2={{\displaystyle \sum }}_{j=1}^s\frac{2j+1}{2}{E}_R^2\left[{\psi }_j\left(r\right)\right].$$

(22)

It should be noted that the infinite amount has been changed to finite. Such a change has advantages: the sum is calculated faster, giving robustness to the projection index. By summing only a finite number of members, the slowly fading “tails” of the projection distributions have a more negligible effect on the value of the projection index. Therefore, it is suggested to choose 4 ≤ s ≤ 7.

There are many methods for finding “interesting” projections. The method used in this research for finding the ‘most interesting’ projection direction is a mixed optimization strategy [55,64,84]. After defining the analytical expression of the projection index, its gradient in the projection direction τ is obtained as follows

$$\frac{\partial I}{\partial \tau }=\frac{2}{\sqrt{2\pi }}{{\displaystyle \sum }}_{j=1}^s\left(2j+1\right)E\left[{\psi }_j\left(r\right)\left]E\right[{\psi }_j^{\prime }\left(r\right)e^{-y^2/2}\left(z-\tau y\right)\right].$$

(23)

Here, the Lagrangian polynomial derivative is calculated by an iterative formula: ${\psi }_1^{\prime }\left(r\right)=1$, then ${\psi }_j^{\prime }\left(r\right)=r{\psi }_{j-1}^{\prime }\left(r\right)+j{\psi }_{j-1}\left(r\right)$, then $j\ge 1$. Initially, an approximate step optimizer is found by searching in the directions of the main components and their combinations so that the initial convergence to the maximum can be achieved quickly. Then, the approximate step optimizer (steepest ascent) quickly selects the projections required to ascend to the (local) maximum of the projection index. The projection index is used to search for ‘interesting’ data projections. However, it is usually not enough to find a single projection to reasonably estimate the multidimensional density. In general, “interesting” directions do not have to be orthogonal and may require more projection directions than the data dimension. Therefore, when estimating density by targeted projection, the so-called deletion of the data structure is applied. A nonlinear scale transformation is performed, found in the projection direction, so the distribution of the transformed data becomes normal. This operation ensures that the same direction as before was not found when searching for another projection direction.

The deletion of the data structure is based on the fact that if the projection of one-dimensional data projection ${\tau }^{\prime }Z$ has a distribution density $f_{\tau }\left(y\right)$ and a corresponding distribution function $F_{\tau }$, then the random variable is equal to

$$\tilde{Y}={\mathsf{\Phi }}^{-1}\left(F_{\tau }\left(Y\right)\right),$$

(24)

where ${\mathsf{\Phi }}^{-1}$ is the inverse of the standard normal distribution. Friedman [55] proposed to calculate the empirical estimate of the distribution function as follows ${\widehat{F}}_{\tau }\left(y\right)=rank\left(Y\right)/n-\frac{1}{2n}$, where $rank\left(y\right)$ is the rank of $Y$ in the whole sample of size n. Unfortunately, this estimate is not accurate and often results in a very uneven density function. By denoting $Z^{\left(0\right)}=Z$, we will discuss how $Z^{\left(k-1\right)}$ is obtained from $Z^{\left(k\right)}$. Based on Equation (14), $Z^{\left(k\right)}$ can be defined as

$$Z^{\left(k\right)}=Z^{\left(k-1\right)}+\left[{\Phi }^{-1}\left(F_{\tau }\left({\tau }^{\prime }{Z}^{\left(k-1\right)}\right)\right)-{\tau }^{\prime }{Z}^{\left(k-1\right)}\right]\tau .$$

(25)

The same procedure is performed to find the ‘most interesting’ projection with $Z^{\left(k\right)}$ searching for a new direction. This sequence is repeated until the multidimensional distribution becomes close to the Gaussian distribution in all directions. It has been observed [55] that gaussianization in one direction disrupts normalcy in the directions previously studied so that their projection index $I\left(\tau \right)$ is no longer zero. However, studies show [54] that the changes in results are minimal. Multidimensional density is calculated from one-dimensional density estimates.

The relationship between the multidimensional densities $Z^{\left(k\right)}$ and $Z^{\left(k-1\right)}$ (where $Z^{\left(k\right)}$ is the structure of the distant data $Z^{\left(k-1\right)}$ along the k-th projection $\tau \left(k\right)$) is $f_{\tau \left(k\right)}\left(z^{\left(k\right)}\right)=\frac{f_{\tau \left(k-1\right)}\left(z^{\left(k-1\right)}\right)}{\left|J_k\left(z^{\left(k-1\right)}\right)\right|}$ and $f_{\tau \left(k-1\right)}\left(z^{\left(k-1\right)}\right)=f_{\tau \left(k\right)}\left(z^{\left(k\right)}\right)\left|J_k\left(z^{\left(k-1\right)}\right)\right|$, here is the Jacobian $J_k\left(z^{\left(k-1\right)}\right)=\frac{\partial z^{\left(k\right)}}{\partial z^{\left(k-1\right)}}=\frac{\partial \left(U{z}^{\left(k\right)}\right)}{\partial \left(U{z}^{\left(k-1\right)}\right)}=\frac{\partial y^{\left(k\right)}}{\partial y^{\left(k-1\right)}}=\frac{f_{\tau \left(k\right)}\left(y^{\left(k-1\right)}\right)}{\varphi \left(y^{\left(k\right)}\right)}=\frac{f_{\tau \left(k\right)}\left({\tau }^{\prime }\left(k\right)z^{\left(k-1\right)}\right)}{\varphi \left({\tau }^{\prime }\left(k\right)z^{\left(k\right)}\right)}\ge 0$.

Starting from the initial multidimensional data $Z^{\left(0\right)}$gaussianization procedure is performed for each “interesting” projection found by $I\left(\tau \right)$. After a certain number, the projections’ multidimensional data $Z^{\left(M\right)}$ differ slightly from the normal distribution. Density $Z^{\left(0\right)}$ can be calculated as follows

$$\begin{array}{c}f\left(z^{\left(0\right)}\right)=f_{\tau \left(1\right)}\left(z^{\left(1\right)}\right)J_1\left(z^{\left(0\right)}\right)=f_{\tau \left(2\right)}\left(z^{\left(2\right)}\right)J_2\left(z^{\left(1\right)}\right)J_1\left(z^{\left(0\right)}\right)=f_{\tau \left(M\right)}\left(z^{\left(M\right)}\right){\displaystyle {\displaystyle \prod }_{k=1}^M}{J}_k\left(z^{\left(k-1\right)}\right)\\ \approx \varphi \left(z^{\left(M\right)}\right){\displaystyle {\displaystyle \prod }_{k=1}^M}{J}_k\left(z^{\left(k-1\right)}\right)=\varphi \left(z^{\left(M\right)}\right){\displaystyle {\displaystyle \prod }_{k=1}^M}\frac{f_{\tau \left(k\right)}\left({\tau }^{\prime }\left(k\right)z^{\left(k-1\right)}\right)}{\varphi \left({\tau }^{\prime }\left(k\right)z^{\left(k\right)}\right)}.\end{array}$$

(26)

The one-dimensional density of the projected data $f_{\tau \left(k\right)}\left({\tau }^{\prime }\left(k\right)z^{\left(k-1\right)}\right)$ is calculated according to Equation (18) or more precisely

$$f_{\tau \left(k\right)}\left({\tau }^{\prime }\left(k\right)z^{\left(k-1\right)}\right)=\varphi \left({\tau }^{\prime }\left(k\right)z^{\left(k-1\right)}\right){\displaystyle \sum }_{j=0}^s\frac{2j+1}{n}{\displaystyle \sum }_{t=1}^n{\psi }_j\left(r_t^{\left(k-1\right)}\right){\psi }_j\left(r^{\left(k-1\right)}\right).$$

(27)

Then, replacing the unknown one-dimensional distribution densities on the right-hand side (26) with their statistical estimates, we obtain

$$\widehat{f}\left(z\right)=\varphi \left(z^{\left(M\right)}\right){\displaystyle \prod }_{k=1}^M\frac{{\widehat{f}}_{\tau \left(k\right)}\left({\tau }^{\prime }\left(k\right)z^{\left(k-1\right)}\right)}{\varphi \left({\tau }^{\prime }\left(k\right)z^{\left(k\right)}\right)}.$$

(28)

The target projection density estimate is calculated relatively quickly because of the shape of the multivariate projection index and the iterative relationship between polynomials.

2.5. Inversion Formula

When examining approximations of parametric methods, it should be emphasized that as the data dimension increases, the number of model parameters increases rapidly, making it more difficult to find accurate parameter estimates. One-dimensional data projections $X_{\tau }={\tau }^{\prime }X$ density $f_{\tau }$ is much easier to find than multidimensional data density $f$ because there exists a mutually unambiguous correspondence, $f\leftrightarrow \left\{f_{\tau }, \tau \in R^d\right\}$. It is quite natural to try to find the multidimensional density $f$ using the density estimates ${\widehat{f}}_{\tau }$ of one-dimensional observational projections. It should be noted that in the case of the mixture, when the distributions are Gaussian, the projections of observations are also distributed according to the (one-dimensional) Gaussian mixture model

$$f_{\tau }\left(x\right)={{\displaystyle \sum }}_{k=1}^q{p}_{k,\tau }{\phi }_{k,\tau }\left(x\right)=f_{\tau }\left(x,{\theta }_{\tau }\right).$$

(29)

Here ${\phi }_{k,\tau }\left(x\right)=\phi \left(x; m_{k,\tau }, {\sigma }_{k,\tau }^2\right)$– one-dimensional Gaussian density. The parameter $\theta$ of the multidimensional mixture. The distribution parameters of the data projections ${\theta }_{\tau }=\left(p_{k,\tau }, m_{k,\tau }, {\sigma }_{k,\tau }^2\right)$, $k=1,\dots ,q$ are related by equations: $p_{j,\tau }=p_j, m_{j,\tau }={\tau }^{\prime }{M}_j \mathrm{and} {\sigma }_{j,\tau }^2={\tau }^{\prime }{R}_j\tau$. Using the inversion formula

$$f\left(x\right)=\frac{1}{{\left(2\pi \right)}^d}\underset{R^d}{\overset{}{{\displaystyle \int }}}{e}^{-i{t}^{\prime }x}\psi \left(t\right)dt,$$

(30)

where $\psi \left(t\right)=E{e}^{i{t}^{\prime }x}$ denotes the characteristic function of the random variable $X$. Marking $u=\left|t\right|$, $\tau =t/\left|t\right|$ and changing the variables to a spherical coordinate system, density is written

$$f\left(x\right)=\frac{1}{{\left(2\pi \right)}^d}{{\displaystyle \int }}_{\tau : \left|\tau \right|=1}^{}ds{{\displaystyle \int }}_0^{\infty }{e}^{-iu{\tau }^{\prime }x}\psi \left(u\tau \right)u^{d-1}du.$$

(31)

Here, the first integral is understood as the surface integral of the unit sphere. After noting the characteristic function of the projection of the observed random variable as ${\psi }_{\tau }\left(u\right)=E{e}^{iu{\tau }^{\prime }X}$. Then equality $\psi \left(u\tau \right)={\psi }_{\tau }\left(u\right)$ holds. By selecting the set $T$ of projection directions evenly spaced on the sphere and replacing the characteristic function with its estimate ($\widehat{f}\left(x\right))$a formula

$$\widehat{f}\left(x\right)=\frac{A\left(d\right)}{\#T}{\displaystyle \sum }_{\tau \in T}\underset{0}{\overset{\infty }{{\displaystyle \int }}}{e}^{-iu{\tau }^{\prime }x}{\widehat{\psi }}_{\tau }\left(u\right)u^{d-1}{e}^{-h{u}^2}du,$$

(32)

is obtained to calculate the estimate [85,86]. Here and # continue to denote the number of elements in the set $T$. Using the d-meter ball volume

$$V_d\left(R\right)=\frac{{\pi }^{\frac{d}{2}}{R}^d}{\mathsf{\Gamma }\left(\frac{d}{2}+1\right)}=\left\{\begin{array}{cc}\frac{{\pi }^{\frac{d}{2}}{R}^d}{\left(\frac{d}{2}\right)!},& \mathrm{then} d \mathrm{mod} 2\equiv 0\\ \frac{2^{\frac{d+1}{2}}{\pi }^{\frac{d-1}{2}}{R}^d}{d!!},& \mathrm{then} d \mathrm{mod} 2\equiv 1\end{array}\right.,$$

(33)

the constant $A\left(d\right)$ depending on the data dimension can be calculated using

$$A\left(d\right)=\frac{{\left(V_d\left(1\right)\right)}_R^{\prime }}{{\left(2\pi \right)}^d}=\frac{d{2}^{-d}{\pi }^{-\frac{d}{2}}}{\mathsf{\Gamma }\left(\frac{d}{2}+1\right)}.$$

(34)

Computer simulation studies have shown that the density estimates obtained using the inversion formula are not smooth. Therefore, in formula (32), an additional multiplier $e^{-h{u}^2}$ is used below the integral sign. This multiplier further smoothes the estimate $\widehat{f}\left(x\right)$ (32) with the Gaussian kernel function. This form of the multiplier allows the value of the integral to be calculated analytically. The number of clusters and Gaussian mixture parameters was selected using the constructive procedure and software developed at the Lithuanian Institute of Mathematics and Informatics, applying the w² type criterion [87]. Formula (32) can be used for various estimates of the characteristic function of the projected data. We will discuss the two methods used in this work.

One of them is based on the density approximation of the Gaussian distribution mixture model. In the present case, after replacing the parameters of the Gaussian mixture with their statistical estimates (${\widehat{p}}_{k,\tau }=p_k, {\widehat{m}}_{k,\tau }={\tau }^{\prime }{M}_k, {\sigma }_{k,\tau }^2={\tau }^{\prime }{R}_k\tau )$ (Page 10), the following parametric estimate

$${\widehat{\psi }}_{\tau }\left(u\right)={\displaystyle \sum }_{k=1}^{{\widehat{q}}_{\tau }}{\widehat{p}}_{k,\tau }{e}^{iu{\widehat{m}}_{k,\tau }-u^2{\widehat{\sigma }}_{k,\tau }^2/2}$$

(35)

of the characteristic function is used, and adding (32) to (35) gives

$$\begin{array}{c}\widehat{f}\left(x\right)=\frac{A\left(d\right)}{\#T}{\displaystyle \sum _{\tau \in T}}{\displaystyle \sum _{k=1}^{{\widehat{q}}_{\tau }}}{\widehat{p}}_{k,\tau }\underset{0}{\overset{\infty }{{\displaystyle \int }}}{e}^{iu\left({\widehat{m}}_{k,\tau }-{\tau }^{\prime }x\right)-u^2\left(h+{\widehat{\sigma }}_{k,\tau }^2/2\right)}{u}^{d-1}du\\ =\frac{A\left(d\right)}{\#T}{\displaystyle \sum _{\tau \in T}}{\displaystyle \sum _{k=1}^{{\widehat{q}}_{\tau }}}{\widehat{p}}_{k,\tau }{I}_{d-1}\left(\frac{{\widehat{m}}_{k,\tau }-{\tau }^{\prime }x}{\sqrt{{\widehat{\sigma }}_{k,\tau }^2+2h}}\right){\left(\sqrt{{\widehat{\sigma }}_{k,\tau }^2+2h}\right)}^{-d}\end{array}$$

(36)

and where $I_j\left(y\right)$ can be written as

$$I_j\left(y\right)=\mathrm{Re}\left[\underset{0}{\overset{\infty }{{\displaystyle \int }}}{e}^{iyz-z^2/2}{z}^{j}dz\right].$$

(37)

It should be noted that only the real part of the expression can be considered here. The sum of the imaginary parts must be equal to zero. Because the density estimate $\widehat{f}\left(x\right)$ can acquire only real values. The chosen form of the smoothing multiplier $e^{-h{u}^2}$ allows relating the smoothing parameter h to the variances of the projection clusters—in the calculations, the variances are increased by $2h$. How to calculate expression (37) is given in Appendix B.

2.6. Modified Density Estimate of the Inversion Formula

One of the disadvantages of the inversion formula method defined in (32) is that the Gaussian distribution mixture model described by this estimate (where $f_k={\phi }_k$) evaluates well only the density of observations close to it. However, when approximating the density under study with a mixture of Gaussian distributions, the estimation of the density of the inversion formula often becomes complicated due to a large number of components with low a priori probabilities. Their number can be reduced by introducing a noise cluster—the modified algorithm based on a multidimensional Gaussian distribution mixture model. Let us use the inversion formula (30). The parametric estimate of the characteristic function of uniform distribution density can be calculated as follows

$$\widehat{\psi }\left(u\right)=\frac{2}{\left(b-a\right)u}\sin \frac{\left(b-a\right)u}{2}·e^{\frac{iu\left(a+b\right)}{2}}.$$

(38)

In uniform distribution density function (38), b is the maximum value, and a is the minimum value. In the density estimate calculation formula (32), construct the estimation of the characteristic function as a union of the characteristic functions of a mixture of Gaussian distributions and a uniform distribution with corresponding a priori probabilities as follows

$${\widehat{\psi }}_{\tau }\left(u\right)={\displaystyle \sum }_{k=1}^{{\widehat{q}}_{\tau }}{\widehat{p}}_{k,\tau }{e}^{iu{\widehat{m}}_{k,\tau }-u^2{\widehat{\sigma }}_{k,\tau }^2/2}+{\widehat{p}}_{0,\tau }\frac{2}{\left(b\left(\tau \right)-a\left(\tau \right)\right)u}\sin \frac{\left(b\left(\tau \right)-a\left(\tau \right)\right)u}{2}·e^{\frac{iu\left(a\left(\tau \right)+b\left(\tau \right)\right)}{2}}.$$

(39)

Here the second term describes a uniform distributed noise cluster and ${\widehat{p}}_0$ is the weight of the noise cluster. Based on the parameters of the uniform distribution and the projected data, we can write

$$a\left(\tau \right)={\left({\tau }^{\prime }x\right)}_{\min }-\frac{{\left({\tau }^{\prime }x\right)}_{\max }-{\left({\tau }^{\prime }x\right)}_{\min }}{2\left(n-1\right)} \mathrm{and}$$

(40)

$$b\left(\tau \right)={\left({\tau }^{\prime }x\right)}_{\max }+\frac{{\left({\tau }^{\prime }x\right)}_{\max }-{\left({\tau }^{\prime }x\right)}_{\min }}{2\left(n-1\right)}.$$

(41)

Using notations such as (36), we can write

$$\begin{array}{c}\widehat{f}\left(x\right)=\frac{A\left(d\right)}{\#T}{\displaystyle \sum _{\tau \in T}}\left[{\displaystyle \sum _{k=1}^{{\widehat{q}}_{\tau }}}{\widehat{p}}_{k,\tau }{I}_{d-1}\left(\frac{{\widehat{m}}_{k,\tau }-{\tau }^{\prime }x}{\sqrt{{\widehat{\sigma }}_{k,\tau }^2+2h}}\right){\left({\widehat{\sigma }}_{k,\tau }^2+2h\right)}^{-\frac{d}{2}}\right.\\ \left.+\frac{2{\widehat{p}}_{0,\tau }}{b\left(\tau \right)-a\left(\tau \right)}{J}_{d-2}\left(\frac{a\left(\tau \right)+b\left(\tau \right)-2{\tau }^{\prime }x}{2\sqrt{2h}},\frac{b\left(\tau \right)-a\left(\tau \right)}{2\sqrt{2h}}\right)·{\left(2h\right)}^{-\frac{d-1}{2}}\right].\end{array}$$

(42)

where the expression $I_j\left(y\right)$ is the same as (37) and its value is $I_j\left(y\right)=C_j\left(y\right)$and

$$J_j\left(y,z\right)=\mathrm{Re}\left[\underset{0}{\overset{\infty }{{\displaystyle \int }}}{e}^{iyu-u^2/2}·\sin zu·u^{j}du\right].$$

(43)

By integrating, we get

$$\begin{array}{c}\underset{0}{\overset{\infty }{{\displaystyle \int }}}{e}^{iyu-u^2/2}·\sin zu·u^{j}du=\underset{0}{\overset{\infty }{{\displaystyle \int }}}\left(\cos yu+i\sin yu\right)·\sin zu·e^{-u^2/2}·u^{j}du\\ =\underset{0}{\overset{\infty }{{\displaystyle \int }}}\left(\frac{\sin \left(y+z\right)u+\sin \left(z-y\right)u}{2}+i\frac{\cos \left(y-z\right)u-\cos \left(y+z\right)u}{2}\right)·e^{-u^2/2}·u^{j}du\\ =\frac{1}{2}{S}_j\left(y+z\right)+\frac{1}{2}{S}_j\left(z-y\right)+i\frac{1}{2}{C}_j\left(y-z\right)-i\frac{1}{2}{C}_j\left(y+z\right).\end{array}$$

(44)

the above formula uses the variables $S_j\left(y\right)$ and $C_j\left(y\right)$, the calculation of which is given in formulas (52) and (53) in Appendix B.

3. Materials and Methods

Density estimation algorithms were presented in the previous section. The Monte Carlo method was used in this study. Such a comparison of algorithms allows us to measure the real observation density values and evaluate algorithms’ efficiency. For the research, we used multidimensional (d = 2, 5, 10, 15) distributions of the Cauchy mixture

$${{\displaystyle \sum }}_{j=1}^q{p}_{j}C\left(x,m_j,u_j\right).$$

(45)

Additionally, $C\left(x,m_j,u_j\right)$ is defined as follows

$$C\left(x,m_j,u_j\right)={{\displaystyle \prod }}_{k=1}^d\frac{u_{jk}}{\pi [u_{jk}^2+\left(x_k-m_{jk}{)}^2\right]}.$$

(46)

Calculations were performed using sample sizes of n = 50, 100, 200, 400, 800 while changing the number of distributions, their weights, and centers (see.

Table 1. Parameters table.

Number of Components	Proportions of Components	Location Parameters	Separation Size of Locations
q = 2	p1 = (1 − p2), p2 = 0.1, 0.3, 0.5	m1 = (0, 0), m2 = (0.5i, 0.5i)	i = 1, 2,…, 6
q = 3	p1 = p2 = (1 − p3)/2, p3 = 0.1, 1/3, 0.8	m1 = (0, 0), m2 = (0.5i, 0.5i), m3 = (0.5i, 0)	i = 1, 2,…, 6
q = 4	p1 = p2 = p3 = (1 − p4)/3, p4 = 0.1, 0.25, 0.7	m1 = (0, 0), m2 = (0.5i, 0.5i), m3 = (0.5i, 0), m4 = (0, 0.5i)	i = 1, 2,…, 6
q = 2	p1 = (1 − p2), p2 = 0.1, 0.2, 0.3, 0.4, 0.5	m1 = (0, 0, 0, 0, 0), m2 = (0.5i, 0.5i, 0.5i, 0.5i, 0.5i)	i = 1, 2,…, 6
q = 3	p1 = p2 = (1 − p3)/2, p3 = 0.1, 0.2, 1/3, 0.4, 0.6, 0.8	m1 = (0, 0, 0, 0, 0), m2 = (0.5i, 0.5i, 0.5i, 0.5i, 0.5i), m3 = (0.5i, 0.5i, 0, 0, 0)	i = 1, 2,…, 6
q = 4	p1 = p2 = p3 = (1 − p4)/3, p4 = 0.1, 0.16, 0.25, 0.4, 0.7	m1 = (0, 0, 0, 0, 0), m2 = (0.5i, 0.5i, 0.5i, 0.5i, 0.5i), m3 = (0.5i, 0.5i, 0, 0, 0), m4 = (0, 0, 0.5i, 0.5i, 0.5i)	i = 1, 2,…, 6

). In each case, 100,000 samples were generated.

In cases of d = 10, 15, the same weights were used as in d = 5. Additionally, centres are located on the apexes of the hypercube.

Algorithms used in the research: AKDE—adaptive kernel, PPDE—projection pursuit, LSDE—logspline, SKDE—semi-parametric kernel, IFDE—inversion formula, MIDE—inversion formula with noise cluster. In IFDE and MIDE methods are used mixture parameters, calculated with a program made in an institute of Mathematics and Informatics (Vilnius) [87].

Selection of parameters in the density estimation procedure. In this study, the Monte Carlo method aimed to perform the accuracy of the non-parametric estimates of distribution density previously described in the methodological sections (AKDE, PPDE, LSDE, SKDE, IFDE, MIDE) comparative analysis. The authors [34] propose to collect the value of the sensitivity parameter (γ, see. AKDE method step 3) used in the AKDE method from the set {0.2; 0.4; 0.6; 0.8}. The specific value of the parameter is determined by a probabilistic cross-check [88,89]. In the SKDE, all possible values of the sub-vector Y dimension s (1≤ s ≤ d − 1, where d is dimensions, see page 5) and their corresponding coordinates were reselected. The most factual errors were used to compare the results with other studied methods. The LSDE method minimizes the Akaike information criterion by selecting the number of baseline spline points [78]. The computer program for calculating this estimate is provided in the R package and was used in the study. Akaike information criterion $AIC=-2l\left(t\right)+aJ\left(t\right)$, $J$—degree of spline, $a=\log \left(n\right)$, $l$—probability function used to select the spline coefficients. The MIDE method has a smoothing parameter, h. The chosen form of the smoothing multiplier $e^{-h{u}^2}$ allows relating the smoothing parameter h to the variances of the projection clusters. Modelling studies have shown that this method is sensitive to parameter selection. If h is set too low, the estimate becomes very slick and has large errors. Excessive smoothing of the density estimate does not greatly affect its quality. In the studies, it was observed that the estimation becomes uneven due to the similarity of the values of the observations projected in some directions, thus distinguishing low-weight components with small dispersions. The smoothing parameter (h) as well as the specific value of the noise cluster weight (probability) from the set {0.05; 0.1; 0.15; 0.2; 0.3; 0.4} are selected by cross-checking the least squares [65]. The vector of the estimate parameters is searched for in such a way that it minimizes the integrated square error

$$\mathsf{\Theta }=\arg \underset{\mathsf{\Theta }}{\min }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}{\left(\widehat{f_{\mathsf{\Theta }}}\left(x\right)-f\left(x\right)\right)}^{2}dx=arg\underset{\mathsf{\Theta }}{\min }\left\{\widehat{\Vert f_{\mathsf{\Theta }}}{\left(x\right)\Vert }_2^2-\frac{2}{n}{\displaystyle \sum }_{t=1}^n\widehat{f_{\mathsf{\Theta }}}\left(X\left(t\right)\right)\right\},$$

(47)

where Θ is the evaluated parameter and F(x) is the observed random variable distribution function. Changing an unknown distribution function to an empirical distribution function yields an expression for the parameter estimate

$$\mathsf{\Theta } ^=\arg \underset{\mathsf{\Theta }}{\min }(\widehat{\Vert f_{\mathsf{\Theta }}}{\left(x\right)\Vert }_2^2-\frac{2}{n}{\displaystyle \sum }_{t=1}^n\widehat{f_{\mathsf{\Theta }}}(X\left(t\right)|t),$$

(48)

where $\widehat{f_{\mathsf{\Theta }}}(x|t)$ is the value of the estimate at point x, which is calculated by subtracting the value of X(t) from the observations. In addition, empirical research suggests that it is better to look for a maximum local minimum point rather than a global minimum [90]. Using PPDE method and following the recommendation of the paper [38], the order of the spread was $4\le s\le 6$ (see Page 9), and the projection directions were chosen to maximize the value of the estimate of the design index (2) recommended by J. H. Friedman

$$I\left(\alpha \right)={{\displaystyle \int }}_{-1}^1{f}_r^2\left(r\right)dr-\frac{1}{2}={{\displaystyle \sum }}_{j=1}^J\frac{2j+1}{2}{E}_r^2\left[{\psi }_j\left(r\right)\right].$$

(49)

4. Results and Discussion

This section presents the main results obtained during the simulations. We calculate the mean absolute error and (50) mean absolute percentage error (51) to evaluate the accuracy.

$${\delta }_1=\frac{1}{n}{{\displaystyle \sum }}_{t=1}^n\left|f\left(x\left(t\right)\right)-\widehat{f}\left(x\left(t\right)\right)\right|\cong \int \left|f\left(x\right)-\widehat{f}\left(x\right)\right|f\left(x\right)dx.$$

(50)

$${\delta }_2=\frac{2}{n}{{\displaystyle \sum }}_{t=1}^n\left|\frac{f\left(x\left(t\right)\right)-\widehat{f}\left(x\left(t\right)\right)}{f\left(x\left(t\right)\right)+\widehat{f}\left(x\left(t\right)\right)}\right|\cong \int \left|f\left(x\right)-\widehat{f}\left(x\right)\right|dx.$$

(51)

The result tables (

Table 2. An example of mean absolute percentage error.

Evaluation Methods		Density
		d = 5; p1 = p2 = p3 = 1/3; n = 100
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	0.8268	0.8257	0.8198	0.8128	0.8066	0.8075
	SD	(0.0760)	(0.0814)	(0.0848)	(0.0827)	(0.0788)	(0.0731)
PPDE	Mean	0.9243	0.9319	0.9303	0.9300	0.9284	0.9250
	SD	(0.0500)	(0.0364)	(0.0375)	(0.0387)	(0.0410)	(0.0433)
LSDE	Mean	0.8043	0.8162	0.8583	0.8611	0.8613	0.8711
	SD	(0.0534)	(0.0540)	(0.0490)	(0.0349)	(0.0434)	(0.0577)
SKDE	Mean	0.7158	0.7144	0.7088	0.7071	0.7179	0.7227
	SD	(0.0260)	(0.0905)	(0.0905)	(0.0830)	(0.0631)	(0.0499)
IFDE	Mean	0.94593	0.8886	0.7857	0.8463	0.8761	0.8312
	SD	(0.0362)	(0.1318)	(0.0706)	(0.0380)	(0.1110)	(0.0538)
MIDE	Mean	0.7389	0.7332	0.7235	0.7149	0.7121	0.7219
	SD	(0.0280)	(0.0221)	(0.0338)	(0.0195)	(0.0208)	(0.0203)

and

Table A1. 10² times upscaled mean absolute error with d = 2; p₁ = 0.5; p₂ = 0.5; n = 100.

Evaluation Methods		Density
		d = 2; p1 = 0.5; p2 = 0.5; n = 100
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	1.74	1.41	1.16	1.08	1.02	0.99
	SD	(0.94)	(0.80)	(0.70)	(0.59)	(0.54)	(0.50)
PPDE	Mean	2.21	1.85	1.52	1.32	1.25	1.21
	SD	(0.49)	(0.41)	(0.40)	(0.44)	(0.37)	(0.31)
LSDE	Mean	0.87	0.71	0.78	0.63	0.69	0.69
	SD	(0.43)	(0.20)	(0.08)	(0.08)	(0.04)	(0.09)
SKDE	Mean	0.63	0.61	0.52	0.52	0.51	0.51
	SD	(0.12)	(0.17)	(0.07)	(0.05)	(0.06)	(0.04)
IFDE	Mean	1.69	1.31	0.97	0.75	0.61	0.53
	SD	(0.06)	(0.10)	(0.08)	(0.01)	(0.04)	(0.06)
MIDE	Mean	0.69	0.66	0.57	0.55	0.51	0.51
	SD	(0.06)	(0.10)	(0.08)	(0.01)	(0.04)	(0.06)

,

Table A2. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.7; p₂ = 0.3; n = 200.

Evaluation Methods		Density
		d=5; p1 = 0.7; p2 = 0.3; n = 200
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	0.979	0.801	0.703	0.664	0.655	0.654
	SD	(0.171)	(0.146)	(0.145)	(0.155)	(0.158)	(0.159)
PPDE	Mean	1.001	0.822	0.722	0.681	0.671	0.669
	SD	(0.174)	(0.152)	(0.151)	(0.160)	(0.163)	(0.163)
LSDE	Mean	5.039	4.185	2.632	0.944	0.665	0.660
	SD	(1.265)	(6.747)	(1.081)	(0.138)	(0.140)	(0.112)
SKDE	Mean	0.857	0.759	0.705	0.658	0.649	0.638
	SD	(0.087)	(0.069)	(0.076)	(0.085)	(0.083)	(0.083)
IFDE	Mean	0.912	0.801	0.721	0.681	0.667	0.666
	SD	(0.133)	(0.149)	(0.151)	(0.160)	(0.162)	(0.163)
MIDE	Mean	0.956	0.788	0.694	0.661	0.657	0.640
	SD	(0.162)	(0.154)	(0.144)	(0.152)	(0.158)	(0.163)

,

Table A3. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.45; p₂ = 0.45; p₃ = 0.1; n= 200.

Evaluation Methods		Density
		d = 5; p1 = 0.45; p2 = 0.45; p3 = 0.1; n = 200
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	0.970	0.725	0.576	0.533	0.504	0.500
	SD	(0.127)	(0.089)	(0.078)	(0.073)	(0.069)	(0.066)
PPDE	Mean	0.992	0.746	0.594	0.528	0.506	0.500
	SD	(0.137)	(0.092)	(0.077)	(0.072)	(0.069)	(0.066)
LSDE	Mean	1.057	0.775	0.652	0.590	0.508	0.503
	SD	(0.164)	(0.203)	(0.650)	(0.491)	(0.067)	(0.081)
SKDE	Mean	0.6245	0.6274	0.6312	0.629	0.630	0.628
	SD	(0.072)	(0.025)	(0.027)	(0.049)	(0.049)	(0.050)
IFDE	Mean	0.990	0.743	0.589	0.525	0.497	0.499
	SD	(0.136)	(0.091)	(0.076)	(0.071)	(0.071)	(0.067)
MIDE	Mean	0.993	0.746	0.574	0.525	0.496	0.490
	SD	(0.137)	(0.092)	(0.077)	(0.072)	(0.069)	(0.066)

,

Table A4. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.4; p₂ = 0.4; p₃ = 0.2; n = 400.

Evaluation Methods		Density
		d = 5; p1 = 0.4; p2 = 0.4; p3 = 0.2; n = 400
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	0.916	0.689	0.527	0.445	0.410	0.396
	SD	(0.099)	(0.068)	(0.048)	(0.044)	(0.048)	(0.052)
PPDE	Mean	0.937	0.709	0.545	0.461	0.423	0.407
	SD	(0.109)	(0.074)	(0.049)	(0.044)	(0.048)	(0.052)
LSDE	Mean	0.815	0.549	0.511	0.443	0.404	0.401
	SD	(0.007)	(0.063)	(0.151)	(0.094)	(0.040)	(0.030)
SKDE	Mean	0.655	0.499	0.413	0.388	0.385	0.384
	SD	(0.064)	(0.049)	(0.034)	(0.031)	(0.028)	(0.027)
IFDE	Mean	0.937	0.709	0.544	0.460	0.423	0.404
	SD	(0.109)	(0.074)	(0.049)	(0.044)	(0.048)	(0.052)
MIDE	Mean	0.757	0.509	0.415	0.391	0.391	0.388
	SD	(0.109)	(0.074)	(0.049)	(0.044)	(0.048)	(0.052)

,

Table A5. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.25; p₂ = 0.25; p₃ = 0.25; p₄ = 0.25; n = 400.

Evaluation Methods		Density
		d = 5; p1 = 0.25; p2 = 0.25; p3 = 0.25; p4 = 0.25; n = 400
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	0.912	0.645	0.447	0.351	0.309	0.290
	SD	(0.128)	(0.068)	(0.029)	(0.019)	(0.026)	(0.030)
PPDE	Mean	0.934	0.665	0.464	0.365	0.321	0.299
	SD	(0.145)	(0.089)	(0.048)	(0.031)	(0.033)	(0.035)
LSDE	Mean	0.934	0.676	0.464	0.365	0.321	0.293
	SD	(0.145)	(0.064)	(0.048)	(0.031)	(0.033)	(0.039)
SKDE	Mean	0.658	0.472	0.372	0.345	0.316	0.290
	SD	(0.071)	(0.031)	(0.020)	(0.017)	(0.019)	(0.018)
IFDE	Mean	0.933	0.665	0.464	0.365	0.321	0.299
	SD	(0.145)	(0.089)	(0.048)	(0.031)	(0.033)	(0.035)
MIDE	Mean	0.889	0.622	0.433	0.341	0.307	0.281
	SD	(0.118)	(0.074)	(0.037)	(0.026)	(0.019)	(0.027)

,

Table A6. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.1; p₂ = 0.1; p₃ = 0.1; p₄ = 0.7; n = 400.

EvaluationMethods		Density
		d = 5; p1 = 0.1; p2 = 0.1; p3 = 0.1; p4 = 0.7; n = 400
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	0.957	0.800	0.678	0.617	0.587	0.571
	SD	(0.131)	(0.127)	(0.103)	(0.090)	(0.087)	(0.087)
PPDE	Mean	0.979	0.821	0.697	0.634	0.603	0.586
	SD	(0.141)	(0.137)	(0.112)	(0.099)	(0.094)	(0.093)
LSDE	Mean	0.979	0.821	0.697	0.634	0.596	0.586
	SD	(0.141)	(0.137)	(0.112)	(0.099)	(0.098)	(0.093)
SKDE	Mean	0.687	0.580	0.514	0.496	0.491	0.489
	SD	(0.076)	(0.070)	(0.058)	(0.058)	(0.058)	(0.056)
IFDE	Mean	0.979	0.820	0.697	0.634	0.602	0.585
	SD	(0.141)	(0.137)	(0.112)	(0.098)	(0.094)	(0.093)
MIDE	Mean	0.924	0.770	0.652	0.597	0.533	0.488
	SD	(0.135)	(0.131)	(0.108)	(0.093)	(0.092)	(0.091)

,

Table A7. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.5; p₂ = 0.5; n = 50.

Evaluation Methods		Density
		d = 5; p1 = 0.5; p2 = 0.5; n = 50
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	1.093	0.828	0.758	0.741	0.739	0.740
	SD	(0.095)	(0.099)	(0.114)	(0.123)	(0.130)	(0.134)
PPDE	Mean	1.147	0.872	0.794	0.770	0.764	0.762
	SD	(0.157)	(0.150)	(0.157)	(0.156)	(0.159)	(0.160)
LSDE	Mean	2.100	1.997	2.002	2.010	2.013	2.014
	SD	(0.078)	(0.028)	(0.017)	(0.019)	(0.024)	(0.025)
SKDE	Mean	1.149	0.875	0.797	0.773	0.765	0.763
	SD	(0.160)	(0.154)	(0.160)	(0.160)	(0.161)	(0.162)
IFDE	Mean	1.145	0.864	0.780	0.765	0.763	0.757
	SD	(0.163)	(0.137)	(0.140)	(0.150)	(0.163)	(0.156)
MIDE	Mean	1.094	0.860	0.759	0.767	0.742	0.752
	SD	(0.142)	(0.167)	(0.160)	(0.156)	(0.163)	(0.160)

,

Table A8. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.6; p₂ = 0.4; n = 50.

Evaluation Methods		Density
		d = 5; p1 = 0.6; p2 = 0.4; n = 50
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	1.138	0.872	0.770	0.748	0.745	0.746
	SD	(0.105)	(0.085)	(0.126)	(0.143)	(0.152)	(0.155)
PPDE	Mean	1.192	0.918	0.808	0.778	0.771	0.769
	SD	(0.150)	(0.136)	(0.172)	(0.178)	(0.181)	(0.182)
LSDE	Mean	2.114	1.995	1.977	1.983	1.986	1.987
	SD	(0.101)	(0.083)	(0.080)	(0.079)	(0.082)	(0.084)
SKDE	Mean	1.195	0.919	0.810	0.780	0.772	0.770
	SD	(0.154)	(0.138)	(0.174)	(0.182)	(0.183)	(0.183)
IFDE	Mean	1.183	0.906	0.802	0.778	0.765	0.769
	SD	(0.142)	(0.125)	(0.163)	(0.185)	(0.176)	(0.185)
MIDE	Mean	1.152	0.882	0.782	0.754	0.747	0.742
	SD	(0.136)	(0.124)	(0.155)	(0.175)	(0.176)	(0.180)

,

Table A9. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.33; p₂ = 0.33; p₃ = 0.33; n = 50.

Evaluation Methods		Density
		d = 5; p1 = 0.33; p2 = 0.33; p3 = 0.33; n = 50
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	1.166	0.828	0.634	0.547	0.512	0.500
	SD	(0.120)	(0.086)	(0.057)	(0.064)	(0.074)	(0.080)
PPDE	Mean	1.224	0.879	0.677	0.581	0.540	0.523
	SD	(0.184)	(0.128)	(0.107)	(0.108)	(0.108)	(0.109)
LSDE	Mean	2.075	1.934	1.921	1.939	1.938	1.937
	SD	(0.127)	(0.094)	(0.058)	(0.054)	(0.048)	(0.044)
SKDE	Mean	1.226	0.881	0.678	0.583	0.542	0.524
	SD	(0.186)	(0.130)	(0.109)	(0.110)	(0.110)	(0.110)
IFDE	Mean	1.215	0.839	0.649	0.554	0.522	0.513
	SD	(0.175)	(0.099)	(0.110)	(0.102)	(0.110)	(0.111)
MIDE	Mean	1.182	0.834	0.638	0.545	0.518	0.501
	SD	(0.167)	(0.124)	(0.097)	(0.101)	(0.106)	(0.106)

and

Table A10. 10⁴ times upscaled mean absolute error with d = 5; p₁ = 0.45; p₂ = 0.45; p₃ = 0.1; n = 50.

Evaluation Methods		Density
		d = 5; p1 = 0.45; p2 = 0.45; p3 = 0.1; n = 50
		i = 1	i = 2	i = 3	i = 4	i = 5	i = 6
AKDE	Mean	1.126	0.838	0.690	0.633	0.618	0.615
	SD	(0.112)	(0.126)	(0.063)	(0.053)	(0.061)	(0.067)
PPDE	Mean	1.182	0.882	0.727	0.660	0.640	0.634
	SD	(0.156)	(0.132)	(0.091)	(0.085)	(0.087)	(0.088)
LSDE	Mean	2.101	2.002	1.985	2.010	2.014	2.015
	SD	(0.105)	(0.063)	(0.012)	(0.029)	(0.028)	(0.025)
SKDE	Mean	1.183	0.885	0.729	0.663	0.642	0.635
	SD	(0.157)	(0.134)	(0.094)	(0.089)	(0.090)	(0.090)
IFDE	Mean	1.170	0.859	0.702	0.649	0.624	0.619
	SD	(0.142)	(0.129)	(0.074)	(0.088)	(0.083)	(0.086)
MIDE	Mean	1.142	0.850	0.696	0.639	0.620	0.618
	SD	(0.141)	(0.125)	(0.080)	(0.084)	(0.083)	(0.087)

) provide 100,000 samples densities mean absolute percentage error. The values in parentheses provide information about the standard deviation of errors. The best results in these tables are bolded and underlined. According to Table 2, it is concluded that when n = 100, d = 5, the best results are obtained by SKDE and MIDE methods. Based on Table A2, it can be observed that when q = 2, n = 200, the best results are obtained using SKDE and MIDE methods. According to Table A3, it is concluded that when q = 3, n = 200, in the case of highly overlapping distributions (i = 1, 2), the best results are obtained by the SKDE method, and in the case of more isolated distributions (i ≥ 3)—by the MIDE method. Based on Table A4, it can be observed that when q = 3, n ≥ 400, the best results are obtained by SKDE, while the second-best method is MIDE. According to Table A5, it is concluded that when q = 4, n = 400, in the case of highly overlapping distributions (i ≤ 3), the best results are obtained by the SKDE method and in the case of more isolated distributions (i ≥ 4)—by the MIDE method. Table A6 shows results of q = 4, n ≥ 400, it can be noticed that, in the case of highly overlapping or average isolated distributions (i ≤ 5), the best results are obtained by the SKDE method and in the case of more isolated distributions (i = 6)—by the MIDE method. Table A7 and Table A8 show results of q = 2 and n = 50, It can be noticed that in all cases highly overlapping or isolated distributions, the best results are obtained by AKDE method and in the case of more isolated distributions (i = 6) with p₁ = 0.6; p₂ = 0.4— by the MIDE method. Table A9 and Table A10 show results q = 3 and n = 50; the best results are obtained by the AKDE method in all cases (highly overlapping or isolated distributions).The LSDE method with huge outliers (|x − m_j| > 100 uj) is grouped with a more significant number of values closer to the centre of the distribution. With the help of the calculated spline coefficients, the density in the outliers is estimated at a value close to 10¹⁰⁰. That is incorrect, and in such cases, the use of this method is not recommended.

The results for the smaller dimensions (d = 2) are presented in Table A1. It can be seen that the best results are obtained using the SKDE method, both in large- and small-scale overlapping cases (i < 4). On the other hand, in the case of isolated distributions (i ≥ 5), good results were obtained by the MIDE method.

In the case of mean absolute percentage error, recommended using the semiparametric kernel when the sample has overlapping distributions. In the case of two dimensions (d ~ 2) and a sample is with overlapping distributions, it is recommended to use the semiparametric kernel method and for isolated distributions, to use the adaptive kernel method.

5. Conclusions

This paper reviewed the most popular and most often used nonparametric density estimation algorithms. The density estimation inversion formula was also presented in this article. It was observed that when a noise cluster is included, the results of the inversion formula improved statistically significantly. Based on the mean absolute error, in the case of higher dimension (d ~ 5) and small samples (n ~ 50), it is recommended to use the adaptive kernel method. If the sample is n ~ 100, then the modified inversion formula method showed the best results. For larger samples with overlapping distributions it is recommended to use a semi-parametric kernel and for more isolated distribution—modified inversion methods. Based on the mean absolute percentage error, it is recommended to use the semiparametric kernel when the sample is with overlapping distributions. In the case of two dimensions (d ~ 2) and a sample is with overlapping distributions, it is recommended to use the semiparametric kernel method. For isolated distributions, it is recommended to use the adaptive kernel method.