Estimating densities with nonlinear support using Fisher-Gaussian kernels

TL;DR

This paper introduces Fisher-Gaussian kernels, a new class of density estimators that better capture data supported on nonlinear manifolds by incorporating curvature, improving performance over traditional Gaussian-based methods.

Contribution

It proposes a novel Fisher-Gaussian kernel class with an added curvature parameter, enabling more accurate density estimation on nonlinear supports.

Findings

Fisher-Gaussian kernels effectively model data on curved manifolds.

The method shows improved density estimation in simulations and real data.

Analytic form facilitates Bayesian implementation with MCMC.

Abstract

Current tools for multivariate density estimation struggle when the density is concentrated near a nonlinear subspace or manifold. Most approaches require choice of a kernel, with the multivariate Gaussian by far the most commonly used. Although heavy-tailed and skewed extensions have been proposed, such kernels cannot capture curvature in the support of the data. This leads to poor performance unless the sample size is very large relative to the dimension of the data. This article proposes a novel generalization of the Gaussian distribution, which includes an additional curvature parameter. We refer to the proposed class as Fisher-Gaussian (FG) kernels, since they arise by sampling from a von Mises-Fisher density on the sphere and adding Gaussian noise. The FG density has an analytic form, and is amenable to straightforward implementation within Bayesian mixture models using Markov chain Monte Carlo. We provide theory on large support, and illustrate gains relative to competitors in simulated and real data applications.

Figures

Tables

Table 1. Table 1: MAD δ for different choices of δ 𝛿 \delta for 3 competing methods in Euler Spiral dataset

$\delta$	KDE	DPMG	FG-mixture
0.005	0.928	0.964	0.924
0.010	2.816	3.242	1.924
0.015	4.512	5.624	2.428
0.020	5.836	7.628	2.930
0.025	6.886	9.906	3.646
0.030	7.714	11.524	4.118
0.035	8.202	13.084	4.660

Table 2. Table 3: The MAD δ for different choices of δ 𝛿 \delta for the ‘Torus’ dataset

$\delta$	KDE	DPMG	FG-mixture
0.1	0.121	0.096	0.117
0.2	0.626	0.552	0.625
0.3	1.255	1.381	1.204
0.4	2.127	2.653	1.982
0.5	2.829	3.991	2.611
0.6	3.703	5.754	3.121
0.7	5.515	11.071	5.340

Table 3. Table 5: Classification accuracies of the competing methods for the Gesture Phase data

D against P		Classification Accuracy			D against S		Classification Accuracy
Train size *	Test size *	KDE	DPMG	FG-mixture	Train size *	Test size *	KDE	DPMG	FG-mixture
50	100	**	0.490	0.585	100	100	0.475	0.690	0.700
100	63	0.508	0.603	0.706	200	100	0.495	0.740	0.700
D against R		Classification Accuracy			P against S		Classification Accuracy
50	100	**	0.455	0.655	50	100	**	0.505	0.590
100	91	0.500	0.725	0.736	100	63	0.500	0.548	0.579
P against R		Classification Accuracy			S against R		Classification Accuracy
50	100	**	0.530	0.585	50	100	0.495	0.475	0.550
100	63	0.516	0.627	0.635	100	91	0.532	0.595	0.611

Table 4. Table 6: Sum of log-likelihoods of test samples for competing methods in Balance scale data

Training size	Test size	KDE	DPMG	FG-mixture
50	100	-645.20	-611.72	-574.93
100	100	-625.78	-571.38	-549.56
150	100	-621.60	-574.87	-550.17
200	100	-605.33	-567.37	-551.23

Equations

$$f(\mathbf{x})=\sum _{k=1}^M{\pi }_k\mathcal{K}\left(\mathbf{x};{\mathbf{\theta }}_k\right),$$

$${\mathbf{x}}_i=\mathbf{c}+r{\mathbf{y}}_i+{\mathbf{ϵ}}_i\text{mbox}where{\mathbf{ϵ}}_i\sim N\left(0,{\sigma }^{2}I\right),$$

$$f(\mathbf{x})$$

$$={\int }_{{\mathbb{R}}^d}\frac{1}{r^d}{f}_{\mathrm{vMF}}\left\{\frac{1}{r}\left(\mathbf{z}-\mathbf{c}\right)\mid \mathbf{\mu },\tau \right\}{\varphi }_{{\sigma }^2}\left(\mathbf{x}-\mathbf{z}\right)d\mathbf{z},$$

$${\mathrm{FG}}_{\sigma }(\mathbf{x}\mid \mathbf{\rho })=\frac{C_d(\tau ){\left(2\pi {\sigma }^2\right)}^{-(d/2)}}{C_d\left(\Vert \tau \mathbf{\mu }+r(\mathbf{x}-\mathbf{c})/{\sigma }^2\Vert \right)}\exp \left\{-\frac{1}{2{\sigma }^2}\left(\Vert \mathbf{x}-\mathbf{c}{\Vert }^2+r^2\right)\right\},$$

$$\frac{1}{r_{s_i}}\left({\mathbf{z}}_i-{\mathbf{c}}_{s_i}\right)\sim f_{\mathrm{vMF}}\left(\cdot \mid {\mathbf{\mu }}_{k_i}^{s_i},{\tau }_{k_i}^{s_i}\right),$$

$$\left({\mathbf{x}}_i\mid s_i=l,k_i=k\right)\sim {\mathrm{FG}}_{\sigma }\left(\cdot \mid {\mathbf{c}}_l,r_l,{\mathbf{\mu }}_k^l,{\tau }_k^l\right),$$

$${\mathbf{x}}_i\sim \sum _{l=1}^{\infty }{\lambda }_l\sum _{k=1}^M{\pi }_k^l{\mathrm{FG}}_{\sigma }\left(\cdot \mid {\mathbf{c}}_l,r_l,{\mathbf{\mu }}_k^l,{\tau }_k^l\right).$$

$$p\left(\mathbf{\mu },\tau \right)\propto {\left\{\frac{{\tau }^{d/2-1}}{I_{d/2-1}(\tau )}\right\}}^a\exp \left(b\tau {\mathbf{\mu }}_0^{\prime }\mathbf{\mu }\right),$$

$$P(s_{i}=l)\propto\begin{cases}n_{-i}^{l}\sum_{k=1}^{M}\pi_{k}^{l}\mathrm{FG}_{\sigma}\left({\bf x}_{i}\mid\boldsymbol{\rho}_{k}^{l}\right)&\mbox{for $l=1,\ldots,L$},\\ \alpha J^{-1}\sum_{k=1}^{M}\pi_{k}^{l}\mathrm{FG}_{\sigma}\left({\bf x}_{i}\mid\boldsymbol{\rho}_{k}^{l}\right)&\mbox{for $l=L+1,\ldots,L+J$},\end{cases}$$

$${\mathbf{\mu }}_{c\mid \cdot }^l={\sigma }_{c\mid \cdot }\times \left\{\frac{1}{{\sigma }^2}\sum _{i=1}^{n}I\left(s_i=l\right)\left({\mathbf{x}}_i-r_l{\mathbf{y}}_i\right)\right\}\text{mbox}and{\sigma }_{c\mid \cdot }^{l2}={\left(\frac{n^l}{{\sigma }^2}+\frac{1}{{\sigma }_c^2}\right)}^{-1}.$$

$${\mu }_{r\mid \cdot }^l={\sigma }_{r|\cdot }\times \left\{\frac{{\mu }_r}{{\sigma }_r^2}+\frac{1}{{\sigma }^2}\sum _{i=1}^{n}I\left(s_i=l\right){\mathbf{y}}_i^{\prime }\left({\mathbf{x}}_i-{\mathbf{c}}_l\right)\right\}\text{mbox}and{\sigma }_{r\mid \cdot }^{l2}={\left(\frac{n^l}{{\sigma }^2}+\frac{1}{{\sigma }_r^2}\right)}^{-1}.$$

$$P(k_i=k\mid s_i=l,{\mathbf{y}}_i,\mathbf{\mu },\tau )=\frac{{\pi }_k^l{\mathrm{FG}}_{\sigma }\left({\mathbf{x}}_i\mid {\mathbf{c}}_l,r_l,{\mathbf{\mu }}_k^l,{\tau }_k^l\right)}{{\sum }_{m=1}^M{\pi }_m^l{\mathrm{FG}}_{\sigma }\left({\mathbf{x}}_i\mid {\mathbf{c}}_l,r_l,{\mathbf{\mu }}_m^l,{\tau }_m^l\right)}.$$

$${\mathbf{\mu }}_{y\mid \cdot }={\tau }_{y\mid \cdot }^{-1}\left\{\frac{r_{s_i}}{{\sigma }^2}\left({\mathbf{x}}_i-{\mathbf{c}}_{s_i}\right)+{\tau }_{k_i}^{s_i}{\mathbf{\mu }}_{k_i}^{s_i}\right\}\text{mbox}and{\tau }_{y\mid \cdot }=\Vert \frac{r_{s_i}}{{\sigma }^2}\left({\mathbf{x}}_i-{\mathbf{c}}_{s_i}\right)+{\tau }_{k_i}^{s_i}{\mathbf{\mu }}_{k_i}^{s_i}\Vert .$$

$${\mathbf{\mu }}_{0\mid \cdot }={\tau }_{\mu \mid \cdot }^{-1}{\mathbf{w}}_k^l\text{mbox}and{\tau }_{\mu \mid \cdot }=\Vert {\mathbf{w}}_k^l\Vert ,\text{mbox}with{\mathbf{w}}_k^l={\tau }_k^l\left\{b{\mathbf{\mu }}_0+\sum _{i=1}^{n}I\left(s_i=l,k_i=k\right){\mathbf{y}}_i\right\}.$$

$$p_k^l\left(\tau \mid \cdot \right)\propto {\left\{\frac{{\tau }^{d/2-1}}{I_{d/2-1}(\tau )}\right\}}^{a+n_k^l}\exp \left[\tau {\left\{b{\mathbf{\mu }}_0+\sum _{i}I\left(s_i=l,k_i=k\right){\mathbf{y}}_i\right\}}^{\prime }{\mathbf{\mu }}_k^l\right].$$

$$f_{P,\sigma }(\mathbf{x}):={\int }_{\mathbf{z}}{\varphi }_{{\sigma }^2}(\mathbf{x}-\mathbf{z})\left\{\int \frac{1}{r^d}\sum _{j=1}^M{\pi }_j{f}_{\mathrm{vMF}}\left(\frac{\mathbf{z}-\mathbf{c}}{r}|\mathbf{c},r,{\mathbf{\mu }}_j,{\tau }_j\right)dP\right\}d\mathbf{z},$$

$$\frac{1}{r^d}{\int }_{{\mathbb{R}}^d}{\varphi }_{{\sigma }^2}(\mathbf{x}-\mathbf{z})f_{\mathrm{vMF}}\left(\frac{\mathbf{z}-\mathbf{c}}{r}|\mathbf{c},r,\mathbf{\mu },\tau \right)d\mathbf{z}$$

$$=\frac{1}{r^d}{C}_d(\tau ){\left(2\pi {\sigma }^2\right)}^{-d/2}{\int }_{\Vert \mathbf{z}-\mathbf{c}\Vert =r}\exp \left\{-\frac{1}{2{\sigma }^2}\Vert \mathbf{x}-\mathbf{z}{\Vert }^2+\tau {\left(\frac{\mathbf{z}-\mathbf{c}}{r}\right)}^{\prime }\mathbf{\mu }\right\}d\mathbf{z}$$

$$=\frac{1}{r^d}[C_d(\tau ){\left(2\pi {\sigma }^2\right)}^{-d/2}{\int }_{\Vert \mathbf{z}-\mathbf{c}\Vert =r}\exp \{-\frac{1}{2{\sigma }^2}\left(\Vert \mathbf{x}-\mathbf{c}{\Vert }^2+r^2\right)$$

$$+{\left(\frac{\mathbf{z}-\mathbf{c}}{r}\right)}^{\prime }\left(\tau \mathbf{\mu }+\frac{r(\mathbf{x}-\mathbf{c})}{2{\sigma }^2}\right)\}d\mathbf{z}]$$

$$=C_d(\tau ){\left(2\pi {\sigma }^2\right)}^{-d/2}{\int }_{\Vert \mathbf{w}\Vert =1}\exp \left\{-\frac{1}{{\sigma }^2}\left(\Vert \mathbf{x}-\mathbf{c}{\Vert }^2+r^2\right)+{\mathbf{w}}^{\prime }\left(\tau \mathbf{\mu }+\frac{r(\mathbf{x}-\mathbf{c})}{2{\sigma }^2}\right)\right\}d\mathbf{w}$$

$$=\frac{C_d(\tau ){\left(2\pi {\sigma }^2\right)}^{-(d/2)}}{C_d\left(\Vert \tau \mu +r(\mathbf{x}-\mathbf{c})/{\sigma }^2\Vert \right)}\exp \left\{-\frac{1}{2{\sigma }^2}\left(\Vert \mathbf{x}-\mathbf{c}{\Vert }^2+r^2\right)\right\}.$$

$${\int }_{{\mathbb{R}}^d}{\mathrm{FG}}_{\sigma }(\mathbf{x}|\mathbf{\rho })d\mathbf{x}={\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\frac{1}{r^d}{\varphi }_{\sigma }(\mathbf{x}-\mathbf{z})f_{\mathrm{vMF}}\left(\frac{\mathbf{z}-\mathbf{c}}{r}|\mathbf{\mu },\tau \right)d\mathbf{z}d\mathbf{x}.$$

$${\int }_{{\mathbb{R}}^d}\frac{1}{r^d}{f}_{\mathrm{vMF}}\left(\frac{\mathbf{z}-\mathbf{c}}{r}|\mathbf{\mu },\tau \right)\left\{{\int }_{{\mathbb{R}}^d}{\varphi }_{\sigma }(\mathbf{x}-\mathbf{z})d\mathbf{x}\right\}d\mathbf{z}={\int }_{{\mathbb{R}}^d}\frac{1}{r^d}{f}_{\mathrm{vMF}}\left(\frac{\mathbf{z}-\mathbf{c}}{r}|\mathbf{\mu },\tau \right)d\mathbf{z}=1.$$

$$\frac{C_d(\tau )}{C_d\left(\Vert \tau \mathbf{\mu }+\mathbf{w}\Vert \right)}=\frac{\Vert \tau \mathbf{\mu }+\mathbf{w}{\Vert }^{-(d/2-1)}{I}_{d/2-1}(\tau \mathbf{\mu }+\mathbf{w})}{(\tau {)}^{-(d/2-1)}{I}_{d/2-1}(\tau )}.$$

$$\exp \left\{-\left|\Vert \tau \mathbf{\mu }+\mathbf{w}\Vert -\tau \right|\right\}\le \frac{(\tau \mathbf{\mu }+\mathbf{w}{)}^{-(d/2-1)}{I}_{d/2-1}(\tau \mathbf{\mu }+\mathbf{w})}{(\tau {)}^{-(d/2-1)}{I}_{d/2-1}(\tau )}\le \exp \left(\left|\Vert \tau \mathbf{\mu }+\mathbf{w}\Vert -\tau \right|\right).$$

$$\left|\tau -\Vert \mathbf{w}\Vert \right|-\tau$$

$$-\Vert \mathbf{w}\Vert$$

$$\Vert \mathbf{w}\Vert -2\tau$$

$$-\Vert \mathbf{w}\Vert$$

$$\exp \left(-2\tau -\Vert w\Vert \right)\le \exp \left(-\Vert \tau \mathbf{\mu }+\mathbf{w}\Vert -\tau \right)$$

$$f_{\epsilon }(\mathbf{c})=\{\begin{array}{ll}{t}_m{f}_0(\mathbf{c})& \Vert \mathbf{c}\Vert <m_{\epsilon },\\ 0& \text{mbox}otherwise,\end{array}$$

$$f_{\epsilon }(\mathbf{x})={\int }_{\Vert \mathbf{c}\Vert <m_{\epsilon }}\sum _{j=1}^M{\pi }_{0,j}\frac{t_m{C}_d({\tau }_j){\left(2\pi {\sigma }_{\epsilon }^2\right)}^{-(d/2)}}{C_d\left(\Vert {\tau }_j{\mathbf{\mu }}_j+r_{\epsilon }(\mathbf{x}-\mathbf{c})/{\sigma }_{\epsilon }^2\Vert \right)}\exp \left\{-\frac{1}{2}\left(\frac{\Vert \mathbf{x}-\mathbf{c}{\Vert }^2}{{\sigma }_{\epsilon }^2}+1\right)\right\}{f}_0(\mathbf{c})d\mathbf{c}.$$

Full text

Estimating densities with nonlinear support using Fisher-Gaussian kernels

Minerva Mukhopadhyay1, Didong Li211footnotemark: 1, and David B Dunson2,3

Department of Mathematics and Statistics, Indian Institute of Technology, Kanpur1,

Department of Mathematics, Duke University2,

Department of Statistical Science, Duke University3 The authors contributed equally to this paper.

Current tools for multivariate density estimation struggle when the density is concentrated near a nonlinear subspace or manifold. Most approaches require choice of a kernel, with the multivariate Gaussian by far the most commonly used. Although heavy-tailed and skewed extensions have been proposed, such kernels cannot capture curvature in the support of the data. This leads to poor performance unless the sample size is very large relative to the dimension of the data. This article proposes a novel generalization of the Gaussian distribution, which includes an additional curvature parameter. We refer to the proposed class as Fisher-Gaussian (FG) kernels, since they arise by sampling from a von Mises-Fisher density on the sphere and adding Gaussian noise. The FG density has an analytic form, and is amenable to straightforward implementation within Bayesian mixture models using Markov chain Monte Carlo. We provide theory on large support, and illustrate gains relative to competitors in simulated and real data applications.

Key Words: Bayesian, kernel density estimation, manifold learning, Markov chain Monte Carlo, mixture model, spherical data, von Mises-Fisher.

1 Introduction

Density estimation is one of the canonical problems in statistics and machine learning. Even when the focus is not directly on the density of the data, density estimation often is a key component of the analysis, arising in clustering, classification, dimensionality reduction, and robust modeling applications among many others. Our focus is on kernel-based approaches to multivariate density estimation, with an emphasis on mixture models. In the simplest setting, one has observed data ${\bf x}_{i}=\left(x_{i,1},\ldots,x_{i,d}\right)\textquoteright,\ i=1,\ldots,n$, which can be assumed to be independent and identically distributed (i.i.d.) draws from an unknown density $f$ that can be expressed as:

[TABLE]

where $\boldsymbol{\pi}=\left(\pi_{1},\ldots,\pi_{M}\right)\textquoteright$ is a vector of probability weights summing to one, $M$ is the number of mixture components, and $\mathcal{K}\left({\bf x};\boldsymbol{\theta}\right)$ is a kernel having parameters $\boldsymbol{\theta}$.

Model (1) is extremely widely used, forming the foundation of model-based clustering, model-based approaches to density estimation and associated non/semiparametric modeling, and Bayesian nonparametrics. There is an immense literature focused on different approaches for choosing the weights $\boldsymbol{\pi}$ and number of components, see Frühwirth-Schnatter, (2006) and Ghosal and van der Vaart, (2017). It is also well known that the choice of kernel $\mathcal{K}({\bf x};\boldsymbol{\theta})$ is very important, not only in applications in which one wants to use the model for clustering and cluster-based inferences (see, e.g., Miller and Dunson, (2018)) but also for accurate density estimation. However, even with this rich literature, the class of kernels that can be used practically in inferences under (1) remains quite restrictive from the perspective of characterizing curvature in the data.

We clarify the problem using two toy examples - the noised spiral and Olympic rings - shown in Figure 1. Both are examples of $d=2$ dimensional data concentrated near a one dimensional curve(s).

These are common toy examples from the literature on nonlinear dimensionality reduction and manifold learning (Lee and Verleysen, (2007)). We use them as motivation, since even in these simple cases problems arise for usual Gaussian kernels due to their intrinsic linear structure. In general for real world datasets, we will have limited prior knowledge about the support of the data, but in general no reason to suspect that it should be close to linear even locally. Indeed, in application areas ranging from imaging to signal processing to astrophysics, there is often strong reason to suggest substantial curvature in the data.

Figure 2 shows results from estimating the density of the data shown in Figure 1 using frequentist kernel density estimation (KDE, see, e.g., Silverman, (1986)) implemented in the ‘ks’ R-package, a Bayesian nonparametric approach relying on Dirichlet process mixtures (DPM, see, e.g., MacEachern and Müller, (1998)) of Gaussian kernels implemented in the ‘dirichletprocess’ R-package, and our proposed approach, which uses (1) within a Bayesian inference approach incorporating a novel class of kernels that accommodate curvature. For easy visualization, we plot new samples from the estimated density in each case. It is clear that the KDE and DPM of Gaussian methods fail to accurately estimate the density. This is due in large part to the fact that the multivariate Gaussian kernel has elliptical contours about a line or plane, leading to challenges in accurately locally approximating the densities of the data shown in Figure 1. The Gaussian kernel provides an accurate local approximation only in a very small region, and hence one needs many kernels and a correspondingly large sample size $n$ for good performance.

A natural question is whether this problem can be solved with the various alternative kernels that have been proposed in the literature. In the multivariate case, the most popular extensions have been to inflate the tails using multivariate-t kernels (McLachlan and Peel, (1998), Lee and McLachlan, (2016)) and/or to generalize the Gaussian distribution to allow skewness (Azzalini and Dalla Valle, (1996), Arellano-Valle and Azzalini, (2009)). Motivated by the problem of more robust clustering, there are also methods available for nonparametrically estimating the kernels subject to unimodality (Rodríguez and Walker, (2014)) and log concavity (Hu et al., (2016)) restrictions. Also to improve robustness of clustering, one can use a mixture of mixtures model employing multiple kernels having similar location parameters to characterize the data within a cluster (Malsiner-Walli et al., (2017)). None of these methods address the fundamental problem we focus on, which is how to better characterize curvature in the support.

An alternative strategy would be to abandon kernel-based methods entirely in attempting to characterize the density of $f$ more accurately. In this regard, one could potentially use flexible non-linear latent structure models developed in the machine learning literature, ranging from Gaussian process latent variable models (GP-LVMs) (Li and Chen, (2016), Doersch, (2016)) to variational auto-encoders (VAEs) leveraging deep neural networks (Lawrence, (2005)). However, such methods are complex black boxes that fail to produce an analytic expression for the density of the data and sacrifice the simple interpretation of (1) in terms of producing data clusters that each have a known distributional form.

With this motivation, we propose to generalize the Gaussian kernel to accommodate curvature in a simple and analytically tractable manner. In particular, we develop a new class of Fisher-Gaussian (FG) kernels by generating observations from a von Mises-Fisher density on a sphere and adding Gaussian noise. The resulting FG distribution has an analytic form, and includes a curvature parameter, in addition to location and scale parameters. Using this kernel, we propose FG mixture models and implement these models within a Bayesian framework using Markov chain Monte Carlo (MCMC).

Section 2 derives the FG class of kernels and studies basic distributional properties and representations. Section 3 includes these kernels within a Bayesian mixture model, and develops a straightforward MCMC algorithm for posterior computation. Section 4 studies asymptotic properties. Section 5 contains a simulation study. Section 6 applies the methods to several datasets, and Section 7 contains a discussion. Proofs are included in the Appendix.

2 The Fisher-Gaussian kernel

Suppose the observed data, ${\bf x}_{1},\ldots,{\bf x}_{n}$, are supported around a sphere with center ${\bf c}$ and radius $r$. Due to noise we do not observe points exactly on the sphere. Therefore we model the data as

[TABLE]

for some unknown noise variance $\sigma^{2}>0$. Here ${\bf c}+r{\bf y}_{i}$ is the part of ${\bf x}_{i}$ supported exactly on the sphere with center ${\bf c}$ and radius $r$, and $\boldsymbol{\epsilon}_{i}$ is a random noise. The random variable ${\bf y}_{i}=\left({\bf x}_{i}-{\bf c}-\boldsymbol{\epsilon}_{i}\right)/r$ has support on the unit sphere, and is assumed to follow a von-Mises Fisher density as ${\bf y}_{i}\sim f_{\mathrm{vMF}}\left(\cdot\mid\boldsymbol{\mu},\tau\right)$, where $\boldsymbol{\mu}$ is the mean direction, and $\tau$ is a concentration parameter. Higher values of $\tau$ indicate higher concentration of the distribution around $\boldsymbol{\mu}$. Thus the density of ${\bf x}$ is

[TABLE]

where $\phi_{\sigma^{2}}(\cdot)$ is the multivariate normal density with mean ${\bf 0}$ and variance $\sigma^{2}I$. We refer to (3) as the Fisher-Gaussian (FG) density.

As the definition suggests, the FG density is essentially a Gaussian density concentrated near the circumference of a sphere. The kernel depends on five parameters: ${\bf c}$ and $r$: the center and radius of the sphere along which the kernel is distributed; $\boldsymbol{\mu}$ and $\tau$: the location on the sphere around which the density has highest concentration and the precision along the sphere; and $\sigma^{2}$ controls how far from the sphere data points concentrate. Figure 3 (a)–(b) demonstrate the effect of the different parameters of the kernel. Throughout this paper, we denote the FG kernel as $\mathrm{FG}_{\sigma}\left(\cdot\mid\boldsymbol{\rho}\right)$ where $\boldsymbol{\rho}=\left({\bf c},r,\boldsymbol{\mu},\tau\right)$.

The FG-kernel has a simple analytic form which makes it theoretically tractable, and easily applicable. By integrating out ${\bf z}$ in (3), the FG kernel can be expressed as follows

[TABLE]

where $C_{d}(\tau)=\tau^{d/2-1}(2\pi)^{-d/2}I_{d/2-1}^{-1}(\tau)$ and $I_{\nu}$ denotes the modified Bessel function of the first kind of order $\nu$. The following lemma proves the above claim.

Lemma 1**.**

1. a.

The density in (3) can be written as in (4). 2. b.

Consequently, for any given $\boldsymbol{\rho}$ and $\sigma^{2}$, $\int_{\mathbb{R}^{d}}\mathrm{FG}_{\sigma}\left({\bf x}\mid\boldsymbol{\rho}\right)d{\bf x}=1$.

The FG kernel provides a simple and flexible generalization of the spherical Gaussian density, commonly used in multivariate density estimation, to allow “curved” densities, as illustrated in Figure 3. The curvature can be modified arbitrarily by changing the radius $r$, with smaller $r$ providing larger curvature and larger $r$ providing densities close to Gaussian. Potentially, one can define a broader class of FG distributions by replacing either or both of the scale parameters $\tau$ and $\sigma^{2}$ with positive semidefinite matrices, obtaining a flexible class of curved multivariate distributions. However, we leave this to future work.

3 Mixtures of Fisher-Gaussian Kernels

Although the FG distribution is useful in other settings as a simple but flexible generalization of the spherical Gaussian density, our focus is on using FG kernels in mixture models. Real data often involve nonlinear relationships among variables, and in such cases mixtures of FG kernels are expected to provide more parsimonious density approximations requiring many fewer mixture components to obtain an accurate approximation.

To make the model more parsimonious, we propose a nested modification to (1). We start by letting ${\bf x}_{i}={\bf z}_{i}+\boldsymbol{\epsilon}_{i}$, as in model (2), with ${\bf z}_{i}$ falling exactly on a sphere. The distribution of ${\bf z}_{i}$ depends on parameters controlling inter (${\bf c}$ and $r$) and intra-sphere ($\boldsymbol{\mu}$ and $\tau$) characteristics. The parameters $({\bf c},r)$ control the location and size of the sphere on which ${\bf z}_{i}$ lies, whereas $(\boldsymbol{\mu},\tau)$ controls the modal point and precision of the distribution of ${\bf z}_{i}$ on the sphere. Within a sphere, the distribution can be unimodal or multimodal. Hence, it is natural to consider a model in which one has a collection of spheres, having varying centers and radii, along with a collection of kernels on each of these spheres. Thus, we can represent the distribution of ${\bf z}_{i}$ as

[TABLE]

where $s_{i}$ indexes the sphere that ${\bf z}_{i}$ lies on, and $k_{i}$ indexes the kernel within the sphere from which ${\bf z}_{i}$ is drawn. The number of spheres, as well as vMF kernels is not known. To keep the representation simple, yet flexible, we let $s_{i}\in\{1,\ldots,\infty\}$ and $k_{i}\in\{1,\ldots,M\}$. Under the proposed specification, we have $P\left(s_{i}=l\right)=\lambda_{l}$ for $l=1,\ldots,\infty$, and $P\left(k_{i}=k\mid s_{i}=l\right)=\pi_{k}^{l}$ for $k=1,\ldots,M$.

Conditionally on $\left(s_{i},k_{i}\right)$ but marginalizing out ${\bf z}_{i}$, we have

[TABLE]

following the FG-kernel described in (4). Further marginalizing out $(s_{i},k_{i})$ leads to the simple two-layer FG kernel mixture

[TABLE]

The proposed representation of ${\bf x}_{i}$ provides a very flexible setup accommodating mixtures of multimodal spherical densities, and is suitable to approximate densities concentrated near an arbitrary curved support. A representative density of $x$, characterized as above, is shown in Figure 3 (c).

A Bayesian model is completed with priors for the two levels of mixture weights and the kernel parameters, as described in the next subsection.

3.1 Prior Specifications

For the mixture distribution over spheres, we use a Dirichlet process prior (DP, Ferguson, (1973), Antoniak, (1974)) with concentration parameter $\alpha$ and base measure $G_{0}=N\left({\bf 0},\sigma^{2}_{c}I\right)\times N_{+}\left(\mu_{r},\sigma_{r}^{2}\right)$, where $\sigma_{c}^{2},\mu_{r},\sigma_{r}^{2}$ are positive constants and $N_{+}$ denotes a Gaussian density truncated to $(0,\infty).$ The DP prior induces a stick-breaking process (Sethuraman,, 1994) on the probability weights $\boldsymbol{\lambda}=\left(\lambda_{1},\lambda_{2},\ldots\right)$. For the $l^{th}$ sphere, the von-Mises kernel weights $\boldsymbol{\pi}^{l}=\left(\pi^{l}_{1},\ldots,\pi^{l}_{M}\right)$ follow a Dirichlet prior with parameters $a_{0},\ldots,a_{0}$ for some fixed constant $a_{0}$. For each sphere and each kernel, the parameters of the von-Mises Fisher density, $\left(\boldsymbol{\mu}^{l}_{m},\tau^{l}_{m}\right)$, are assigned the conjugate prior proposed in Nuñez Antonio and Gutiérrez-Peña, (2005),

[TABLE]

with hyper-parameters $\boldsymbol{\mu}_{0}$, $a$ and $b$ for some $a>b>0$. Given $\tau$, the marginal density of $\boldsymbol{\mu}$ follows a von-Mises Fisher density. However, the marginal density of $\tau$ does not have a standard form. The variance of measurement error $\sigma^{2}$ is assigned an inverse-gamma prior with parameters $a_{\sigma},b_{\sigma}>0$, where $a_{\sigma}$ and $b_{\sigma}$ are constants. The choices of the hyperparameters are discussed in Section 5.

3.2 Posterior Computation

Throughout this section we use the following conventions: (i) The sphere and kernel indicators of the $i^{th}$ observation are denoted by $s_{i}$ and $k_{i}$, respectively, and the sphere and kernel indices are denoted by $l$ and $k$, $l=1,2,\ldots$, $k=1,\ldots,M$; (ii) The number of observations in the $l^{th}$ sphere is denoted by $n^{l}$, and that in the $k^{th}$ kernel of the $l^{th}$ sphere is indicated by $n^{l}_{k}$; (iii) We let $\boldsymbol{\rho}_{k}^{l}=\left({\bf c}_{l},r_{l},\boldsymbol{\mu}_{k}^{l},\tau_{k}^{l}\right)$ denote the parameters of the $k^{th}$ kernel on the $l^{th}$ sphere.

We develop a Metropolis within Gibbs sampler using the characterization of the density in (3). The sampler is straightforward to implement, involving the following simple updating steps from standard distributions and no need for algorithm tuning.

Step 1

Updates on sphere labels. We follow Neal, (2000, Algorithm 8) to update the sphere labels, $s_{1},\ldots,s_{n}$. Excluding the $i^{th}$ observation, suppose there are $L$ spheres represented in the data having parameters $\boldsymbol{\rho}_{k}^{l}$, $\pi_{k}^{l}$ and $\sigma^{2}$ for $l=1,\ldots,L$ and $k=1,\ldots,M$. Fix a number $J$, and create $J$ additional empty sphere labels, $L+1,\ldots,L+J$. If the $i^{th}$ observation was the only data point allocated to sphere $L^{\prime}$ in the last iteration, then attach the corresponding $\boldsymbol{\rho}_{k}^{L^{\prime}}$ and $\pi_{k}^{L^{\prime}}$ to the first empty sphere label, and associate the remaining empty sphere labels with values of $\boldsymbol{\rho}_{k}^{l}$ and $\pi_{k}^{l}$ generated randomly from the base distribution. We draw a new value of $s_{i}$ with

[TABLE]

where $n_{-i}^{l}=\sum_{j\neq i}1\left(s_{j}=l\right)$ and $\mathrm{FG}_{\sigma}\left({\bf x}_{i}\mid\boldsymbol{\rho}\right)$ is as given in (4).

Step 2

Conditional posteriors of centers and radii. The conditional posterior of the center of the $l^{th}$ sphere ${\bf c}_{l}$ follows a Gaussian density with location parameter $\boldsymbol{\mu}_{c\mid\cdot}^{l}$ and scale parameter $\sigma^{l~{}2}_{c\mid\cdot}I$, where

[TABLE]

Similarly, the conditional posterior of the radius of the $l^{th}$ sphere $r_{l}$ follows a positive Gaussian density with location parameter $\mu_{r\mid\cdot}^{l}$ and scale parameter $\sigma^{l~{}2}_{r\mid\cdot}$,

[TABLE]

Step 3

Updates on von-Mises Fisher kernel allocations and weights. We first update the kernel allocation $k_{i}$, for $i=1,\ldots,n$, according to

[TABLE]

Then, the posterior distribution of $\left(\pi_{1}^{l},\ldots,\pi_{M}^{l}\right)$ is Dirichlet with parameters $a_{0}+n^{l}_{1},\ldots,a_{0}+n^{l}_{M}$.

Step 4

Update on spherical coordinates. The latent spherical location, ${\bf y}_{i}$, corresponding to the $i^{th}$ observation ${\bf x}_{i}$, is updated from its conditional posterior density, which is a von-Mises Fisher density with parameters $\boldsymbol{\mu}_{y\mid\cdot}$ and $\tau_{y\mid\cdot}$, where

[TABLE]

Step 5

Update on $\sigma^{2}$. The inverse-gamma prior of $\sigma^{2}$ is conjugate yielding an inverse-gamma posterior with parameters $a_{\sigma\mid\cdot}$ and $b_{\sigma\mid\cdot}$ where

$a_{\sigma\mid\cdot}=a_{\sigma}+\displaystyle\frac{nd}{2}\quad\mbox{and}\quad b_{\sigma\mid\cdot}=b_{\sigma}+\displaystyle\frac{1}{2}\displaystyle\sum_{i=1}^{n}\left({\bf x}_{i}-{\bf c}_{s_{i}}-r_{s_{i}}{\bf y}_{i}\right)^{\prime}\left({\bf x}_{i}-{\bf c}_{s_{i}}-r_{s_{i}}{\bf y}_{i}\right).$

Step 6

Update on parameters of von-Mises Fisher distributions. Finally, for each sphere, $l$, and each kernel, $k$, the vMF parameters are updated based on the observations associated with that sphere-kernel combination. The conditional prior on $\boldsymbol{\mu}^{l}_{k}$ given $\tau^{l}_{k}$ is von-Mises Fisher, and is conjugate yielding a vMF posterior with parameters $\boldsymbol{\mu}_{0\mid\cdot}$ and $\tau_{\mu\mid\cdot}$ where

[TABLE]

The conditional posterior of $\tau^{l}_{k}$ given other parameters is

[TABLE]

We apply a Metropolis-Hastings algorithm with an independence sampler, considering a $\mathrm{Gamma}\left(2,u^{l}_{k}\right)$ proposal distribution. The details on choice of $u^{l}_{k}$ are given in Section 5.

4 Asymptotic Properties

In this section, we show the Kullback-Leibler (KL) support property of the proposed prior, which implies weak posterior consistency (see Schwartz, (1965)).

Let $f_{0}$ be the true density function. Under the mixture of FG-kernel approach, the proposed density function $f$ is described in (6). For simplicity we show the KL property holds under a restrictive version of our model in which the same vMF parameters are re-used across the different spheres. Under this assumption, we have

[TABLE]

where the mixing distribution of $({\bf c},r,\boldsymbol{\pi})$, $P$, follows a DP prior with base measure $N\left(0,\sigma^{2}_{c}I_{d}\right)\times N_{\delta}\left(\mu_{r},\sigma^{2}_{r}\right)\times\mathrm{Dir}\left(p,\ldots,p\right)$. Here $\{\boldsymbol{\mu}_{1},\ldots,\boldsymbol{\mu}_{M}\}$ and $\{\tau_{1},\ldots,\tau_{M}\}$ are prefixed quantities. The hyperparameters $\sigma_{c}^{2}$, $\mu_{r}$, $\delta$ and $\sigma_{r}^{2}$ are assumed to be elicited a priori.

First, we state a lemma which makes the mixture density in (6) more theoretically tractable.

Lemma 2**.**

For any $\boldsymbol{\rho}=({\bf c},r,\boldsymbol{\mu},\tau)$ and $\sigma$, let the FG kernel be as in (4), then the following holds:

1. a.

For $d>2$, $\displaystyle(2\pi\sigma^{2})^{-d/2}\exp\left\{-\displaystyle\frac{1}{2\sigma^{2}}\left(\|{\bf x}-{\bf c}\|+r\right)^{2}\right\}\leq\mathrm{FG}_{\sigma}\left({\bf x}\mid\boldsymbol{\rho}\right)$

$\leq 2(2\pi\sigma^{2})^{-d/2}\exp\left\{-\displaystyle\frac{1}{2\sigma^{2}}\left(\displaystyle\|{\bf x}-{\bf c}\|-r\right)^{2}\right\}.$** 2. b.

For $d=2$, the above inequality holds with left and right hand sides multiplied by $\exp(-2\tau)$ and $\exp(2\tau)$, respectively.

Let the true density be $f_{0}$, and $KL\left(f_{0},f\right)=\int f_{0}\log\left(f_{0}/f\right)$ be the KL divergence between $f_{0}$ and any density $f$. Let the KL neighborhood of $f_{0}$ of size $\varepsilon$ be $\mathcal{K}_{\varepsilon}(f_{0})=\left\{f:KL\left(f_{0},f\right)<\varepsilon\right\}$. We say that $f_{0}$ is in KL support of $\Pi$, written as $f_{0}\in KL\left(\Pi\right)$, if $\Pi\left\{\mathcal{K}_{\varepsilon}\left(f_{0}\right)\right\}>0$, for every $\varepsilon>0$.

Next we consider the assumptions under which our results hold.

1. A.

The true density $f_{0}({\bf x})<M_{0}$ for some constant $M_{0}>0$ for all ${\bf x}\in\mathbb{R}^{d}$. 2. B.

The true density $f_{0}$ satisfies: $|f_{0}({\bf x}+{\bf y})-f_{0}({\bf x})|\leq L({\bf x})\exp\left(s\|{\bf y}\|^{2}\right)\|{\bf y}\|^{r_{0}}$, for any ${\bf x},{\bf y}\in\mathbb{R}^{d}$, some $s,r_{0}\geq 0$, and a non-negative function $L$ in $\mathbb{R}^{d}$ such that $L({\bf x})<M^{\prime}f_{0}({\bf x})$, for some large constant $M^{\prime}$. 3. C.

The $q^{th}$ order moments of ${\bf x}$ are finite, for some $q>2$. 4. D.

There exists a fixed point ${\bf x}_{0}$ and some fixed $\delta>0$ such that $\inf_{\|{\bf x}-{\bf x}_{0}\|<\delta}f_{0}({\bf x})=\phi_{\delta}>0$.

Assumption A is common in the literature (see, for e.g., Wu and Ghosal, (2008)). Assumption B is weaker than the usual Hölder continuity with $\alpha=r_{0}$. For $r_{0}<1$, $f_{0}$ satisfying assumption B belongs to a locally $r_{0}$-Hölder class with envelope $L$ (see Shen et al., (2013)). Assumptions C and D are weak conditions.

Theorem 1**.**

Suppose the true density $f_{0}$ satisfies assumptions A-D, and the prior $\Pi$ is as given in (8). Then $\Pi\left\{\mathcal{K}_{\varepsilon}(f_{0})\right\}>0$ for any $\varepsilon>0$.

Remark 1**.**

In Theorem 1, we consider the mixing distribution of $({\bf c},r,\boldsymbol{\pi})$ is $P$, $P\sim DP$, and $(\boldsymbol{\mu},\tau)$ varies on a fixed set. One can trivially extend Theorem 1 to the case where the mixing distribution of $({\bf c},r,\boldsymbol{\pi},\boldsymbol{\mu}_{1},\cdots,\boldsymbol{\mu}_{M},\tau_{1},\cdots,\tau_{M})$ is $P$, if $\tau$ is restricted to a compact set $[0,\tau_{\max}]$.

5 Simulation Experiments

In this section we carry out a variety of simulation experiments to validate our method. Our focus is on datasets which are concentrated near a manifold. We compare our method with frequentist kernel density estimation (KDE) and Dirichlet process mixtures of Gaussians (DPMG).

Our method: The algorithm is described in Section 3.2. We keep the choices of precision parameter $\alpha$, and number of empty spheres $J$ fixed, and consider two choices of $\alpha/J$, $\{1,2\}$. The number of von-Mises Fisher (vMF) kernels within each sphere is set to $5,10$ or $15$ based on the complexity of the data. The parameter $a_{0}$ in the prior for the vMF kernel weights $\boldsymbol{\pi}$ is set to $1$. The hyperparameters of the inverse-gamma prior of $\sigma^{2}$ are set to $\mathrm{rate}=1$ and $\mathrm{scale}=0.01$ to make the prior weakly informative. The hyperparameters of the centers and radii are set to $\sigma_{c}^{2}=1$, $\mu_{r}=1$ and $\sigma_{r}=5$. The $\boldsymbol{\mu}_{j}$s are given a vMF prior with parameters $\boldsymbol{\mu}_{0}=\left(1/\sqrt{d},\ldots,1/\sqrt{d}\right)$ and $a=1$ with $b=0.1$. A small value of $b$ makes the prior less informative. For the independence sampler applied to update $\tau^{l}_{m}$, we use a gamma proposal distribution with shape $2$ and rate $2/\hat{\tau}^{l}_{m}$. Here $\hat{\tau}^{l}_{m}$ is an approximation to the maximum likelihood (ML) estimate of $\tau^{l}_{m}$, $\hat{\tau}^{l}_{m}=\bar{t}\left(d-\bar{t}\right)/\left(1-\bar{t}^{2}\right)$, where $\bar{t}=\|\sum_{j}{\bf y}^{l}_{j,m}\|/n^{l}_{m}$ (see Sra, (2012) for details on ML estimation of vMF parameters). Here ${\bf y}_{j,m}^{l}$ denotes the latent spherical coordinates of the $j^{th}$ observation belonging to the $l^{th}$ sphere and $m^{th}$ kernel. This proposal provides acceptance rates between $20\%$-$40\%$ in our experiments.

Other methods: For comparison we implement DMPG using the ‘DirichletProcessMvnormal’ function in the ‘dirichletprocess’ R-package with all default settings. Multivariate Gaussian kernels are used, with the base measure of the DP corresponding to a multivariate normal-inverse Wishart distribution and with the concentration parameter $\alpha$ given a gamma prior. Algorithm 8 in Neal, (2000) is used for sampling. For kernel density estimation (KDE) we use the ‘kde’ function with default settings in the ‘ks’ R-package. This corresponds to using multivariate Gaussian kernels with covariance/bandwidth matrix $H$. By default ‘ks’ selects the bandwidth $H$ by direct plug-in methodology, as described in Wand and Jones, (1994) and Chacón and Duong, (2010) for univariate and multivariate cases, respectively. Most of the packages including ‘ks’ do not allow data of dimension greater than 6 (see, for e.g., Deng and Wickham, (2011)). Thus for cases with dimension greater than 6, we use the ‘kdevine’ package. This package estimates marginal densities using KDE, and then uses the VINE copula to obtain the multivariate joint density estimate.

We split the validation into two categories: density estimation and density-based classification. Below we briefly describe the measures used to assess performance.

(i) Histogram estimate. We generate $n$ samples, say ${\bf z}_{1},\ldots,{\bf z}_{n}$, from the predictive distribution, $n$ being the training size; and choose $N$ points, ${\bf x}_{1},\ldots,{\bf x}_{N}$, randomly with replacement from the training data. Fix a neighborhood size $\delta$. Let $W_{j,\delta}^{z}=\sum_{i=1}^{n}I(\|{\bf z}_{i}-{\bf x}_{j}\|<\delta)$ be the number of predictive points and $W_{j,\delta}^{x}=\sum_{i=1}^{n}I(\|{\bf x}_{i}-{\bf x}_{j}\|<\delta)$ be the number of training points in the $\delta$-neighborhood of ${\bf x}_{j}$, $j=1,\ldots,N$. Then the histogram estimate is the mean absolute difference $\mathrm{MAD}_{\delta}=\sum_{j=1}^{N}|W_{j,\delta}^{z}-W_{j,\delta}^{x}|/N$.

(ii) Likelihood estimate: Splitting the data into training and test sets, we estimate the density $f$ using the training set. We find the predictive likelihood of each test point ${\bf x}$, $\hat{f}({\bf x})$. Out-of-sample goodness of fit is summarized via the estimated log-likelihood, $\sum_{j}\log\hat{f}({\bf x}_{j})$, or boxplots of $\hat{f}({\bf x}_{j})$.

(iii) Classification Accuracy. Using the training data, we separately estimate the density of the features within each class. Letting $\hat{f}_{j}({\bf x})$ denote the density within class $j$, we assign the test data ${\bf x}^{*}$ to the class having the highest $\hat{f}_{j}\left({\bf x}^{*}\right)$. We then report the classification accuracy as the proportion of correct classifications in the test data.

Below we consider three density estimation examples along with a classification example.

Euler Spiral

An Euler spiral is a curve whose curvature changes linearly with its curve length. We consider a sample of size $500$ from the Euler spiral, and add a Gaussian noise of variance $0.001^{2}I$. The mean absolute differences MADδ for $N=2000$ and different values of $\delta$ are given in Table 2. Box-plots of estimated likelihoods of $300$ test points are given in Figure 4. A Scatter plot of an Euler spiral data and that of $500$ predictive points from each method are shown in Figure 5.

Olympic Rings

We generate $1500$ datapoints from 5 circles forming Olympic rings. The five circles have centers at $(0,0)$, $(2.125,0)$, $(4.25,0)$, $(1.125,-1)$ and $(3.25,-1)$, respectively, and radii $1$. We generate $i\times 100$ points from the $i^{th}$ circle, and add a Gaussian noise component with sd $0.01$ for each component. The values of MADδ for $N=2000$ and different values of $\delta$ are given in Table 2, and box-plots of estimated likelihoods are shown in Figure 5.

Torus

A torus (or torus of revolution) is a surface of revolution generated by revolving a circle in three-dimensional space about an axis coplanar with the circle. The axis of revolution does not touch the circle. We generate a dataset of $n=1500$ points on a torus using the ‘torus’ function in the ‘geozoo’ R-package. The radius of the larger circle is set to $3$, that of the smaller circle is set to $1$, and a Gaussian noise with variance $0.1$ is added. The values of MADδ for $N=1000$ and different values of $\delta$ are given in Table 4, figures of the Torus dataset and $1500$ predictive points are shown in Figure 4, and box-plots of estimated likelihoods are in Figure 5.

Two-Spiral

The last example is that of two spirals with varying curvature. The spirals intersect only at the center at the point of maximum curvature. As the curves move outward, the curvature of the spirals decreases simultaneously. The classification accuracies for different training test splits are given in Table 4, and the scatter plot of the actual dataset and $750$ predictive points from each spiral are shown in Figure 6.

The simulation results show almost uniformly better performance for the FG-mixture approach, with gains very dramatic in some cases.

6 Applications

In this section, we consider three datasets, Galaxy, Gesture Phase and Balance scale, for classification, and two datasets, User knowledge modeling and Balance scale, for density estimation. The results are given below.

Galaxy dataset

is available at https://data.galaxyzoo.org/. In the original Galaxy Zoo project (see Lintott et al., (2008) for details), volunteers classified images of Sloan Digital Sky Survey galaxies as belonging to one of six categories. We consider the first dataset (Table 2), where the galaxies are classified as spiral, elliptical or merger (i.e., uncertain). This dataset has information on 10,000 galaxies, among which 2888 are spiral and 931 are elliptical. Among these samples we consider $100$ spiral and $100$ elliptical galaxies at random as the test set, and vary the size of the training set. We consider 5 predictors $\mathrm{P_{EL}}$, $\mathrm{P_{CW}}$, $\mathrm{P_{EDGE}}$, $\mathrm{P_{DK}}$, $\mathrm{P_{MG}}$ for classification. The classification accuracy is given in Figure 7, showing that it is improved significantly using FG mixtures.

Gesture Phase dataset

is available at the UCI machine learning repository. It is composed of features extracted from 7 videos with people gesticulating, aimed at studying gesture phases. We use the processed table a1_va3, which has $32$ attributes, and $5$ labels, viz., D, S, H, P, R. Among these $5$ labels we discard H as it has fewer instances, and compare each pair of labels D, S, P and R. The results are shown in Table 5.

Balance scale dataset

is generated to model psychological experimental results by Klahr and Siegler, (1978), and is obtained from the UCI repository. Each example is classified as having the balance scale tip to the right, left, or balanced. The predictors are left-weight, left-balance, right-weight and right-balance. There are 288 instances of each of left and right tipped balance scales, and we consider classification between left and right tipped balance scales. The classification accuracies are given in Figure 7.

We further use Balance scale data for density estimation. Towards that we consider training samples of various sizes from the unlabeled population, and a test data of size 100. The sum of estimated log-likelihoods of the test data is provided in Table 7.

User knowledge modeling dataset

is about the students’ knowledge status on the subject of Electrical DC Machines, and the data are obtained from the UCI repository. There are 5 attributes measuring the study time, repetition of different aspects and knowledge level, and 258 instances. The sum of estimated log likelihoods of test samples for varying training-test splits are given in Table 7

From these applications it is evident that FG mixtures clearly improve performance in both density estimation and classification, particularly in smaller samples. Also some of the improvements are quite dramatic.

7 Discussion

The main contribution of this article has been to introduce a simple generalization of the Gaussian distribution to allow the density to be concentrated near a curved support. The resulting Fisher-Gaussian distribution has a simple analytic form, and we have demonstrated that it can lead to dramatically better performance in density estimation when used in place of usual Gaussian kernels. This is particularly true when the data are concentrated near a non-linear subspace or manifold and when there are non-linear relationships among the variables being modeled, as is very commonly the case in applications. There have been many different multivariate distributions proposed in the literature, motivated largely by limitations of the multivariate Gaussian distribution in terms of symmetry and light tails, but essentially no consideration of the curved support problem. Multivariate distributions that can characterize data concentrated near a curved support tend to have a very complex form that is not analytically tractable, and hence there are considerable computational and interpretability hurdles in their implementation.

There are a number of interesting directions for future research. Natural extensions include generalizing the Fisher-Gaussian distribution to accommodate a more flexible covariance structure. One possibility is to include an arbitrary covariance in the Gaussian residual density instead of focusing on the spherical covariance case. Another is to use a more flexible generalization of the von Mises-Fisher density on the sphere; several such generalizations are available in the literature (Bingham, (1974),Kent, (1982)). We have made initial attempts along these directions and have thus far been unable to obtain an analytically tractable form for the density of the data ${\bf x}_{i}$ marginalizing out the coordinates on the sphere ${\bf y}_{i}$. It is possible to conduct computation without such a closed form, but it will be unwieldy.

Another direction is to start with a density on a more flexible non-linear manifold than the sphere before adding Gaussian noise. Possibilities include an ellipse, quadratic surface or hyperbolic space. Unfortunately, it is not so straightforward to define tractable densities on such spaces, and it seems daunting to maintain analytic tractability in doing so. One promising direction is to start with a Gaussian on a tangent plane to a manifold and then apply the exponential map to induce a density on the target manifold. For a number of manifolds, it becomes possible to define the Jacobian and obtain a tractable form for the density on the manifold. For certain cases, it may also be possible to add Gaussian noise and marginalize out the coordinates on the manifold, as we did for the FG distribution. Otherwise, one can always rely on data augmentation algorithms for computation.

As a final fascinating direction, it will be interesting to develop theory directly justifying the use of ‘curved support’ kernels in density estimation. We have focused on showing weak posterior consistency but it would be appealing to show that a faster rate of posterior concentration can be obtained by using FG kernels instead of spherical Gaussian kernels, at least when the true density is concentrated near a non-linear manifold with sufficiently large curvature.

Acknowledgements

This research was partially supported by grant N00014-14-1-0245/N00014-16-1-2147 of the United States Office of Naval Research (ONR) and 5R01ES027498-02 of the United States National Institutes of Health (NIH).

8 Appendix

Proof of Lemma 1

Proof of part (a)..

Observe that

[TABLE]

Proof of part (b). Given $\boldsymbol{\rho}=({\bf c},r,\boldsymbol{\mu},\tau)$ we can write

[TABLE]

Next we interchange the integrals of ${\bf x}$ and ${\bf z}$ as the integrand is non-negative and the integral is finite. After the changing the integral we get:

[TABLE]

∎

Proof of Lemma 2

Proof.

Let $r({\bf x}-{\bf c})/\sigma^{2}={\bf w}$. Note that

[TABLE]

By (Baricz,, 2010, Equation (2.2)), for $d>2$,

[TABLE]

Next, observe that

[TABLE]

as $\|\boldsymbol{\mu}\|=1$ and $|\|{\bf a}\|-\|{\bf b}\||\leq\|{\bf a}+{\bf b}\|\leq\|{\bf a}\|+\|{\bf b}\|$. This reduces to

[TABLE]

Finally, when $\tau<\|{\bf w}\|$ then $\|{\bf w}\|-2\tau\geq-\|{\bf w}\|$, and when $2\tau-\|{\bf w}\|<\|{\bf w}\|$. Thus, $-\|{\bf w}\|\leq\left|\|\tau\boldsymbol{\mu}+{\bf w}\|-\tau\right|\leq\|{\bf w}\|$ for all choices of $\tau$ and ${\bf w}$. The result for $d>2$ immediately follows by inserting the above inequality in the kernel.

For $d=2$, i.e., $(d/2-1)=0$, observe that for $x\geq 0$, $I_{0}(x)=\frac{1}{\pi}\int_{0}^{\pi}\exp\left(x\mathrm{cos}\theta\right)d\theta\implies e^{-x}\leq I_{0}(x)\leq e^{x}$. Therefore,

[TABLE]

This completes the proof for $d=2$. ∎

8.1 Proof of Theorem 1

Proof.

The proof uses ideas of similar proofs from Wu and Ghosal, (2008). The proof is done in three parts shown in the following lemma, which is same as Theorem 1 in Wu and Ghosal, (2008).

Lemma 3**.**

Let $\sigma^{2}\in S$, $\boldsymbol{\theta}\in\Theta$, $P$ be the mixing distribution on $\Theta$, i.e., $f_{P,\sigma^{2}}(\cdot)=\int\mathrm{FG}_{\sigma}(\cdot\mid\Theta)dP(\Theta)$ where $\mathrm{FG}_{\sigma}(\cdot\mid\Theta)$ is as described in (4). Further, let $\mathcal{M}\left(\Theta\right)$ be the space of probability measures on $\Theta$. Suppose there exists $P_{\varepsilon}$, $\sigma^{2}_{\varepsilon}$, $S\in\mathbb{R}^{+}$ and $\mathcal{W}\in\mathcal{M}\left(\Theta\right)$, such that $p(\mathcal{W})>0$ and $p(S)>0$, where $p$ denotes the prior distribution, and the followings hold:

1. I.

$\displaystyle\int_{\mathbb{R}^{d}}f_{0}({\bf x})\log\displaystyle\frac{f_{0}({\bf x})}{f_{P_{\varepsilon},\sigma_{\varepsilon}^{2}}({\bf x})}d{\bf x}<\varepsilon$, 2. II.

$\displaystyle\int_{\mathbb{R}^{d}}f_{0}({\bf x})\log\displaystyle\frac{f_{P_{\varepsilon},\sigma_{\varepsilon}^{2}}({\bf x})}{f_{P_{\varepsilon},\sigma^{2}}({\bf x})}d{\bf x}<\varepsilon$* for every $\sigma^{2}\in S$,* 3. III.

$\displaystyle\int_{\mathbb{R}^{d}}f_{0}\log\displaystyle\frac{f_{P_{\varepsilon},\sigma^{2}}({\bf x})}{f_{P,\sigma^{2}}({\bf x})}d{\bf x}<\varepsilon$* for every $P\in\mathcal{W}$ and $\sigma^{2}\in S$,*

then $f_{0}\in KL\left(\Pi\right)$, where $KL(\Pi)$ is as described in Section 4.

The proof is same as the proof of (Wu and Ghosal,, 2008, Theorem 1), as $\sigma^{2}$ is independent of $P$.

The proof of Theorem 1 relies on showing that the above conditions hold for our proposed prior. Without loss of generality we assume that the fixed point ${\bf x}_{0}$ in assumption C is ${\bf 0}$.

Showing condition (I) holds: We choose $P_{\varepsilon}$ and $\sigma_{\varepsilon}^{2}$ as follows. For any $\varepsilon>0$ we fix an $m_{\varepsilon}$ such that

[TABLE]

where $t_{m}^{-1}=\int_{\|{\bf c}\|<m_{\varepsilon}}f_{0}({\bf c})dc$, and $r=m_{\varepsilon}^{-\eta}$ for some $\eta>0$ with probability 1. The weight $\boldsymbol{\pi}_{\varepsilon}=\boldsymbol{\pi}_{0}$ for any $\boldsymbol{\pi}_{0}\in\Delta_{M}$, $\Delta_{M}$ being the $M$-dimensional simplex, with probability 1. Let $F_{\varepsilon}$ be the probability measure corresponding to $f_{\varepsilon}$, then $P_{\varepsilon}=F_{\varepsilon}\times\delta\left(m_{\varepsilon}^{-\eta}\right)$, where $\delta$ denotes the degenerate distribution. Also let $\sigma_{\varepsilon}=r_{\varepsilon}=m_{\varepsilon}^{-\eta}$, with probability 1. Here $(m_{\varepsilon},\eta)$ is such that $E\left(\|{\bf x}\|^{2(\eta+1)}\right)<\infty$ (see assumption C), and $m_{\varepsilon}^{2\eta}>4s$, where $s$ is as in assumption B.

Therefore, $f_{\varepsilon}$ is as follows:

[TABLE]

Next taking a transformation $({\bf c}-{\bf x})/\sigma_{\varepsilon}={\bf w}$, we get $f_{\varepsilon}({\bf x})$ as follows

[TABLE]

where $I(\cdot)$ is the indicator function. Observe that, as $m_{\varepsilon}\rightarrow\infty$, $t_{m}\rightarrow 1$ and the integrand above converges to $\sum_{j}\pi_{0,j}h_{1}\left({\bf w},({\bf 0},1),(\boldsymbol{\mu}_{j},\tau_{j})\right)f_{0}({\bf x})$.

Observe that $t_{m}\leq t_{1}$ for $m>1$, and by Lemma 2 the above integrand is no bigger than

[TABLE]

for $d>2$ and a suitable constant $M_{1}>0$ due to assumption A. The last term is ${\bf w}$-integrable. Using polar transformation one can show that the above integral is equal to $M_{1}t_{1}2^{-(d-3)/2}\pi^{-1/2}/\Gamma(d/2)$. For $d=2$, the above integrand is the same as (9) with a multiplier $\sum_{j}\pi_{0,j}\exp\left(2\tau_{j}\right)\leq\exp\left(2\tau_{\max}\right)$, and therefore ${\bf w}$-integrable. Here $\tau_{\max}=\max\{\tau_{1},\ldots,\tau_{M}\}$. Therefore by the dominated convergence theorem (DCT), $f_{\varepsilon}({\bf x})\rightarrow f_{0}({\bf x})$, as $m_{\varepsilon}\rightarrow\infty$.

Next we show that $\left|f_{0}({\bf x})\log\left(f_{0}({\bf x})/f_{\varepsilon}({\bf x})\right)\right|$ is bounded by an integrable function so that the DCT can be applied to the integrand in (I). Observe that

[TABLE]

where the inequality in the second line holds by Lemma 2, and that in the last line by assumption B. Without loss of generality assume that $m_{\varepsilon}>M^{*}$ for some $M^{*}$ such that $\exp\{2\tau_{\max}\}t_{m}<t_{0}$ for some positive constant $t_{0}$. Therefore the above integral reduces to

[TABLE]

Observe that both the integrals are finite as long as $s\sigma_{\varepsilon}^{2}<1/3$ (i.e., $s<m^{2\eta}/3)$. To see this, observe that $\sigma_{\varepsilon}<1$, and $s\sigma_{\varepsilon}^{2}\eqqcolon\xi<1/3$. Hence the second integral is no bigger than

[TABLE]

Now it is easy to see that the above integral is finite for $\xi<1/3$. In particular, one can use polar transformation to evaluate the same. Similarly one can also show that the first integral is also finite. Therefore

[TABLE]

as $\log(1+u)<u$ and $E_{0}\left\{L({\bf x})/f_{0}({\bf x})\right\}<\infty$, i.e., $L$ is integrable by assumption B.

Next we find a lower-bound $g({\bf x})$ of $f_{\varepsilon}({\bf x})$. Towards that from Baricz, (2010), $y^{-\nu}I_{\nu}(y)$ is a strictly increasing function over $(0,\infty)$, and therefore $y^{\nu}I^{-1}_{\nu}(y)$ is strictly decreasing over $(0,\infty)$. Therefore $C_{d}(\tau_{0})/C_{d}(\tau)>1$ for any $(\tau,\tau_{0})$ such that $\tau>\tau_{0}$. For $\tau<\tau_{0}$, from Baricz, (2010, Equation (2.2)) we have $C_{d}(\tau_{0})/C_{d}(\tau)>\exp\left(\tau-\tau_{0}\right)>\exp\left(-\tau_{0}\right)$. Therefore,

[TABLE]

Now let $\|{\bf x}\|\geq m_{\varepsilon}/2$. Then the above expression is no less than

[TABLE]

Note that $g(m_{\varepsilon})$ is a decreasing function of $m_{\varepsilon}$ for sufficiently large $m_{\varepsilon}$. When $m_{\varepsilon}$ is large, then $g(m_{\varepsilon})\leq g(2\|{\bf x}\|)$, i.e., the above expression is no less than

[TABLE]

for some suitable constants $K_{1}$ and $K_{2}$.

Let $\|{\bf x}\|<m_{\varepsilon}/2$, $0<\delta<m_{\varepsilon}/3$ fixed, and $\phi_{m}({\bf x})=\inf_{\|{\bf t}-{\bf x}\|<\delta\sigma_{\varepsilon}}f_{0}({\bf t})$, then $\{{\bf c}:\|{\bf c}\|<m_{\varepsilon}\}\cap\{\|{\bf c}-{\bf x}\|\leq\delta\sigma^{2}_{\varepsilon}\}=\{{\bf c}:\|{\bf c}-{\bf x}\|\leq\delta\sigma^{2}_{\varepsilon}\}$ for some small enough $\delta>0$ and

[TABLE]

We bound the last term inside the bracket by

[TABLE]

Let $\sigma_{\varepsilon}$ be sufficiently small (i.e., $m_{\varepsilon}$ be sufficiently large) such that $0<s\sigma_{\varepsilon}^{2}<1/4$, which implies that $1/2<1-2s\sigma_{\varepsilon}^{2}<1$. Then the above term is no bigger than

[TABLE]

Combining the above, noting that $t_{m}\geq 1$ and by assumption B. we get:

[TABLE]

for some suitable constant $0<K^{*}<1$. The last inequality is true of we take $m_{\varepsilon}$ large enough and $\delta$ small enough. Therefore

[TABLE]

for any $R<m_{\varepsilon}/2$. Consequently,

[TABLE]

Recalling that $1/K^{*}\geq 1$, we have

[TABLE]

where $A\coloneqq\left\{{\bf x}:\|{\bf x}\|\geq R,~{}~{}\mbox{and}~{}~{}f_{0}({\bf x})\geq K_{1}\|{\bf x}\|^{d\eta}\exp\left(-K_{2}\|{\bf x}\|^{2(\eta+1)}\right)\right\}$. To see that the second term above is finite, observe that as $f_{0}({\bf x})\geq K_{1}\|{\bf x}\|^{d\eta}\exp\left(-K_{2}\|{\bf x}\|^{2(\eta+1)}\right)\geq 0$ for all ${\bf x}$ in $A$, therefore

[TABLE]

If we show that the RHS of (13) has a lower-bound and the LHS of (13) has an upper-bound then the second term of (12) is also finite. Towards that,

[TABLE]

for some constant $C>0$, $K_{3}=K_{1}/d\eta$ as $E\left(\|{\bf x}\|^{2(\eta+1)}\right)<\infty$, and $\|{\bf x}\|\geq R$ in $A$. Again $\int f_{0}({\bf x})\log(f_{0})d{\bf x}<\log M$ as $f_{0}({\bf x})<M$ for all ${\bf x}\in\mathbb{R}^{d}$. Thus (12) is finite.

As $|a|=\max\{a,-a\}$, from (11) and the above arguments $f_{0}\log(f_{0}/f_{\varepsilon})\leq|f_{0}\log(f_{0}/f_{\varepsilon})|$ and $\int|f_{0}\log(f_{0}/f_{\varepsilon})|$ is finite, and therefore the DCT holds. Thus (I) is satisfied.

Showing condition (II) holds: The proof follows in 4 steps:

1. i.

We will first show that $h_{\sigma}\left\{{\bf x},({\bf c},r)\right\}\rightarrow h_{\sigma_{\varepsilon}}\left\{{\bf x},({\bf c},r)\right\}$ for any given ${\bf x}$ and $({\bf c}^{\prime},r)^{\prime}$, as $\sigma\rightarrow\sigma_{\varepsilon}$. 2. ii.

We use the DCT to show that $f_{P_{\varepsilon},\sigma}\rightarrow f_{P_{\varepsilon},\sigma_{\varepsilon}}$ for $\sigma$ in the neighborhood of $\sigma_{\varepsilon}$, $\sigma\in N(\sigma_{\varepsilon})$, by choosing an appropriate dominating function. 3. iii.

Then we apply (i) to show that $\displaystyle\int_{\mathbb{R}^{d}}f_{0}({\bf x})\log\displaystyle\frac{f_{P_{\varepsilon},\sigma_{\varepsilon}}({\bf x})}{f_{P_{\varepsilon},\sigma}({\bf x})}d{\bf x}\rightarrow 0$ using DCT as $\sigma\rightarrow\sigma_{\varepsilon}$, by finding an appropriate dominating function on $\boldsymbol{\theta}\in D$ where $\{({\bf c},r,\boldsymbol{\pi}_{0}):\|{\bf c}\|<m_{\varepsilon},r=m_{\varepsilon}^{-\eta}\}\subseteq D$. 4. iv.

The proof is completed by taking $S=\{\sigma:|\sigma-\sigma_{\varepsilon}|<\delta\}\cap N(\sigma_{\varepsilon})$, showing that $\pi(S)>0$.

The proof of (i) follows from the fact that $h_{\sigma}\left\{{\bf x},({\bf c},r)\right\}$ is a (pointwise) continuous function of $\sigma$ on $(0,\infty)$ given ${\bf x}$ and $({\bf c},r)$. To see this observe that $\sigma^{-d}$, $\exp\left\{-\left(\|{\bf x}-{\bf c}\|^{2}+r^{2}\right)/\left(2\sigma^{2}\right)\right\}$ and $\|\tau\boldsymbol{\mu}_{0}+r({\bf x}-{\bf c})/\sigma^{2}\|$ are continuous in $(0,\infty)$. Further, for any fixed $\tau\in\{\tau_{1},\ldots,\tau_{M}\}$, $C_{d}^{-1}(\tau)\leq C_{d}^{-1}(\tau_{\max})\approx(2\pi)^{d/2}\sum_{m=0}^{\infty}(\tau_{\max}/2)^{2m}/\left\{m!\Gamma(d/2+m)\right\}$ is finite. These together with the facts that product and convolution of continuous functions are continuous, shows that $h_{\sigma}\left\{{\bf x},({\bf c},r)\right\}\rightarrow h_{\sigma_{\varepsilon}}\left\{{\bf x},({\bf c},r)\right\}$ as $\sigma\rightarrow\sigma_{\varepsilon}$ pointwise.

To show (ii) we first fix a neighborhood of $\sigma_{\varepsilon}$. Recall that $\sigma_{\varepsilon}=m_{\varepsilon}^{-\eta}$. We fix the neighborhood as $N(\sigma_{\varepsilon})=(\underline{\sigma},\bar{\sigma})$, such that $\underline{\sigma}<\sigma_{\varepsilon}<\bar{\sigma}$. Next we define $g_{\sigma}\left\{{\bf x},({\bf c},r)\right\}=(2\pi\sigma^{2})^{-d/2}\exp\left\{-\left(\|{\bf x}-{\bf c}\|-r\right)^{2}/(2\sigma^{2})\right\}$. Note that $g_{\sigma}\left\{{\bf x},({\bf c},r)\right\}$ is an increasing function of $\sigma$ if $\sigma<\left|\|{\bf x}-{\bf c}\|-r\right|/\sqrt{d}$. To see this consider the following:

[TABLE]

Let $A=\{({\bf c},r):|\|{\bf x}-{\bf c}\|-r|>\sqrt{d}\bar{\sigma}\}$, then $h_{\sigma}\left\{{\bf x},({\bf c},r)\right\}\leq e^{2\tau_{\max}}g_{\bar{\sigma}}\left\{{\bf x},({\bf c},r)\right\}$ for $({\bf c},r)\in A$ by Lemma 2. For $A^{c}=\{({\bf c},r):|\|{\bf x}-{\bf c}\|-r|<\sqrt{d}\bar{\sigma}\}$ we bound the kernel $h_{\sigma}\left\{{\bf x},({\bf c},r)\right\}$ by $g_{\underline{\sigma}}\left\{{\bf x},({\bf c},r)\right\}=(2\pi\underline{\sigma}^{2})^{-d/2}\exp\left\{-\left(\|{\bf x}-{\bf c}\|-r\right)^{2}/(2\bar{\sigma}^{2})\right\}$. Now it is easy to see that $g_{\underline{\sigma}}\left\{{\bf x},({\bf c},r)\right\}$ and $g_{\bar{\sigma}}\left\{{\bf x},({\bf c},r)\right\}$ are $P_{{\bf c},r}$-integrable and $g({\bf x},\boldsymbol{\theta})=\max\left[g^{\prime}\left\{{\bf x},({\bf c},r)\right\},g_{\bar{\sigma}}\left\{{\bf x},({\bf c},r)\right\}\right]<C_{0}$, for some appropriate constant $C_{0}$, depending on $(\underline{\sigma},\bar{\sigma})$. Thus (ii) holds.

Define $D=\{({\bf c},r,\boldsymbol{\pi}):\|c\|<R,\underline{r}<r<\bar{r};\boldsymbol{\pi}\in\Delta_{M}\}$. To prove (iii) we will show that for any $\sigma\in N(\sigma_{\varepsilon})$:

[TABLE]

For any $\sigma\in N(\sigma_{\varepsilon})$, and by Lemma 2

[TABLE]

Now $\sup_{\boldsymbol{\theta}\in D}\sum_{j}\pi_{j}^{({\bf c},r)}h_{\sigma}({\bf x},\boldsymbol{\theta})\leq\sup_{\boldsymbol{\pi}\in\Delta_{M}}\sum_{j}\pi_{j}^{({\bf c},r)}\sup_{({\bf c},r)\in D}h_{\sigma}({\bf x},\boldsymbol{\theta})$. Note that $\sup_{({\bf c},r)\in D}h_{\sigma}({\bf x},\boldsymbol{\theta})$ is no bigger than

[TABLE]

and therefore

[TABLE]

which is $f_{0}$-integrable. Similarly $\left|\log\inf_{\boldsymbol{\theta}\in D}\sum_{j}\pi_{j}h_{\sigma}({\bf x},\boldsymbol{\theta})\right|$ is no bigger than

[TABLE]

The last expression is also $f_{0}$-integrable as $E_{0}\left(\|{\bf x}\|^{2}\right)<\infty$. Thus we dominate

[TABLE]

which is $f_{0}$-integrable for $\sigma\in N(\sigma_{\varepsilon})$. Therefore by the DCT

[TABLE]

Therefore, for any given $\varepsilon>0$, there exists $\delta>0$ such that

[TABLE]

whenever $|\sigma-\sigma_{\varepsilon}|<\delta$. Thus we choose $S=\left\{\sigma:\underline{\sigma},\bar{\sigma}\right\}\cap\left\{\sigma:\sigma_{\varepsilon}-\delta,\sigma_{\varepsilon}+\delta\right\}=\left\{\sigma:\underline{\sigma}^{*},\bar{\sigma}^{*}\right\}$, say. The proof of (II) is complete noticing that $\pi(S)>0$, $S$ is nonempty (it contains a neighborhood of $\sigma_{\varepsilon}$) as $\sigma^{2}$ follows an inverse-gamma prior.

Showing condition (III) holds: To check the last condition we apply the following Lemma 3 of Wu and Ghosal, (2008), which states that (III) holds if the following three conditions hold:

1. i.

For any $\sigma\in S$, $\displaystyle{\int_{\mathbb{R}^{d}}\left|\log\frac{f_{P_{\varepsilon},\sigma}({\bf x})}{\inf_{\boldsymbol{\theta}\in D}\sum_{j}\pi_{j}^{({\bf c},r)}h_{\sigma}({\bf x},\boldsymbol{\theta})}\right|f_{0}({\bf x})d{\bf x}<\infty}$. 2. ii.

Define $c^{*}\coloneqq\inf_{{\bf x}\in C}\inf_{\boldsymbol{\theta}\in D}\sum_{j}\pi_{j}^{({\bf c},r)}h_{\sigma}({\bf x},\boldsymbol{\theta})$ then $c^{*}>0$, for any compact $C\in\mathcal{X}$. 3. iii.

For any $\phi\in S$ and compact $C\in\mathcal{X}$, there exists $E$ containing $D$ in its interior such that the family of maps $\{\boldsymbol{\theta}\mapsto\sum_{j}\pi_{j}^{({\bf c},r)}h_{\sigma}({\bf x},\boldsymbol{\theta}),{\bf x}\in C\}$ is uniformly equicontinuous on $E\subset\Theta$, and $\sup\left\{\sum_{j}\pi_{j}^{({\bf c},r)}h_{\sigma}({\bf x},\boldsymbol{\theta}):{\bf x}\in C,\theta\in E^{c}\right\}<k\varepsilon/4$, for some $k>0$.

To check i., we show that

[TABLE]

and

[TABLE]

The second part is already shown in the proof of condition (II). To see the first inequality, we proceed in a similar way as in the proof of (I). Observe that for fixed $\sigma\in S$

[TABLE]

where $C_{\sigma}=(2\pi\sigma^{2})^{-d/2}$. Again observe that

[TABLE]

by Lemma 2. Next note that for $\theta\in D$, $\underline{r}\leq m^{-\eta}_{\varepsilon}\leq\bar{r}$. Therefore,

[TABLE]

By assumption D, there exists some $\delta>0$ such that $\inf_{\|{\bf x}\|<\delta}f_{0}({\bf x})=\phi_{\delta}>0$. Further recall that $t_{m}\geq 1$, therefore

[TABLE]

For the chosen fixed $\delta>0$ in assumption D, the last term is positive. Therefore (i) follows noting that $E_{0}\left(\|{\bf x}\|^{2}\right)<\infty$.

From (14) we get:

[TABLE]

Note that $g({\bf x})$ is positive, bounded and continuous on ${\bf x}\in C$, where $C$ is any compact set. Hence (ii) holds.

It remains to show that (iii) is satisfied. Let $C\in\mathcal{X}$ be a given compact set. We need to show that $h_{\sigma}({\bf x},\boldsymbol{\theta})$ is uniformly equicontinuous as a family of functions of $\boldsymbol{\theta}=({\bf c},r)$ on a set $E\supset D$. We choose $E=\{({\bf c},r,\boldsymbol{\pi}):\|{\bf c}\|<a;~{}r<\bar{r}^{*};~{}\boldsymbol{\pi}\in\Delta_{M}\}$ where $0<\bar{r}<\bar{r}^{*}<\infty$, $a>R$. Clearly, $E$ contains $D$ and is compact. By the definition equicontinuouity, for all $\boldsymbol{\theta}_{1}=({\bf c}_{1},r_{1},\boldsymbol{\pi}_{1}),\boldsymbol{\theta}_{2}=({\bf c}_{2},r_{2},\boldsymbol{\pi}_{2})\in E$ with $\|({\bf c}_{1},r_{1},\boldsymbol{\pi}_{1})-({\bf c}_{2},r_{2},\boldsymbol{\pi}_{2})\|<\delta_{\epsilon}$, we have $\left|\sum_{j}\pi_{1,j}^{({\bf c}_{1},r_{1})}h_{\sigma}({\bf x},\boldsymbol{\theta}_{1})-\sum_{j}\pi_{2,j}^{({\bf c}_{2},r_{2})}h_{\sigma}({\bf x},\boldsymbol{\theta}_{2})\right|<\epsilon$.

For any $\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}$,

[TABLE]

Thus it is enough to show that $\left|h_{\sigma}({\bf x},\boldsymbol{\theta}_{1})-h_{\sigma}({\bf x},\boldsymbol{\theta}_{2})\right|\leq\epsilon_{1}/M$ for some appropriate $\epsilon_{1}<\epsilon$ and $h_{\sigma}({\bf x},\boldsymbol{\theta}_{2})$ is bounded for ${\bf x}\in C$ and $\boldsymbol{\theta}\in E$. For the first part, we will show that the ratio $h_{\sigma}({\bf x},\boldsymbol{\theta}_{1})\left/h_{\sigma}({\bf x},\boldsymbol{\theta}_{2})\right.$ goes to $1$ uniformly, as

[TABLE]

Now

[TABLE]

by the multivariate mean value theorem where ${\bf c}^{*}=t{\bf c}_{1}+(1-t){\bf c}_{2}$, $0<t<1$. The constant $M$ depends on the upper bound of $\|{\bf x}\|$ and $a$ only. Similarly, $\left|r_{1}^{2}-r_{2}^{2}\right|=(r_{1}+r_{2})|r_{1}-r_{2}|<2\bar{r}^{*}\delta_{\epsilon}$.

Let $\nu\coloneqq d/2-1$ and observe that

[TABLE]

Note that $C_{d}^{-1}(\tau)$ is a continuous, strictly increasing function of $\tau$ (see Baricz, (2010)). Further,

[TABLE]

Next note that $g(\cdot)$ is a continuous function on a compact set, and hence the image of $g(\cdot)$ is compact. Thus, $C_{d}(\cdot)$ is uniformly continuous on the image of $g(\cdot)$, and we can write

[TABLE]

where $\xi_{\delta}\rightarrow 0$ as $\delta_{\epsilon}\rightarrow 0$. Therefore, $h_{\sigma}({\bf x},\boldsymbol{\theta}_{1})/h_{\sigma}({\bf x},\boldsymbol{\theta}_{2})\leq\exp\left(M^{*}\delta_{\epsilon}\right)$ for some suitable constant $M^{*}$. Similarly, it can be shown that $h_{\sigma}({\bf x},\boldsymbol{\theta}_{1})/h_{\sigma}({\bf x},\boldsymbol{\theta}_{2})\geq\exp\left(-M^{*}\delta_{\epsilon}\right)$. Hence the above uniformly converges to 1 as $\delta_{\epsilon}\rightarrow 0$.

Again, from Lemma 2  $\exp\left(-U/\sigma^{2}\right)\leq\left(2\pi\sigma^{2}\right)^{d/2}h_{\sigma}({\bf x},\boldsymbol{\theta}_{2})\leq\exp\left(U/\sigma^{2}\right),$  where $U=\left(\|{\bf x}\|^{2}+a^{2}+\bar{r}^{*2}+2\tau_{m}\right)$. Hence $h_{\sigma}({\bf x},\boldsymbol{\theta}_{2})$ is bounded away from zero and bounded above in ${\bf x}\in C$ and $\boldsymbol{\theta}_{2}\in E$. This completes the proof of equicontinuouity.

The proof is complete if we show that $\sup\left\{\sum_{j}\pi_{j}^{({\bf c},r)}h_{\sigma}({\bf x},\boldsymbol{\theta}):{\bf x}\in C,\boldsymbol{\theta}\in E^{c}\right\}<k\varepsilon/4$. For that we just show $\sup\left\{h_{\sigma}({\bf x},\boldsymbol{\theta}):{\bf x}\in C,F\right\}<k\varepsilon/4$, where $F=\{({\bf c},r):\|c\|>a;~{}\mbox{or}~{}r>\bar{r}^{*}\}$. As under $E$, $\boldsymbol{\pi}$ is any point in $\Delta_{M}$, $E^{c}=\left(F\cup\{\phi\}\right)$ and the result follows.

Let $\|{\bf x}\|<m^{*}$ if ${\bf x}\in C$. Also, we can choose $a$ such that whenever $\|{\bf c}\|>a$ then $\|{\bf x}-{\bf c}\|\geq M_{\varepsilon}$. We split the space $\{({\bf x},\boldsymbol{\theta}):{\bf x}\in C;({\bf c},r)\in F\}$ into two parts:

[TABLE]

Here $C_{\sigma}=e^{-2\tau_{\max}}(2\pi\sigma)^{d/2}$. Further we consider $\varepsilon$ to be small enough so that $C_{\sigma}k\varepsilon<1$. Observe that

[TABLE]

For $B$, first note that when $\boldsymbol{\theta}\in E^{c}$, either $\|{\bf c}\|>a$ or $r>\bar{r}^{*}$, or both happen. Let $\tau_{m}\geq\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|$. Note that $C_{d}^{-1}(\tau_{m})$ is a strictly increasing function of $\tau_{m}$. Therefore, $C_{d}^{-1}(\tau_{m})\geq C_{d}^{-1}\left(\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|\right),$ implying $C_{d}(\tau_{m})/C_{d}\left(\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|\right)\leq 1.$ Further, as $\left({\bf c},r\right)\in F$, for fixed $\sigma^{2}$ and suitable choices of $a$ and $\bar{r}^{*}$, $\left(\|{\bf x}-{\bf c}\|^{2}+r^{2}\right)/2\sigma^{2}>-\log\left(C_{\sigma}k\varepsilon\right)$, $\sup_{({\bf x},\boldsymbol{\theta})\in A}h_{\sigma}({\bf x},\boldsymbol{\theta})\leq k\varepsilon$.

Next consider $\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|>\tau_{m}$. By Baricz, (2010, Equation (2.10)) for $0<\tau_{m}<\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|$,

[TABLE]

Therefore, $\sup_{({\bf x},\boldsymbol{\theta})\in B}h_{\sigma}({\bf x},\boldsymbol{\theta})$ is no bigger than

[TABLE]

For sufficiently large $M_{\varepsilon}$ and $\bar{r}^{*}$ we can write the bracketed portion in the last expression as $\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|+d/2-1\geq c^{*}r\|{\bf x}-{\bf c}\|/\sigma^{2}$, for some suitable constant $c^{*}$. Note that under $B$, $\left(r-\|{\bf x}-{\bf c}\|\right)^{2}<-2\sigma^{2}\log\left(C_{\sigma}k\varepsilon\right)=\delta_{\varepsilon}$, say, implying $r\|{\bf x}-{\bf c}\|>\left(r^{2}+\|{\bf x}-{\bf c}\|^{2}-\delta_{\varepsilon}\right)/2$. Therefore, we can expand $E$ by choosing $a$ and $\bar{r}^{*}$ such that the last expression is less than $k\varepsilon$ for $d>2$.

Finally consider the case with $d=2$. Fix a small $\epsilon=1/4$. By (Baricz,, 2010, Theorem 2.2) $I_{0}(\tau)/I_{-0.5+\epsilon}(\tau)$ is a strictly decreasing function of $x$. So for any $\tau>1$,

[TABLE]

Therefore, for $d=2$,

[TABLE]

where $M=c^{*}I_{-0.5+\epsilon}(\tau_{m})/I_{0}(\tau_{m})>0$ as without loss of generality for $({\bf c},r)\in B$, we can assume that $\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|>1$. Next by Joshi and  , (1991), as $\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|>\tau_{m}$, $I_{-0.5+\epsilon}\left\{\|\tau_{m}\boldsymbol{\mu}_{m}+r({\bf x}-{\bf c})/\sigma^{2}\|\right\}/I_{-0.5+\epsilon}(\tau_{m})$ is no bigger than

[TABLE]

Thus $\sup_{({\bf x},\boldsymbol{\theta})\in B}h_{\sigma}({\bf x},\boldsymbol{\theta})$ is no bigger than

[TABLE]

for an appropriate constant $M^{*}$. As $\epsilon=1/4$, we can enlarge $E$ to show that last expression is less than $k\varepsilon/4$ for all $\tau_{m}$, $m=1,\ldots,M$. This completes the proof. ∎