# Goodness-of-Fit and Clustering of Spherical Data: the QuadratiK package in R
and Python

Giovanni Saraceno  Department of Biostatistics, University at Buffalo, NY, USA
Marianthi Markatou  Department of Biostatistics, University at Buffalo, NY,
USA  Raktim Mukhopadhyay  Department of Biostatistics, University at Buffalo,
NY, USA  CDSE Program
University at Buffalo  Mojgan Golzy  Biomedical Informatics, Biostatistics and
Medical Epidemiology (BBME), University of Missouri

###### Abstract

We introduce the QuadratiK package that incorporates innovative data analysis
methodologies. The presented software, implemented in both R and Python,
offers a comprehensive set of goodness-of-fit tests and clustering techniques
using kernel-based quadratic distances, thereby bridging the gap between the
statistical and machine learning literatures. Our software implements one, two
and k-sample tests for goodness of fit, providing an efficient and
mathematically sound way to assess the fit of probability distributions.
Expanded capabilities of our software include supporting tests for uniformity
on the $d$-dimensional Sphere based on Poisson kernel densities, and
algorithms for generating random samples from Poisson kernel densities.
Particularly noteworthy is the incorporation of a unique clustering algorithm
specifically tailored for spherical data that leverages a mixture of Poisson-
kernel-based densities on the sphere. Alongside this, our software includes
additional graphical functions, aiding the users in validating, as well as
visualizing and representing clustering results. This enhances
interpretability and usability of the analysis. In summary, our R and Python
packages serve as a powerful suite of tools, offering researchers and
practitioners the means to delve deeper into their data, draw robust
inference, and conduct potentially impactful analyses and inference across a
wide array of disciplines.

Keywords: kernel-based quadratic distances, goodness-of-fit, k-sample tests,
two-sample test, clustering on sphere

##  1 Introduction

††footnotetext: ∗Corresponding author

The Goodness-of-Fit (GoF) problem is crucial in statistical research for
several reasons. Firstly, GoF tests constitute a classical tool in deciding
the compatibility of data with an assumed probability distribution. Secondly,
comparing two or more samples allows researchers to examine differences or
similarities between groups or conditions. This can help identify factors that
may contribute to certain outcomes or behaviours. Furthermore, such
comparisons facilitate the validation of findings and ensure the
generalizability of research results, thus enhance the robustness and
reliability of research findings. The standard way of constructing a GoF test
procedure involves the computation of some distance-like functional between
the null distribution and the given observations, where the null hypothesis is
rejected if the obtained distance is greater than a critical value. Some of
the commonly used tests in this category are the Kolmogorov-Smirnov test,
performed through the R function ks.test in the stats package, the Cramér-von
Mises, in the R package goftest [Faraway et al., 2019], and Anderson-Darling
test, available in the R package goft [Gonzalez-Estrada and Villasenor-Alva,
2017]. These tests are also available in the Python library scipy [Virtanen et
al., 2020]. The R package vsgoftest performs GoF tests of various families of
distributions (uniform, normal, lognormal, exponential, gamma, Weibull,
Pareto, Fisher, Laplace and Beta) based on Shannon entropy and the Kullback-
Leibler divergence. These tests were developed by Vasicek [1976] and Song
[2002]. The Vasicek-Song (VS) tests are intended to be applied to univariate
data obtained as realizations of continuous random variables, and, in case of
multivariate data, only the multivariate normality test is included. The
dbEmpLikeGOF package developed by Miecznikowski et al. [2013] proposes GoF
normality and uniformity tests, and two-sample test, based on empirical
likelihood ratio for univariate samples. These tests are closely related to VS
tests. The GoFKernel package contains an implementation of the Fan’s test
[Fan, 1994], which is based on the L2-distance, and a practical approximation,
using the L1-norm, to compare the discrepancy between a theoretical density
function and a sample kernel estimate of the density function [Cao and Lugosi,
2005].

The two-sample problem has a significant body of literature, and many
univariate GoF testing procedures, have been suggested. The R package
robnptests proposed by Abbas et al. [2023] offers a compilation of univariate
two-sample tests with respect to differences in location and dispersion.
However, relatively less attention has been given to the multivariate two-
sample case. Some of the available methods for multivariate testing utilize
graph-based multivariate ranking strategies for the extensions of univariate
two-sample GoF tests in high dimensional settings. One such approach was
introduced by Friedman and Rafsky [1979], who employed the minimal spanning
tree (MST) as a multivariate ordered list. They extended the univariate run-
based test introduced by Wald and Wolfowitz [1940], the univariate two-sample
Kolmogorov-Smirnov test and a modified Kolmogorov-Smirnov test against scale
alternatives. These methods are available in the R package GSAR. Similarly,
Biswas et al. [2014] adopted the shortest Hamiltonian path instead of MST as
ranking strategy for multivariate observations. Chen and Friedman [2017]
presented a test statistic based on a similarity graph. In addition to graph-
based methods, general distance measures have been widely adopted to handle
the multivariate two-sample GoF testing problem. Rosenbaum [2005] introduced a
two-sample test based on interpoint distances, available in the R package
crossmatch. The R package energy [Rizzo and Szekely, 2022], proposed by Rizzo
and Székely [2016], presents a collection of test statistics for multivariate
and univariate inference based on the Energy-statistics. This includes
normality tests, two- and multi-sample tests, a test for independence and
hierarchical clustering. From the community of Machine Learning, Gretton et
al. [2012] constructed the test statistic called Maximum Mean Discrepancy
(MMD), utilizing the properties of the kernel mean embedding. This test can be
performed through the function kmmd in the R package kernlab. The Friedman-
Rafsky, MMD and energy tests are also provided in the Python package hyppo
[Panda et al., 2019a]. The operating characteristics of several of the tests
mentioned above, were studied in Chen and Markatou [2020].

The two-sample problem can be extended to the $k$-sample testing problem,
where more than two groups are compared. The Anderson-Darling test, initially
introduced for the case of testing equality of two distributions, was extended
by Scholz and Stephens [1987] and it is available in the R package kSamples
[Scholz and Zhu, 2019]. This package includes several nonparametric Rank Score
$k$-sample tests, that is Kruskal-Wallis [Kruskal and Wallis, 1952], van der
Waerden scores and normal scores [Sidak et al., 1999]. An alternative to
nonparametric test construction is given by the general framework of
permutation tests. The R package coin [Hothorn et al., 2008] provides an
implementation of such tests tailored against location and scale alternatives,
and for survival distributions. Rizzo and Székely [2010] proposed a
nonparametric extension of the standard analysis of variance procedure, called
Distance Components (DISCO), using the Energy statistics. This is available in
the package energy. Kernel-based methods have gained popularity also for the
$k$-sample problem due to their ability to handle high-dimensional data.
Balogoun et al. [2022] proposed a $k$-sample test for functional data by
introducing a generalization of the MMD. Recently, Panda et al. [2019b]
highlighted the connection between $k$-sample testing and independence testing
problems showing that independence tests, such as the MMD and the Energy
statistics, can be used for consistent $k$-sample testing. These methods are
available in the Python package hyppo.

In many applications of interest, data can be represented as points in some
non-Euclidean space or manifold. This is the case of directional data, which
naturally arise in problems where the directions of the data and not the
magnitudes are of interest. Directional observations may be regarded as points
on the surface of a unit hypersphere
$\mathcal{S}^{d-1}=\\{x\in\mathbb{R}^{d}:||x||^{2}=1\\}$. Furthermore, the
analysis of many non-directional data sets may be enhanced by re-expressing
the data in the form of directions [Golzy and Markatou, 2020, Dortet-Bernadet
and Wicker, 2008], and determining uniformity is a critical first step when
analyzing such data. There is a considerable amount of work on testing for
uniformity for circular data. For example the R package circular [Agostinelli
and Lund, 2023] provides Rayleigh’s test, Rao’s Spacing test, Kuiper’s and the
Watson’s tests of uniformity. Fernández-Durán and Gregorio-Domínguez [2016]
propose a test for uniformity based on nonnegative trigonometric sums,
available in the R package CircNNTSR [Fernández-Durán and Gregorio-Domínguez,
2016]. For dimensions $d>2$, the class of Sobolev tests, introduced by Giné
[1975], contains most of the nonparametric tests on $\mathcal{S}^{d-1}$. Some
notable examples belonging to this class are Rayleigh’s test [Rayleigh, 1919],
Watson’s test, Ajné test [Ajne, 1968], Bingham’s test [Bingham, 1974]. The
Rayleigh’s test can be performed using the function rayleigh in the R package
Directional [Tsagris et al., 2023]. Since Sobolev tests depend on the choice
of a sequence of coefficients, Jupp [2008], Bogdan et al. [2002] were
motivated to introduce data-driven Sobolev tests. Cuesta-Albertos et al.
[2009] proposed a test based on random projection for assessing uniformity on
$\mathcal{S}^{d-1}$. The implementation of this test is available in the
package Directional. It is worth mentioning the R package sphunif [García-
Portugués and Verdebout, 2023], which includes a collection of almost all the
available methods for testing uniformity on the Circle, Sphere and
Hypersphere. See García-Portugués and Verdebout [2018] and Section 5.1 in
Pewsey and García-Portugués [2021] for a review of recent advances on tests of
uniformity.

Our interest focuses on tests that are based on the family of kernel-based
quadratic distances (KBQDs). Quadratic distances are central to the study of
goodness of fit, and important statistics that belong to this family are the
Pearson’s Chi-squared and Cramer-von Mises statistics. Lindsay et al. [2008]
introduced a unified framework for the study of quadratic distance measures,
and Lindsay et al. [2014] studied the power properties of kernel-based
quadratic distance goodness-of-fit tests introducing the midpower analysis and
the concept of degrees of freedom for the kernel function. The authors, then,
investigate the power properties in case of testing for multivariate normality
by using the normal kernel. Markatou et al. [2021] offer a comprehensive
review of the quadratic distance-based methodologies, while Markatou et al.
[2017] offer interesting interpretations of distances and discuss their role
in robustness. Additionally, Markatou et al. [2023] introduce a unified
framework for constructing two- and $k$-sample hypothesis tests using kernel-
based quadratic distances. The kernels we use are diffusion kernels [Ding et
al., 2023], that is, probability distributions that depend on a tuning
parameter and satisfy the convolution property. We also note here that the
tuning parameter of the kernel is intimately connected with the degrees of
freedom (DOF) of a kernel and hence with the power characteristics of the
tests based on these kernels. By choosing the Poisson kernel, Ding et al.
[2023] proposed the Poisson kernel-based tests for uniformity on the
$d$-dimensional sphere. Indeed, the Poisson kernel appears to be the natural
generalization of the normal distribution to the sample space
$\mathcal{S}^{d-1}$. Additionally, through a normalizing constant, the Poisson
kernel can be considered as a density on the sphere. Applications of the
Poisson kernel-based densities to clustering can be found in Golzy and
Markatou [2020] and Golzy et al. [2016].

In this article, we present the package QuadratiK, a comprehensive package of
methods for multivariate analysis and inference that are based on quadratic
distances. These include the one-sample, two-sample and $k$-sample goodness-
of-fit tests, implemented through the functions pk.test and kb.test, and the
clustering procedure introduced by Golzy and Markatou [2020] with the function
pkbc. The implemented methods are practical and accessible, and they can
handle high dimensional data sets. In particular, the test procedures employ
optimized C++ versions which speed up the calculations. The implementation of
the quadratic distance tests involves the following steps. First, the kernels
need to be centered. Centering of the kernel is done both parametrically and
nonparametrically. Second, the computation of the critical value,
incorporating various sampling algorithms, and the selection of the tuning
parameter of the kernel are important aspect for computing the tests. Indeed,
the software offers the possibility to select the optimal value of the tuning
parameter according to the midpower analysis. By employing the use of parallel
computing, the algorithm exhibits a feasible computational time. The
clustering algorithm is based on the Poisson kernel-based densities. The
software includes the functions dpkb and rpkb, which are used for computing
the density value and generating random observations from the Poisson kernel-
based distribution. Additionally, in order to facilitate the interpretability
and usability of the clustering results, the package includes graphical
functions, which aid the users in validating and visualizing clustering
results. The presented software is implemented as R and Python packages in
order to enhance their adoption into pre-existing workflows by users from
statistics and computer science communities thus enhancing the dissemination
and use of these tools. Finally, we propose a dashboard application, a
graphical user interface to the implemented methods, with the aim to
facilitate the usage of the presented software among practitioners from
different than the above two mentioned domains.

The rest of the article is organized as follows. In Section 2, the one-sample
kernel-based quadratic distance test statistic is presented. In the
corresponding subsections, the derived test statistic is used for constructing
the multivariate normality test [Lindsay et al., 2014] and for testing
uniformity on the Sphere [Ding et al., 2023], under appropriate choices of the
kernel function. Each method in the presented package is introduced together
with a synthetic example with codes both in R and Python. Subsection 2.2
introduces the non-parametric two-sample and $k$-sample KBQD test statistic
for comparing two or more multivariate distributions [Markatou et al., 2023].
This includes the algorithm proposed for selecting the bandwidth parameter $h$
of the normal kernel based on the power analysis of the test. The clustering
algorithm appropriate for use with data on the sphere [Golzy and Markatou,
2020] provided in the software is described in Section 4, and its usage,
together with the additional graphical functions, is shown through a real data
application. In addition, section 4.1 presents the implementation of different
algorithms for generating random samples following the Poisson kernel-based
density. Section 5 illustrates the Dashboard application, while a summary is
provided in Section 6.

##  2 Kernel-based Quadratic Distances

Given two probability distributions $F$ and $G$, the KBQD $d_{K}(F,G)$ is
defined as

$d_{K}(F,G)=\iint
K_{G}(\boldsymbol{s},\boldsymbol{t})d(F-G)(\boldsymbol{s})d(F-G)(\boldsymbol{t}),$
(1)

where $K_{G}(\boldsymbol{s},\boldsymbol{t})$ denotes a non-negative definite
kernel function that possibly depends on $G$. Uniqueness considerations of the
value of distance associated with the kernel on which the distance is based
dictate the use of centered kernels. Thus, the $G$-centered kernel, where $G$
is the probability distribution to be fitted to the data is defined as

$K_{\mathrm{cen}(G)}(\boldsymbol{s},\boldsymbol{t})=K(\boldsymbol{s},\boldsymbol{t})-K(\boldsymbol{s},G)-K(G,\boldsymbol{t})+K(G,G),$

with $K(\boldsymbol{s},G)=\int
K(\boldsymbol{s},\boldsymbol{t})dG(\boldsymbol{t})$, $K(G,\boldsymbol{t})=\int
K(\boldsymbol{s},\boldsymbol{t})dG(\boldsymbol{s})$, $K(G,G)=\iint
K(\boldsymbol{s},\boldsymbol{t})dG(\boldsymbol{s})dG(\boldsymbol{t})$. A
cornerstone in constructing and using KBQDs is the choice of the kernel
function. We concentrate on the special class of diffusion kernels, introduced
in Ding et al. [2023]. A family of symmetric kernels
$K_{t}(\boldsymbol{s},\boldsymbol{t})$ with parameter $t\in(0,\infty)$ is
called a diffusion kernel family with respect to a measure $\gamma$ if the
following properties are satisfied.

  * •

The kernel is nonnegative, i.e. $K_{t}(\boldsymbol{s},\boldsymbol{t})\geq 0$;

  * •

The kernel satisfies the diffusion equation

$K_{t_{1}+t_{2}}(\boldsymbol{s},\boldsymbol{t})=\int
K_{t_{1}}(\boldsymbol{s},\boldsymbol{r})K_{t_{2}}(\boldsymbol{r},\boldsymbol{t})dr;$
(2)

  * •

The kernel is a probability density under the measure $\gamma$ with respect to
both arguments, that is

$\int K_{t}(\boldsymbol{s},\boldsymbol{t})d\gamma(\boldsymbol{s})=1\qquad\int
K_{t}(\boldsymbol{s},\boldsymbol{t})d\gamma(\boldsymbol{t})=1.$

This class of kernels allows easy computation. Furthermore, the kernels are
”tunable”, since they depend on a tuning parameter $t$, connected with their
degrees of freedom (DOF). The DOF of the kernel, in turn, drive the power
characteristics of tests based on KBQDs (see Lindsay et al. [2008, 2014],
Markatou et al. [2023] and Ding et al. [2023] for further information).

###  2.1 Goodness-of-fit Tests

If $G$ denotes a distribution whose goodness of fit we wish to assess, the
kernel-based quadratic distance can be employed for constructing a test
statistic by measuring the distance between the sample data and the target
distribution. This is known as the one-sample case. Consider $H_{0}:F=G$
versus $H_{1}:F\not=G$ and let
$\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n}\sim F$ independent identically
distributed (i.i.d.) observations, with cumulative distribution function
(c.d.f.) $\hat{F}$. By centering $K$ with respect to the distribution $G$, the
KBQD can be written as [Lindsay et al., 2014]

$d_{K}(\hat{F},G)=\iint
K_{cen(G)}(\boldsymbol{s},\boldsymbol{t})d\hat{F}(\boldsymbol{s})d\hat{F}(\boldsymbol{t}).$
(3)

Then, the corresponding $V$\- and $U$-statistics are given by

$V_{n}=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=i}^{n}K_{cen(G)}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\qquad\mbox{and}\qquad
U_{n}=\frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j\not=i}^{n}K_{cen(G)}(\boldsymbol{x}_{i},\boldsymbol{x}_{j}).$
(4)

Under the simple null hypothesis, the asymptotic distribution of $U_{n}$ is

$nU_{n}\longrightarrow\sum_{j}\lambda_{j}(Z_{j}^{2}-1)\mbox{ as }n\to\infty,$

where $\lambda_{j}$s are the nonzero eigenvalues of the centered kernel
$K_{G}$ under the null distribution G, and $Z_{j}$ are independent standard
normal random variables.

The two-sample test is presented here as a spacial case of the $k$-sample
test.

###  2.2 $k$-Sample Tests

The goodness-of-fit tests based on a kernel-based quadratic distance can be
extended to the $k$-sample problem via the employment of a matrix distance,
introduced by Markatou et al. [2023].

Consider the null hypothesis $H_{0}:F_{1}=F_{2}=\ldots=F_{k}$ against the
alternative $H_{1}:F_{i}\not=F_{j}$, for some $1\leq i\not=j\leq k$. Assume
that under the null hypothesis $F_{1}=\ldots=F_{k}=\bar{F}$, then the $ij$-th
element of the matrix distance $D=(D_{ij})$ is defined as

$D_{ij}=\iint
K(\boldsymbol{x},\boldsymbol{y})d(F_{i}-\bar{F})(\boldsymbol{x})d(F_{j}-\bar{F})(\boldsymbol{y}).$

Consider random samples of i.i.d. observations
$\boldsymbol{x}^{(i)}_{1},\boldsymbol{x}^{(i)}_{2},\ldots,\boldsymbol{x}^{(i)}_{n_{i}}\sim
F_{i}$, let $\hat{F}_{i}$ be the corresponding c.d.f. and
$n=\sum_{i=1}^{K}n_{i}$. By centering the kernel function with respect to a
chosen distribution $\bar{F}$, we have that

$\hat{D}_{ij}=\iint
K_{cen(\bar{F})}(\boldsymbol{x},\boldsymbol{y})d\hat{F}_{i}(\boldsymbol{x})d\hat{F}_{j}(\boldsymbol{y}).$

Then, the empirical version of the matrix distance $\hat{D}_{n}$ has as
elements the $U$-statistics

$\hat{D}_{ij}=\frac{1}{n_{i}n_{j}}\sum_{l=1}^{n_{i}}\sum_{r=1}^{n_{j}}K_{cen(\bar{F})}(\boldsymbol{x}^{(i)}_{l},\boldsymbol{x}^{(j)}_{r})\qquad\mbox{
for }i\not=j$ (5)

and

$\hat{D}_{ii}=\frac{1}{n_{i}(n_{i}-1)}\sum_{l=1}^{n_{i}}\sum_{r\not=l}^{n_{i}}K_{cen(\bar{F})}(\boldsymbol{x}^{(i)}_{l},\boldsymbol{x}^{(i)}_{r})\qquad\mbox{
for }i=j.$ (6)

If we want to test the equality of $k$ samples against the general alternative
$H_{1}:F_{i}\not=F_{j}$, for some $1\leq i\not=j\leq k$, it is possible to
construct an omnibus test considering all the possible equalities included in
the $k$-sample null hypothesis. This can be derived from the matrix distance
as

$T_{n}=\mathrm{trace}(\hat{D}_{n})-\frac{2}{K-1}\sum_{i=1}^{K}\sum_{j>i}^{K}\hat{D}_{ij}.$
(7)

In our implementation, we provide the values of $T_{n}$ and $(k-1)T_{n}$, with
the associated critical values.

For the special case $K=2$, the matrix distance gives the test statistic

$\displaystyle D_{{n_{1}},n_{2}}=$
$\displaystyle\frac{1}{{n_{1}}({n_{1}}-1)}\sum_{i=1}^{n_{1}}\sum_{j\not=i}^{n_{1}}K_{cen(\bar{F})}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})-\frac{2}{{n_{1}}{n_{2}}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}K_{cen(\bar{F})}(\boldsymbol{x}_{i},\boldsymbol{y}_{j})$
$\displaystyle+\frac{1}{{n_{2}}({n_{2}}-1)}\sum_{i=1}^{n_{2}}\sum_{j\not=i}^{n_{2}}K_{cen(\bar{F})}(\boldsymbol{y}_{i},\boldsymbol{y}_{j}).$

For further theoretical information see Markatou et al. [2023].

##  3 Multivariate KBQD tests

In this section, we deal with the computational aspects of the kernel-based
quadratic distance tests. These include the choice of the kernel function, the
selection of the corresponding tuning parameter, and the computation of the
critical values for assessing the hypothesis testing problem. Then, we
introduce the main functions offered by the QuadratiK package which perform
the presented multivariate KBQD test procedures. We will show their usage
through simulated data examples.

###  3.1 Kernel computation

Our kernels belong in the general class of diffusion kernels. This is because
tunable diffusion kernels allow for easy computation. Furthermore, we take
into consideration the type of sample space our data lie, that is Euclidean
versus a general manifold.

For multivariate data in Euclidean space, the KBQD are constructed considering
the multivariate normal kernel. For
$\boldsymbol{s},\boldsymbol{t}\in\mathbb{R}^{d}$ and covariance matrix
$\boldsymbol{\Sigma}_{h}$, depending on the tuning parameter $h$, it is given
as

$K_{{h}}(\boldsymbol{s},\boldsymbol{t})=(2\pi)^{-d/2}\left(\det{\boldsymbol{\Sigma}_{h}}\right)^{-\frac{1}{2}}\exp\left\\{-\frac{1}{2}(\boldsymbol{s}-\boldsymbol{t})^{\top}\boldsymbol{\Sigma}_{h}^{-1}(\boldsymbol{s}-\boldsymbol{t})\right\\}.$

In our current implementation, we employ a covariance matrix
$\boldsymbol{\Sigma}_{h}=h^{2}\boldsymbol{I}_{d}$, where $\boldsymbol{I}_{d}$
represents the $d$-dimensional identity matrix. It is worth noting that more
complex variance-covariance structures can be explored. For instance, one
potential extension involves considering $\boldsymbol{\Sigma}_{h}$ as a
diagonal matrix with diagonal elements $(h_{1},h_{2},\ldots,h_{d})$. This
approach introduces additional parameters, $h_{1}$ through $h_{d}$, which
allow for capturing different scales along each dimension of the data. Such
flexibility can be valuable in scenarios where the underlying data exhibits
varying degrees of correlation or dispersion across dimensions. While our
current implementation focuses on the simpler
$\boldsymbol{\Sigma}_{h}=h^{2}\boldsymbol{I}_{d}$ case, the introduced
extension is under consideration for future versions of our package, for
potentially enhancing the modeling capabilities and adaptability.

For testing normality, given the sample of observations
$\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n}\sim F$, we want to assess if
$F=G$ with $G\sim N_{d}(\boldsymbol{\mu}_{G},\boldsymbol{\Sigma}_{G})$ for
some mean vector $\boldsymbol{\mu}_{G}$ and covariance matrix
$\boldsymbol{\Sigma}_{G}$, and center the kernel function with respect to $G$.
The normal kernel centered with respect to $G$ can be computed as

$\displaystyle K_{G}(\boldsymbol{s},\boldsymbol{t})=$ $\displaystyle
K_{\boldsymbol{\Sigma_{h}}}(\boldsymbol{s},\boldsymbol{t})-K_{\boldsymbol{\Sigma_{h}}+\boldsymbol{\Sigma}_{G}}(\boldsymbol{\mu}_{G},\boldsymbol{t})$
$\displaystyle-
K_{\boldsymbol{\Sigma_{h}}+\boldsymbol{\Sigma}_{G}}(\boldsymbol{s},\boldsymbol{\mu}_{G})+K_{\boldsymbol{\Sigma_{h}}+2\boldsymbol{\Sigma}_{G}}(\boldsymbol{\mu}_{G},\boldsymbol{\mu}_{G}).$

The package offers also a method for analysis of spherical data, i.e.
$\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n}\in S^{d-1}$. In this case, the
Poisson kernel appears to be the natural generalization of the normal
distribution on the $d$-dimensional sphere. For
$\boldsymbol{u},\boldsymbol{v}\in{S}^{d-1}$ and concentration parameter
$\rho$, with $0<\rho<1$, the Poisson kernel is given as

$K_{{\rho}}(\boldsymbol{u},\boldsymbol{v})=\frac{1-\rho^{2}}{(1+\rho^{2}-2\rho(\boldsymbol{u}\cdot\boldsymbol{v}))^{\frac{d}{2}}}.$
(8)

The Poisson kernel is the basis for the development of test statistics that
assess the goodness-of-fit of $H_{0}:F=G$, where $G=\mathcal{U}(S^{d-1})$
denotes the Uniform distribution on the sphere. It belongs to the class of
diffusion kernels; centered with respect to the uniform distribution
$\mathcal{U}(S^{d-1})$ is computed as

$K_{cen(G)}(\boldsymbol{u},\boldsymbol{v})=K_{\rho}(\boldsymbol{u},\boldsymbol{v})-1.$
(9)

For the $k$-sample and two-sample problems, we would prefer to not make any
distributional assumption on the (unknown) common distribution $\bar{F}$, that
is $F_{1}=\ldots,F_{K}=\bar{F}$. In order to construct non-parametric KBQD
test statistics, the kernel function is centered with respect to the average
distribution $\bar{F}=\frac{1}{n}\sum_{i=1}^{K}n_{i}F_{i}$, with
$n=\sum_{i=1}^{K}n_{i}$. Let $\boldsymbol{z}$ denote the pooled sample, then
the non-parametric centered kernel can be computed as

$\displaystyle K_{cen(\bar{F})}(\boldsymbol{x},\boldsymbol{y})=$
$\displaystyle
K(\boldsymbol{x},\boldsymbol{y})-\frac{1}{n}\sum_{i=1}^{n}K(\boldsymbol{x},\boldsymbol{z}_{i})-\frac{1}{n}\sum_{i=1}^{n}K(\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle+\frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j\not=i}^{n}K(\boldsymbol{z}_{i},\boldsymbol{z}_{j}).$

To be able to obtain the critical value of the test is not necessary to know
the eigen-decomposition of the centered kernel. The package provides the non-
parametric calculation of the critical value, with the possibility to choose
from subsampling, bootstrap and permutation sampling procedures. Algorithm 1
presents the critical value computation.

Let $G$ denote the true (common) distribution under the null hypothesis. Let
$B$ denote the number of sampling repetitions.

1 Generate i.i.d. observations
$\boldsymbol{z}_{1},\ldots,\boldsymbol{z}_{n_{B}}$ from $G$, or from $\bar{F}$
following the chosen sampling method (subsampling/permutation/bootstrap),
according to $H_{0}$ ;

2 Compute the KBQD test statistic;

3 Repeat lines 1-2 $B$ times;

4 Select the 95th quantile.

Algorithm 1 Computation of Critical Value

Bootstrap and subsampling are two distinct methods for sampling data.
Bootstrap involves drawing data points from the original dataset with
replacement, creating new samples that are of the same size as the original
dataset. On the other hand, subsampling creates new samples without
replacement and typically these samples have smaller sample size, that is
$n_{B}=b*n$ and $b\in(0,1]$. There is no clear guidance for the choice of the
”optimal” subsample size $n_{B}$ and, the literature investigates this aspect
according to optimal subsampling probabilities formulated by minimizing some
function of the asymptotic distribution. For a more detailed discussion see
Politis et al. [2012] and Yao and Wang [2021]. On the other hand, the
permutation algorithm generates new samples of the same sample size drawing
from the observations in the original data set. Notice that this procedure
coincides with the subsampling algorithm with $b=1$, where new samples with
the same sample size of the original data set are generated by sampling
without replacement. Markatou et al. [2023] investigated the algorithm for
computing the critical value with respect to the chosen sampling algorithm and
the number of repetitions $B$. The permutation and bootstrap sampling methods
showed higher and more stable performance, in terms of level and power, for
high dimensions and for low sample size. On the other hand, the subsampling
algorithm is preferred in case of large sample size, since generating samples
with smaller sample sizes requires significantly less computational time. It
is worth noting that the choice of the number of
bootstrap/subsampling/permutation replications $B$ can have an evident
influence in terms of performance. Previous simulation studies indicate that
varying the value of $B$ does not consistently impact the kernel-based
quadratic distance tests, in terms of both level and power. While smaller
values of $B$ may achieve slightly higher power, they can also introduce
instability in the results, potentially leading to different conclusions in
the goodness-of-fit problem. A larger $B$ can significantly increases the
computational time. In our implementation, we suggest $B=150$ since it
provides a balance between satisfactory statistical power while maintaining
result stability and practical computational time. However, note that users
are given the flexibility to adjust the value of $B$ based on the specific
needs and computational resources.

###  3.2 The function kb.test

In this section we describe the usage of the function kb.test. The datasets
used in the rest of the section have been generated in R and then saved.

The kb.test function performs the kernel-based quadratic distance Goodness-of-
Fit tests using the Gaussian kernel with tuning parameter $h$. In particular,
the introduced Goodness-of-Fit tests are performed depending on the given
input y. If y is not given, the function computes the test for normality of
sample x, while if y is a vector indicating the group assignment, the
$k$-sample test is performed. The special case of the two-sample problem can
be additionally assessed providing the two samples to be compared are x and y.
The function takes the following arguments.

  * x

numeric matrix or vector of data values.

  * y

NULL or numeric matrix or vector of data values or vector of group
assignments.

  * h

: tuning parameter.

  * method

: sampling algorithm (”subsampling”,”bootstrap”,”permutation”)

  * B

the number of iterations to use for critical value estimation, B = 150 as
default.

  * b

subsample size in the subsampling algorithm.

  * Quantile

the quantile to use for critical value estimation, 0.95 is the default value.

  * mu_hat

Mean vector for the reference distribution.

  * Sigma_hat

Covariance matrix of the reference distribution.

  * centeringType

String indicating the method used for centering the normal kernel, ”Param” or
”Nonparam”. (default: ”Nonparam”). For $k$-sample tests only ”Nonparam” is
available.

  * alternative

Family of alternatives chosen for selecting $h$, between location, scale and
skewness (only if $h$ is not provided).

Except for the data entry x (and y), and the bandwidth parameter h, the other
arguments do not need to be specified. In such case, the implementation
follows the default setting as specified in the list above. The function
returns a class object with the results and summary of the performed test. In
the case that the value of the tuning parameter $h$ is unknown, the provided
function select_h selects the ”optimal” value of $h$ according to the mid-
power analysis algorithm in [Markatou et al., 2023]. The description and
implementation of this algorithm are treated in section 3.4.

In the following examples, we show the usage of the kb.test function in
different contexts in which the available tests, that is the test for
normality, the two-sample test or the $k$-sample test, are needed. For each
example the obtained results are displayed, together with the needed
computational time. The value of $h$ has been fixed in advance in these
examples. Recall that, in the general setting, we recommend that the optimal
value of $h$ is obtained by the use of the function select_h. See the
simulation study in Markatou et al. [2023] for more details about the choice
of the default setting here.

######  Example 3.1 (Test for normality).

Here, we illustrate the usage of the introduced function for the normality
test. We generate one sample from a multivariate standard normal distribution,
that is $\boldsymbol{x}=(\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})\sim
N_{d}(\boldsymbol{0},\boldsymbol{I}_{d})$ with dimension $d=4$, and sample
size $n=500$. Listing 1 and 2shows the usage in R and Python, including the
provided output.

Listing 1: R Code for the Normality test.

⬇

> library(QuadratiK)

> library(mvtnorm)

> n <- 500

> d <- 4

> set.seed(2468)

> dat_norm <- rmvnorm(n, sigma = diag(d))

> h = 0.4

> set.seed(2468)

> system.time(norm_test <- kb.test(x=dat_norm, h=h, centering="Param"))

user system elapsed

0.57 0.11 0.68

> norm_test

Kernel-based quadratic distance Normality test

H0 is rejected: FALSE

Test Statistic: 2.685621e-06

Critical value (CV): 0.0002145414

CV method: Empirical

Selected tuning parameter h: 0.4

Listing 2: Python code for the Normality test.

⬇

>>> import numpy as np

>>> from QuadratiK.kernel_test import KernelTest

>>> normality_test = KernelTest(h=0.4, centering_type="param",
num_iter=150,random_state=42).test(data)

>>> print("Test : {}".format(normality_test.test_type_))

>>> print("Execution time: {:.3f}".format(normality_test.execution_time))

>>> print("H0 is Rejected : {}".format(normality_test.h0_rejected_))

>>> print("Test Statistic : {}".format(normality_test.test_statistic_))

>>> print("Critical Value (CV) : {}".format(normality_test.cv_))

>>> print("CV Method : {}".format(normality_test.cv_method_))

>>> print("Selected tuning parameter : {}".format(normality_test.h))

Test : Kernel-based quadratic distance Normality test

Execution time: 1.863

H0 is Rejected : False

Test Statistic : 2.6856207169298032e-06

Critical Value (CV) : 0.00016048191978266087

CV Method : Empirical

Selected tuning parameter : 0.4

The output of the test shows the provided value of $h$ as Selected tuning
parameter. Recall that, if such a value is not provided, the algorithm for the
selection of $h$ is automatically performed and the chosen value is displayed.
Additionally, the package provides the summary method for the kb.test output
object. Together with the results of the performed test, this function
provides the qq-plots of each variable with a table of standard descriptive
statistics. Listings 3 and 4 shows the usage of the summary function, while
the generated qq-plots are displayed in Figure 1.

Listing 3: R Code for the summary results of the normality test.

⬇

> summary_norm <- summary(norm_test)

Kernel-based quadratic distance Normality test

Test_Statistic Critical_Value Reject_H0

1 2.685621e-06 0.0002145414 FALSE

Listing 4: Python Code for the summary results of the normality test.

⬇

>>> print(normality_test.summary())

Time taken for execution: 1.863 seconds

Test Results

\-------------- \----------------------------------------------

Test Type Kernel-based quadratic distance Normality test

Test Statistic 2.6856207169298032e-06

Critical Value 0.00016048191978266087

Reject H0 False

\-------------- \----------------------------------------------

Summary Statistics

Feature 0 Feature 1 Feature 2 Feature 3

\------- \----------- \----------- \----------- \-----------

Mean -0.0133 -0.005 0.0523 -0.0264

Std Dev 0.9861 1.0612 0.9795 1.0231

Median -0.0343 0.0223 0.0296 -0.0457

IQR 1.3974 1.4321 1.341 1.3385

Min -2.6755 -3.2222 -2.86 -3.6937

Max 2.5604 3.1177 2.6906 2.7107

Figure 1: Figure automatically generated by the summary function on the result
of the normality test. It displays the normal qq-plots (left) with a table of
the standard descriptive statistics (right) for each variable.

######  Example 3.2 ($k$-sample test).

We generated three samples, with $n=200$ observations each, from a
2-dimensional Gaussian distributions with mean vectors
$\boldsymbol{\mu}_{1}=(0,\frac{\sqrt{3}}{3})$,
$\boldsymbol{\mu}_{2}=(-\frac{1}{2},-\frac{\sqrt{3}}{6})$ and
$\boldsymbol{\mu}_{2}=(\frac{1}{2},-\frac{\sqrt{3}}{6})$, and the Identity
matrix as covariance matrix. In this situation, the generated samples are well
separated, following different Gaussian distributions, i.e.
$\boldsymbol{X}_{1}\sim N_{2}(\boldsymbol{\mu}_{1},\boldsymbol{I})$,
$\boldsymbol{X}_{2}\sim N_{2}(\boldsymbol{\mu}_{2},\boldsymbol{I})$ and
$\boldsymbol{X}_{3}\sim N_{2}(\boldsymbol{\mu}_{3},\boldsymbol{I})$. Listings
5 and 6 illustrates the usage, in R and Python, of the kb.test function for
the $k$-sample problem. It requires the vector y indicating the membership to
groups.

Listing 5: R Code for the k-sample test.

⬇

> sizes <- rep(200,3)

> eps = 1

> set.seed(2468)

> x1 <- rmvnorm(sizes[1], mean = c(0,sqrt(3)*eps/3))

> x2 <- rmvnorm(sizes[2], mean = c(-eps/2,-sqrt(3)*eps/6))

> x3 <- rmvnorm(sizes[3], mean = c(eps/2,-sqrt(3)*eps/6))

> x <- rbind(x1, x2, x3)

> y <- as.factor(rep(c(1,2,3), times=sizes))

> h=1.5

> set.seed(2468)

> system.time(k_test <- kb.test(x=dat_k, y=y_k, h=h))

user system elapsed

0.71 0.08 0.78

> k_test

Kernel-based quadratic distance k-sample test

H0 is rejected: TRUE

Test Statistic: 0.02452253 0.01226127

Critical value (CV): 0.001164279 0.0005821397

CV method: subsampling

Selected tuning parameter h: 1.5

Listing 6: Python code for the k-sample test.

⬇

>>> from QuadratiK.kernel_test import KernelTest

>>> k_sample_test = KernelTest(h = 1.5, method = "permutation",
random_state=42).test(X,y)

>>> print("Test : {}".format(k_sample_test.test_type_))

>>> print("Execution time: {:.3f}
seconds".format(k_sample_test.execution_time))

>>> print("H0 is Rejected : {}".format(k_sample_test.h0_rejected_))

>>> print("Test Statistic : {}".format(k_sample_test.test_statistic_))

>>> print("Critical Value (CV) : {}".format(k_sample_test.cv_))

>>> print("CV Method : {}".format(k_sample_test.cv_method_))

>>> print("Selected tuning parameter : {}".format(k_sample_test.h))

Test : Kernel-based quadratic distance K-sample test

Execution time: 0.400 seconds

H0 is Rejected : True

Test Statistic : [0.02452253 0.01226127]

Critical Value (CV) : [0.00074708 0.00037354]

CV Method : permutation

Selected tuning parameter : 1.5

Recall that the computed test statistics correspond to $(k-1)T_{n}$ and
$T_{n}$, given in equation (7). When the $k$-sample test is performed, the
summary method on the kb.test object returns the results of the tests together
with the standard descriptive statistics for each variable computed, overall
and with respect to the provided groups. The usage of the summary function is
shown in listings 7 and 8.

Listing 7: R Code for the summary results of the $k$-sample test.

⬇

> summary_ktest <- summary(k_test)

Kernel-based quadratic distance k-sample test

Test_Statistic Critical_Value Reject_H0

1 0.02452253 0.0011642795 TRUE

2 0.01226127 0.0005821397 TRUE

> summary_ktest$summary_tables

[[1]]

Group 1 Group 2 Group 3 Overall

mean -0.005959147 -0.5370127 0.5442058 0.0004113282

sd 0.997319811 0.9583059 1.0374834 1.0900980006

median -0.028244038 -0.5477108 0.5297478 -0.0239486027

IQR 1.478884929 1.4105832 1.4234532 1.5377418198

min -2.860006689 -3.1869808 -2.2119189 -3.1869807848

max 2.151784802 2.0647648 3.1580700 3.1580700259

[[2]]

Group 1 Group 2 Group 3 Overall

mean 0.4935364 -0.4042219 -0.2461729 -0.05228613

sd 1.0449582 1.0411639 1.0474989 1.11391575

median 0.5281635 -0.4325995 -0.2950922 -0.09520111

IQR 1.4001089 1.4662111 1.2867345 1.48444495

min -2.6448703 -2.8786352 -3.4932849 -3.49328492

max 3.0792766 2.6788424 2.8290722 3.07927659

Listing 8: Python code for the summary results of the $k$-sample test.

⬇

>>> print(k_sample_test.summary())

Time taken for execution: 0.400 seconds

Test Results

\-------------- \---------------------------------------------

Test Type Kernel-based quadratic distance K-sample test

Test Statistic [0.02452253 0.01226127]

Critical Value [0.00074708 0.00037354]

Reject H0 True

\-------------- \---------------------------------------------

Summary Statistics

Group 1 Group 2 Group 3 Overall

\------------------------ \--------- \--------- \--------- \---------

(’Feature 0’, ’Mean’) -0.006 -0.537 0.5442 0.0004

(’Feature 0’, ’Std Dev’) 0.9973 0.9583 1.0375 1.0901

(’Feature 0’, ’Median’) -0.0282 -0.5477 0.5297 -0.0239

(’Feature 0’, ’IQR’) 1.4789 1.4106 1.4235 1.5377

(’Feature 0’, ’Min’) -2.86 -3.187 -2.2119 -3.187

(’Feature 0’, ’Max’) 2.1518 2.0648 3.1581 3.1581

(’Feature 1’, ’Mean’) 0.4935 -0.4042 -0.2462 -0.0523

(’Feature 1’, ’Std Dev’) 1.045 1.0412 1.0475 1.1139

(’Feature 1’, ’Median’) 0.5282 -0.4326 -0.2951 -0.0952

(’Feature 1’, ’IQR’) 1.4001 1.4662 1.2867 1.4844

(’Feature 1’, ’Min’) -2.6449 -2.8786 -3.4933 -3.4933

(’Feature 1’, ’Max’) 3.0793 2.6788 2.8291 3.0793

######  Example 3.3 (Non-parametric two-sample test).

The two-sample test can be additionally performed by providing the two sample
to be compared as x and y. We generate the sample
$\boldsymbol{y}=(\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{n})$ from a skew-
normal distribution
$SN_{d}(\boldsymbol{0},\boldsymbol{I}_{d},\boldsymbol{\lambda})$, where $d=4$,
$n=200$ and $\boldsymbol{\lambda}=\lambda(1,\ldots,1)$ with $\lambda=0.5$.
Listings 9 and 10 shows the codes in R and Python using the kb.test for
testing if $\boldsymbol{y}$ and $\boldsymbol{x}$, from Example 3.1, follow the
same distribution.

Listing 9: R Code for the two-sample test.

⬇

> library(sn)

> n <- 200

> d <- 4

> skewness_y <- 0.5

> set.seed(2468)

> x_2 <- rmvnorm(n, mean = rep(0,d))

> y_2 <- rmsn(n=n, xi=0, Omega = diag(d), alpha=rep(skewness_y,d))

> h = 2

> set.seed(2468)

> system.time(two_test <- kb.test(x=x_2, y=y_2, h=h))

user system elapsed

0.25 0.00 0.25

> two_test

Kernel-based quadratic distance two-sample test

H0 is rejected: TRUE

Test Statistic: 5.497909e-05

Critical value (CV): 9.957034e-06

CV method: subsampling

Selected tuning parameter h: 2

Listing 10: Python code for the two-sample test.

⬇

>>> from QuadratiK.kernel_test import KernelTest

>>> two_sample_test = KernelTest(h = 2,num_iter = 150,
random_state=42).test(X,Y)

>>> print("Test : {}".format(two_sample_test.test_type_))

>>> print("Execution time: {:.3f}
seconds".format(two_sample_test.execution_time))

>>> print("H0 is Rejected : {}".format(two_sample_test.h0_rejected_))

>>> print("Test Statistic : {}".format(two_sample_test.test_statistic_))

>>> print("Critical Value (CV) : {}".format(two_sample_test.cv_))

>>> print("CV Method : {}".format(two_sample_test.cv_method_))

>>> print("Selected tuning parameter : {}".format(two_sample_test.h))

Test : Kernel-based quadratic distance two-sample test

Execution time: 0.219 seconds

H0 is Rejected : True

Test Statistic : 5.497908846050398e-05

Critical Value (CV) : 7.10139087862447e-06

CV Method : subsampling

Selected tuning parameter : 2

For the two-sample case, the summary function provides the results from the
test and the descriptive statistics per variable and per group, as similarly
described for the $k$-sample test. Additionally, it generates the qq-plots
comparing the quantiles of the two groups for each variable. Listings 11 and
12 show the R and Python code for the summary function, while Figure 2 shows
the generated qq-plots for the example data set.

Listing 11: R Code for the summary results of the two-sample test.

⬇

> summary_two <- summary(two_test)

Kernel-based quadratic distance two-sample test

Test_Statistic Critical_Value Reject_H0

1 5.497909e-05 9.957034e-06 TRUE

> summary_two$summary_tables

[[1]]

Group 1 Group 2 Overall

mean 0.021762263 0.3799990 0.2008806

sd 1.014655344 0.9498167 0.9977884

median -0.005110155 0.3833061 0.2125618

IQR 1.471877262 1.1310211 1.3666010

min -2.675477796 -2.2219439 -2.6754778

max 2.300153117 3.1690406 3.1690406

[[2]]

Group 1 Group 2 Overall

mean -0.03347117 0.2216529 0.09409085

sd 1.06408749 1.0304067 1.05383755

median 0.02476594 0.1717768 0.09272994

IQR 1.52458343 1.3739349 1.45668193

min -3.22222061 -2.6162342 -3.22222061

max 2.96751758 2.3300745 2.96751758

[[3]]

Group 1 Group 2 Overall

mean -0.06473408 0.3312699 0.1332679

sd 0.93818786 0.9868499 0.9818422

median -0.07044427 0.4006745 0.1382735

IQR 1.37135831 1.2185714 1.3854150

min -2.86000669 -3.0246026 -3.0246026

max 2.56476485 2.7590501 2.7590501

[[4]]

Group 1 Group 2 Overall

mean -0.1658894 0.2065222 0.02031639

sd 1.0175325 0.9718613 1.01104987

median -0.2371959 0.1427746 0.04889195

IQR 1.3802070 1.2320445 1.32957715

min -2.5899601 -2.0159679 -2.58996007

max 2.7066430 2.6637589 2.70664302

Listing 12: Python code for the summary results of the two-sample test.

⬇

>>> print(two_sample_test.summary())

Time taken for execution: 0.219 seconds

Test Results

\-------------- \-----------------------------------------------

Test Type Kernel-based quadratic distance two-sample test

Test Statistic 5.497908846050398e-05

Critical Value 7.10139087862447e-06

Reject H0 True

\-------------- \-----------------------------------------------

Summary Statistics

Group 1 Group 2 Overall

\------------------------ \--------- \--------- \---------

(’Feature 0’, ’Mean’) 0.0218 0.38 0.2009

(’Feature 0’, ’Std Dev’) 1.0147 0.9498 0.9978

(’Feature 0’, ’Median’) -0.0051 0.3833 0.2126

(’Feature 0’, ’IQR’) 1.4719 1.131 1.3666

(’Feature 0’, ’Min’) -2.6755 -2.2219 -2.6755

(’Feature 0’, ’Max’) 2.3002 3.169 3.169

(’Feature 1’, ’Mean’) -0.0335 0.2217 0.0941

(’Feature 1’, ’Std Dev’) 1.0641 1.0304 1.0538

(’Feature 1’, ’Median’) 0.0248 0.1718 0.0927

(’Feature 1’, ’IQR’) 1.5246 1.3739 1.4567

(’Feature 1’, ’Min’) -3.2222 -2.6162 -3.2222

(’Feature 1’, ’Max’) 2.9675 2.3301 2.9675

(’Feature 2’, ’Mean’) -0.0647 0.3313 0.1333

(’Feature 2’, ’Std Dev’) 0.9382 0.9868 0.9818

(’Feature 2’, ’Median’) -0.0704 0.4007 0.1383

(’Feature 2’, ’IQR’) 1.3714 1.2186 1.3854

(’Feature 2’, ’Min’) -2.86 -3.0246 -3.0246

(’Feature 2’, ’Max’) 2.5648 2.7591 2.7591

(’Feature 3’, ’Mean’) -0.1659 0.2065 0.0203

(’Feature 3’, ’Std Dev’) 1.0175 0.9719 1.011

(’Feature 3’, ’Median’) -0.2372 0.1428 0.0489

(’Feature 3’, ’IQR’) 1.3802 1.232 1.3296

(’Feature 3’, ’Min’) -2.59 -2.016 -2.59

(’Feature 3’, ’Max’) 2.7066 2.6638 2.7066

Figure 2: Figure automatically generated by the summary function on the result
of the two-sample test. It displays the qq-plots between the two samples
(left) with a table of the standard descriptive statistics for each variable
(right), computed per group and overall.

###  3.3 The function pk.test

The function pk.test performs the kernel-based quadratic distance test for
uniformity on the sphere using the Poisson kernel developed by Ding et al.
[2023]. Using the centered Poisson kernel given in equation (9), the following
two test-statistics are computed:

$T_{n}=\frac{U_{n}}{\sqrt{Var(U_{n})}},\mbox{ with
}Var(U_{n})=\frac{2}{n(n-1)}\left[\frac{1+\rho^{2}}{(1-\rho^{2})^{d-1}}-1\right],$

and

$S_{n}=\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}K_{cen(G)}(\mathbf{x}_{i},\mathbf{x}_{j}).$

The asymptotic distribution of the statistic $S_{n}$ is given by

$S_{n}\sim c\cdot\chi^{2}_{DOF}\mbox{ as }n\to\infty,$

where

$c=\frac{trace(K^{2}_{cen(G)})}{trace(K_{cen(G)})}=\frac{(1+\rho^{2})-(1-\rho^{2})^{d-1}}{(1+\rho)^{d}-(1-\rho^{2})^{d-1}}$

and

$DOF(K_{cen(G)})=\frac{[trace(K_{cen(G)})]^{2}}{trace(K^{2}_{cen(G)})}=\left(\frac{1+\rho}{1-\rho}\right)^{d-1}\left\\{\frac{\left(1+\rho-(1-\rho)^{d-1}\right)^{2}}{1+\rho^{2}-(1-\rho^{2})^{d-1}}\right\\}.$
(10)

Then, the cutoff value of the $S_{n}$ statistic is determined by multiplying
$c$ with the $95^{th}$ quantile of the chi-squared distribution
$\chi^{2}_{\mathrm{DOF}}$. The critical value of statistic $T_{n}$ is obtained
empirically, following Algorithm 1. The function takes the following
arguments.

  * x

numeric matrix of data values.

  * rho

: concentration parameter.

  * B

the number of iterations to use for critical value estimation, B = 300 as
default.

  * Quantile

the quantile to use for critical value estimation, 0.95 is the default value.

######  Example 3.4 (Uniformity test on the Sphere).

We generated $n=200$ observations from the uniform distribution on $S^{d-1}$,
with $d=3$, and $\rho=0.7$. Listings 13 and 14 illustrates the codes in R and
Python using the pk.test for testing uniformity of the generated sample.

Listing 13: R Code for the Uniformity test on the Sphere.

⬇

> library(QuadratiK)

n <- 200

d <- 3

set.seed(2468)

z <- matrix(rnorm(n * d), n, d)

dat_sphere <- z/sqrt(rowSums(z^2))

> rho = 0.7

> set.seed(2468)

> system.time(res_unif <- pk.test(x=dat_sphere, rho=rho))

user system elapsed

0.8 0.0 0.8

> res_unif

Poisson Kernel-based quadratic distance test of Uniformity on the Sphere

Selected consentration parameter rho: 0.7

U-statistic:

H0 is rejected: FALSE

Statistic Un: -0.9756673

Critical value: 0.02660107

V-statistic:

H0 is rejected: FALSE

Statistic Vn: 14.89598

Critical value: 23.22949

Listing 14: Python Code for the Uniformity test on the Sphere.

⬇

>>> from QuadratiK.poisson_kernel_test import PoissonKernelTest

>>> unif_test = PoissonKernelTest(rho = 0.7, random_state=42).test(data)

>>> print("Execution time: {:.3f} seconds".format(unif_test.execution_time))

>>> print("U Statistic Results")

>>> print("H0 is rejected : {}".format(unif_test.u_statistic_h0_))

>>> print("Un Statistic : {}".format(unif_test.u_statistic_un_))

>>> print("Critical Value : {}".format(unif_test.u_statistic_cv_))

>>> print("V Statistic Results")

>>> print("H0 is rejected : {}".format(unif_test.v_statistic_h0_))

>>> print("Vn Statistic : {}".format(unif_test.v_statistic_vn_))

>>> print("Critical Value : {}".format(unif_test.v_statistic_cv_))

Execution time: 0.088 seconds

U Statistic Results

H0 is rejected : False

Un Statistic : -0.9756672826478275

Critical Value : 0.02486113429716727

V Statistic Results

H0 is rejected : False

Vn Statistic : 14.895983386045055

Critical Value : 23.229486935225513

The summary method for the pk.test output object provides the results of the
performed test, and generates a figure showing the qq-plots against the
uniform distribution of each variable with a table of standard descriptive
statistics. Listings 15 and 16 shows the usage of the summary function, while
the generated qq-plots are displayed in Figure 3.

Listing 15: R Code for the Uniformity test on the Sphere.

⬇

> summary_unif <- summary(res_unif)

Poisson Kernel-based quadratic distance test of Uniformity on the Sphere

Test_Statistics Critical_Value Reject_H0

1 -0.9756673 0.02660107 FALSE

2 14.8959834 23.22948694 FALSE

Listing 16: Python Code for the Uniformity test on the Sphere.

⬇

>>> print(unif_test.summary())

Time taken for execution: 0.088 seconds

Test Results

\-------------------------- \-------------------

Test Type Poisson Kernel-based quadratic

distance test of Uniformity on the Sphere

U Statistic Un -0.9756672826478275

U Statistic Critical Value 0.02486113429716727

U Statistic Reject H0 False

V Statistic Vn 14.895983386045055

V Statistic Critical Value 23.229486935225513

V Statistic Reject H0 False

\-------------------------- \-------------------

Summary Statistics

Feature 0 Feature 1 Feature 2

\------- \----------- \----------- \-----------

Mean -0.0196 -0.0289 -0.0093

Std Dev 0.5457 0.6014 0.5867

Median -0.0514 -0.0075 -0.0186

IQR 0.8386 1.1036 1.0335

Min -0.9836 -0.9999 -0.9719

Max 0.9993 0.9722 0.9872

Figure 3: Figure automatically generated by the summary function on the result
of the test for uniformity. It displays the qq-plots between the given samples
and the uniform distribution (left) with a table of the standard descriptive
statistics for each variable (right).

###  3.4 Selection of the tuning parameter

In this section, we present the algorithm for selecting the tuning parameter h
of the normal kernel. This procedure, proposed by Markatou et al. [2023],
properly adapts the strategy based on the mid-power analysis introduced by
Lindsay et al. [2014] to the $k$-sample problem.

The power of the test statistics can be approximated using asymptotic
normality under the alternative through the power function

$\beta(H_{1})=\mathbb{P}(T_{n}>z_{\alpha}\sigma_{0})=\mathbb{P}\left(Z>\frac{z_{\alpha}\sigma_{0}-\mathbb{E}_{H_{1}}(T_{n})}{\sigma_{1}}\right),$

where $z_{\alpha}$ is the $\alpha$-quantile of the normal distribution with
level $\alpha$, $\sigma_{0}^{2}$ and $\sigma_{1}^{2}$ are the variances of the
test statistic $T_{n}$ under the null and the alternative hypotheses,
respectively. We focus on the mid-power manifold $\beta(H_{1})=0.5$, which
allows for a simpler analytic calculation. The algorithm used to select the
tuning parameter $h$ is described in the following steps.

  * •

Identify the family of target alternatives $\\{F_{\delta}\\}$, $\delta\geq 0$.
The case $\delta=0$ corresponds to $H_{0}$, while larger values of $\delta$
indicate more evident departures from the null hypothesis.

  * •

For each value of $h$, define the smallest value of $\delta$ for which the
test achieves power $0.5$ as

$\delta_{mid}(h)=\arg\min\\{\delta:\beta_{h}(\delta)=0.5\\}.$

  * •

The optimal tuning parameter $h^{\ast}$ is defined as the value of $h$ with
the smallest mid-power sensitivity, that is

$h^{\ast}=\arg\min_{h}\\{\delta_{mid}(h)\\}.$

The implementation of the described steps is reported in Algorithm 2.
A fundamental point is the specification of an alternative family of
distributions. Let this be defined by $F_{\delta}(\mu,\Sigma,\lambda)$ where
$\mu,\Sigma,\lambda$ indicate the location, covariance and skewness parameters
of the distribution. Examples of targeted alternatives are provided after the
description of Algorithm 2.

1 Compute the estimates of mean $\hat{\mu}$, covariance matrix $\hat{\Sigma}$
and skewness $\hat{\lambda}$ from the pooled sample.

2 Choose the family of alternatives
$F_{\delta}=F_{\delta}(\hat{\mu},\hat{\Sigma},\hat{\lambda})$.

for  _$\delta$_ do

for  _$h$_ do

3 Generate $\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{K-1}\sim F_{0}$, for
$\delta=0$;

4 Generate $\boldsymbol{X}_{K}\sim F_{\delta}$;

5 Compute the $k$-sample test statistic between
$\boldsymbol{X}_{1},\boldsymbol{X}_{2},\ldots,\boldsymbol{X}_{K}$ with kernel
parameter $h$;

6 Repeat lines 3-5 $N$ times.

7 Compute the power of the test. If it is greater than 0.5, select $h$ as
optimal value.

8If an optimal values has not been selected, choose the $h$ which corresponds
to maximum power.

Algorithm 2 Selection of tuning parameter

Step 2 of the algorithm requires that the family of alternatives
$\\{F_{\delta}\\}$ must be identified. In order to properly handle the data in
hand, the proposed implementation of the algorithm includes the following
target alternatives, with corresponding default values of $\delta$:

  * (i)

location alternatives
$F_{\delta}=SN_{d}(\hat{\mu}+\delta,\hat{\Sigma},\hat{\lambda})$,with
$\delta=0.2,0.3,0.4$;

  * (ii)

scale alternatives
$F_{\delta}=SN_{d}(\hat{\mu},\hat{\Sigma}*\delta,\hat{\lambda})$,
$\delta=0.1,0.3,0.5$;

  * (iii)

skewness alternatives,
$F_{\delta}=SN_{d}(\hat{\mu},\hat{\Sigma},\hat{\lambda}+\delta)$, with
$\delta=0.2,0.3,0.6$.

The following values of $h=0.6,1,1.4,1.8,2.2$ and $N=50$ are set as default.
The algorithm is implemented through the function select_h. The values of
$\delta$ and $h$ can be set by the user, if a more extensive search is
required. Markatou et al. [2023] conducted an extensive simulation study
considering location and skewness alternatives for the two-sample and
$k$-sample problems. Default values here are chosen considering the results
from this simulation study. In case of data following different distributions,
we suggest a more detailed grid search providing more values for $h$, and
$\delta$. This is the recommended practice whenever the computational
resources are available. The function select_h takes as inputs the data matrix
x, vector of labels y, and the type of alternatives (one of ”location”,
”scale” or ”skewness”). Similarly to the function kb.test, the algorithm for
the selection of $h$ for the two-sample test can be performed providing the
two samples as x and y. select_h returns not only the selected value of $h$,
but also the power plot versus the considered list of $h$ values for each
tested value of $\delta$. This provides a set of possible values of $h$ with
high power performance. The following examples show the usage of these
functions considering the data sets previously generated for the two-sample
and k-sample problems.

######  Example 3.5 ($k$-sample test – Continued).

We consider the three 2-dimensional samples generated in Example 3.2, with
$n=200$ observations each. Listings 17 and 18 show the codes in R and Python
for the usage of the function select_h for the $k$-sample tests. This function
needs the input x and y as the function kb.test for the $k$-sample problem.
Figure 4 shows the generated power plot.

Listing 17: R Code for the the selection of $h$.

⬇

> y <- as.numeric(y_k)

> set.seed(2468)

> time_h_k <- system.time( h_k <- select_h(x=dat_k, y=y,

alternative="skewness"))

> time_h_k

user system elapsed

0.90 0.44 93.28

> h_k$h_sel

[1] 1.2

Listing 18: Python Code for the the selection of $h$.

⬇

>>> from QuadratiK.kernel_test import select_h

>>> h_selected, all_powers = select_h(x=X,y=y, alternative="skewness",
power_plot=False, random_state=42,n_jobs = 8)

>>> print("Selected h is: {}".format(h_selected))

Selected h is: 1.2

Figure 4: Plot generated by the select_h function on the result of the
selection of $h$ algorithm for the $k=3$ two-dimensional samples, with size
$n=200$, in Example 3.5. It displays the obtained power versus the considered
$h$, for each value of skewness alternative $\delta$ considered.

As it is seen from Figure 4, when the alternative distribution $F_{\delta}$
with $\delta=0.2$ is considered, there are no values of $h$ within the
indicated range which achieve power greater than or equal to 0.5. For the
second value of $\delta=0.3$, $h=1.6$ is chosen as optimal value since it is
the smallest value with power greater than 0.5. Additionally, it gives a
possible set of values with high power performance.

######  Example 3.6 (Two-sample test – Continued).

We consider the two samples generated in Example 3.3, with sample size $n=200$
and dimension $d=4$. Listings 19 and 20 show the codes in R and Python for the
usage of the function select_h for the two-sample tests. This function needs
the input x and y as the function kb.test for the two-sample problem. Figure 5
shows the generated power plot.

Listing 19: R Code for the the selection of $h$ for the two-sample test.

⬇

> set.seed(2468)

> time_h_2 <- system.time( h_test2 <- select_h(x=x_2, y=y_2,
> alternative="location"))

> time_h_2

user system elapsed

0.67 0.20 17.55

> h_test2$h_sel

[1] 0.8

Listing 20: Python Code for the the selection of $h$ for the two-sample test.

⬇

>>> from QuadratiK.kernel_test import select_h

>>> h_selected, all_powers = select_h(x=X,y=Y, alternative="location",
power_plot=False, random_state=42)

>>> print("Selected h is: {}".format(h_selected))

Selected h is: 2.0

Figure 5: Plot generated by the select_h function on the result of the
selection of $h$ algorithm on the two-sample data set in Example 3.6, with
sample size $n=200$ and dimension $d=4$. It displays the obtained power versus
the considered $h$, for the alternative $\delta=0.2$.

For the two-sample data set, the power analysis with respect to the
alternative distribution with $\delta=0.2$ returns a set of possible values
for $h$ which achieve high power, as showed in Figure 5. Then, different
values of $\delta$ are not explored.

Usually, there are more than one unique value of $h$ that achieve power
greater than or equal to 0.5. In this case, depending on the number of runs
$N$, it is possible for the algorithm to identify one or more values of $h$
that are appropriate to use since they offer power greater than 0.5.

##  4 Clustering on the Sphere

In this section, we consider the problem of clustering data that reside on the
$d$-dimensional sphere. Examples of data that live on such spaces are
directional data and data that can be transformed into directions.
Computational methods for clustering directional data have been proposed in
literature. The non-parametric approaches include algorithms such as K-means
clustering [Ramler, 2008, Maitra and Ramler, 2010], spherical K-means [Dhillon
and Modha, 2001] and online spherical K-means [Zhong, 2005]. In particular,
spherical K-means corresponds to a version of K-means using the cosine
similarity in place of the Euclidean distance. Spherical K-means clustering
can be performed by using the R function skmeans [Hornik et al., 2012] or the
Python package spherecluster.

Banerjee et al. [2005] proposed a clustering method based on a finite mixture
of von Mises–Fisher (vMF) distributions for clustering genomic and text data.
Note that, the spherical K-means algorithm has been shown to be a special case
of a generative model based on a mixture of vMF distributions with equal
concentration parameters and priors for the components [Banerjee and Ghosh,
2002]. Mixture of vMF clustering can be performed using the R function movMF
[Hornik and Grün, 2014] or the Python package spherecluster, and selecting the
approximation for the estimation of the concentration parameters given in
Banerjee et al. [2005].

While most of the clustering methods in literature perform well on small and
medium dimension data sets, high-dimensional data are less frequently
addressed. A promising alternative for dealing with high-dimensional data is
kernel-based methods. In this context, Ng et al. [2001] introduced the
spectral clustering, which relies on the top eigenvectors of an affinity
matrix. Their algorithm is implemented in R using the function specc, within
the kernlab package. Furthermore, Kim et al. [2020] recently proposed
improvements on the spherical $k$-means to overcome the issues related to
document clustering, available in the Python package soyclustering.

As shown in section 3.3, the Poisson kernel density offers a natural way of
assessing goodness of fit for spherical (or spherically transformed) data.
Golzy and Markatou [2020] presented a clustering algorithm based on mixtures
of Poisson kernel-based distributions (PKBD).

###  4.1 Poisson Kernel-Based Density

The $d$-dimensional Poisson kernel in equation 8 is a density function with
respect to the uniform measure on the sphere, written as

$f(\boldsymbol{x}|\rho,\boldsymbol{\mu})=\frac{1-\rho^{2}}{\omega_{d}(1+\rho^{2}-2\rho\boldsymbol{x}\boldsymbol{\mu})^{d/2}},$
(11)

where $\omega_{d}=2\pi^{d/2}\\{\Gamma(\frac{d}{2})\\}^{-1}$ is the surface
area of the unit sphere, $||\boldsymbol{\mu}||=1$ is the vector orienting the
center of the distribution and the concentration parameter $0<\rho<1$ is
related to the variance of the distribution. When $\rho\to 0$, the Poisson
kernel-based density tends to the uniform density on the sphere. Recall that
the parameter $\rho$ is directly connected to the DOF of the corresponding
kernel-based test given in equation (10) [Ding et al., 2023]. In the
univariate case, the Poisson-kernel based density is the circular density also
known as wrapped Cauchy density, while in dimension $d$ it is obtained as the
density of the “exit” distribution on the sphere [Kato and Jones, 2013].

Our QuadratiK package includes the following innovations. First, the function
dpkb evaluates the density function and the function rpkb generates random
observations from the PKB distribution, given the vector $\boldsymbol{\mu}$
and the parameter $\rho$. The function rpkb has the additional input method
with the following options.

  * •

"rejvmf": implementation of the algorithm based on the von Mises-Fisher
distribution using the acceptance-rejection method [Golzy and Markatou, 2020].

  * •

"rejacg": random generator via rejection sampling simulating with angular
central Gaussian envelopes.

  * •

"rejpsaw": random generator simulating with adaptive rejection sampling from
the projected Saw distributions.

The ”rejvmf” sampling scheme is not uniformly bounded and presents a
limitation when $\rho$ approaches $1$. Sablica et al. [2023], in an
interesting paper, introduced two algorithms, implemented here, using
alternative distributions. Notice that, the "rejpsaw" algorithm is built on
the R package Tinflex [Hörmann and Leydold, 2022, Botts et al., 2013]. We
acknowledge the possibility of using the Tinflex package in Python through the
rpy2 package, however, the current implementation of our package in Python
does not provide the "rejpsaw" alternative. We note that rpy2 is an interface
to R running embedded in a Python process. Users, however, will face
additional complexities, such as the manual installation of R and required R
libraries, before utilizing them in Python via rpy2. Due to these additional
complexities and the fact that the Python implementation would cease to be a
standalone implementation and rather become dependent on R, we have not
currently included the “rejpsaw” alternative in the Python implementation.

######  Example 4.1 (Random sampling from PKBD).

We generated observations from the Poisson kernel-based distribution on
$S^{d-1}$, with mean direction $\boldsymbol{\mu}=(0,0,1)$, $\rho=0.85$, and
$d=3$ for illustration purposes. Listing 21 contains the codes in R which
illustrate the usage of the rpkb function for randomly generating observations
from the PKBD. We sampled $n=1000$ observations for each method available.
Figure 6 shows the generated points on the 3-dimensional sphere.

Listing 21: R Code for random sampling from the Poisson kernel-based density
on the Sphere.

⬇

> library(QuadratiK)

> rho = 0.85

> set.seed(2468)

> dat1 <- rpkb(1000, rho=rho, mu=c(0,0,1), method="rejvmf")$x

> dat2 <- rpkb(1000, rho=rho, mu=c(0,0,1), method="rejacg")$x

> dat3 <- rpkb(1000, rho=rho, mu=c(0,0,1), method="rejpsaw")$x

Figure 6: Spherical observations sampled from the PKBD with $\rho=0.85$ using
the three available algorithms, "rejvmf", "rejacg" and "rejpsaw". $n=1000$
points are generated for each method. Figure 7: Computational time needed for
generating $n=1000,10000$ observations from the PKBD using the available
algorithms, with $\rho=0.1,\ldots,0.9,0.95$. The plot on the right shows a
zoom of the performance of the figure on the left, focusing on time less than
10 seconds.

We also compare the available algorithms in terms of computational time. We
generate $n=1000,10000$ observations from the PKBD with $d=3$,
$\boldsymbol{\mu}=(0,0,1)$ and different values of $\rho$, that is
$\rho=0.1,\ldots,0.9,0.95$. Figure 7 displays the computational time needed by
each algorithm for increasing value of concentration parameter $\rho$. Since
the algorithm "rejvmf" is not uniformily bounded, the computational time
needed to generate observations increases when the concentration parameter
$\rho$ is close to 1. The algorithm "rejacg" based on angular central gaussian
distributions performs similarly requiring a small amount of time. The
"rejpsaw" algorithm, while showing competitive computational time, generates
observations with lower variance, more concentrated around the mean direction
as shown in Figure 6. This last algorithm may be more appropriate when values
of $\rho$ close to $1$ are needed.

###  4.2 Clustering Algorithm

The clustering algorithm, presented as Algorithm 3 below, proposed by Golzy
and Markatou [2020] consists of a parametric mixture model approach based on
Poisson kernel-based distributions on the unit sphere.

Given $\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n}$ points on the Sphere and
$M$ number of clusters

0 Initialize $\Theta$.

E –step: for  _k = 1:M_ do

For each observation $\boldsymbol{x}_{i}$, compute the kernel
$K_{\rho_{k}}(\boldsymbol{x}|\boldsymbol{\mu}_{k})$;

Compute the posterior probabilities
$p(k|\boldsymbol{x}_{i},\theta_{k})=\frac{\alpha_{k}f_{k}(\boldsymbol{x}_{i}|\theta_{k})}{f(\boldsymbol{x}_{i}|\Theta)}$;

Compute the weights
$w_{ik}=\frac{p(k|\boldsymbol{x}_{i},\theta_{k})}{1+\rho_{k}^{2}-2\rho_{k}\boldsymbol{x}_{i}\boldsymbol{\mu}_{k}}$;

M–step: for  _k = 1:M_ do

Estimate parameters as:

  * •

$\alpha_{k}=1/n\sum_{i=1}^{n}p(k|\boldsymbol{x}_{i},\theta_{k})$;

  * •

$\mu_{k}=\frac{\sum_{i=1}^{n}w_{ik}\boldsymbol{x}_{i}}{||\sum_{i=1}^{n}w_{ik}\boldsymbol{x}_{i}||}$;

  * •

$\rho_{k}$ is the solution of $g_{k}(y)=0$ in the interval $(0,1)$, where
$g_{k}(y)=\frac{-2ny\alpha_{k}^{(t-1)}}{1-y^{2}}+d\left||\sum_{i=1}^{n}w_{ik}^{(t-1)}\boldsymbol{x}_{i}\right||-dy\sum_{i=1}^{n}w_{ik}^{(t-1)}.$

Perform E–M steps until convergence.

Algorithm 3 Poisson Kernel-Based Clustering on the Sphere

A mixture of $M$ Poisson kernel-based densities given by

$f(\boldsymbol{x}|\Theta)=\sum_{j=1}^{M}\alpha_{j}f_{j}(\boldsymbol{x}|\rho_{j},\boldsymbol{\mu}_{j}),$

where $M$ is the number of clusters, $\alpha_{j}$ are the mixing proportions
with $\sum_{j=1}^{M}\alpha_{j}=1,$ and
$\Theta=\\{\alpha_{1},\ldots,\alpha_{M},\rho_{1},\ldots,\rho_{M},\boldsymbol{\mu}_{1},\ldots,\boldsymbol{\mu}_{M}\\}$,
is the model that is considered.

###  4.3 The function pkbc

In this section we describe the function pkbc for performing clustering on the
sphere. The function takes the following arguments.

  * dat

data matrix of observations to be clustered.

  * nClust

Number of clusters. It can be single value or a numeric vector.

  * maxIter

Maximum number of iterations before a run is terminated. (default: 300)

  * stoppingRule

string describing the stopping rule to be used within each run, one of ’max’,
’membership’, or ’loglik’. (default: ’loglik’)

  * initMethod

string describing the method used for initialization of the centroids.
Currently must be ’sampleData’.

  * numInit

Number of initializations. (default: 10)

To execute the PKBD model based clustering algorithm we need to initialize the
parameters of the model. The sampleData method initializes the centroids using
randomly chosen observations. Then, the final estimates are those with
corresponding highest likelihood. The concentration parameters are set to 0.5
and the algorithm starts with equal mixing proportions. The available stopping
criteria are: ’max’: until the change in the log-likelihood is less than a
given threshold $(1e-7)$; ’membership’: until the membership is unchanged;
’loglik’: based on a maximum number of iterations. The estimation of the
number of clusters is an important issue in clustering. Instead of providing a
single value for nClust, it can be given as a vector of possible number of
clusters. For each value in the vector, the function estimates the parameters
of the PKBD mixture model which are used for assigning memberships as maximum
argument of posterior probabilities.

The function returns a class object pkbc which contains input, a list of input
information, and res_k, a list of the clustering results for the provided
possible number of clusters. Each element of this list is an object with the
following attributes:

  * •

postProbs Matrix of posterior probabilities for each observation and indicated
number of clusters.

  * •

LogLik Maximum value of log-likelihood function.

  * •

wcss Values of within-cluster sum of squares computed with the Euclidean
distance and the cosine similarity.

  * •

params List of estimated paremeters of the mixture model

    * –

mu Estimated centroids;

    * –

rho Estimated concentration parameters;

    * –

alpha Estimated mixing proportions.

  * •

finalMemb Vector of assigned final memberships.

  * •

runInfo List of information of the EM algorithm runs:

    * –

logLikVec Vector of the obtained log-likelihood values;

    * –

numIterPerRun Number of EM iterations per run.

The Quadratik package also provides additional tools which aid the exploration
of the data under investigation. These include methods of cluster validation,
data visualization and descriptive analysis. They are implemented through the
functions validation and summary_stat. In the following section we illustrate
the usage of the presented pkbc function, as well as the supplementary
functions, in a real data application.

###  4.4 The Wireless Indoor Localization Data

We consider the Wireless Indoor Localization Data Set [Bhatt, 2017], publicly
available in the UCI Machine Learning Repository’s website. This data set is
used to study the performance of different indoor localization algorithms
[Rohra et al., 2017, Wang, 2023]. It is also available within the QuadratiK
package as wireless. Listing 22 shows the first rows of the wireless data set.

Listing 22: The wireless data set.

⬇

> library(QuadratiK)

> head(wireless)

V1 V2 V3 V4 V5 V6 V7 V8

1 -64 -56 -61 -66 -71 -82 -81 1

2 -68 -57 -61 -65 -71 -85 -85 1

3 -63 -60 -60 -67 -76 -85 -84 1

4 -61 -60 -68 -62 -77 -90 -80 1

5 -63 -65 -60 -63 -77 -81 -87 1

6 -64 -55 -63 -66 -76 -88 -83 1

The Wireless Indoor Localization data set contains the measurements of the Wi-
Fi signal strength in different indoor rooms. It consists of a data frame with
2000 rows and 8 columns. The first 7 variables report the values of the Wi-Fi
signal strength received from 7 different Wi-Fi routers in an office location
in Pittsburgh (USA). The last column indicates the class labels, from 1 to 4,
indicating the different rooms. Notice that, the Wi-Fi signal strength is
measured in dBm, decibel milliwatts, which is expressed as a negative value
ranging from -100 to 0. In total, we have 500 observations for each room.
Given that the Wi-Fi signal strength takes values in a limited range, it is
appropriate to consider the spherically transformed observations, by $L_{2}$
normalization, and consequently perform the clustering algorithm on the
7-dimensional sphere.

Figure 8 and Figure 9 show a set of plots displayed with respect to the given
labels of the wireless data set and after that data points have been
normalized, respectively. Each figure displays the density plot of individual
variables and pair-wise scatter plots colored by the given labels.
Additionally, the boxes in the upper triangular matrix show the pair-wise
correlations, overall and by group. Notice that the overall structure in the
data set is preserved with normalization.

Listing 23: The pkbc for the wireless data set.

⬇

> wire <- wireless[,-8]

> labels <- wireless[,8]

> wire_norm <- wire/sqrt(rowSums(wire^2))

> set.seed(2468)

> system.time(res_pk <- pkbc(as.matrix(wire_norm),2:10))

user system elapsed

62.20 1.61 64.12

We perform the clustering algorithm on the wireless data set. Listing 23 shows
the R codes which illustrates the usage of the pkbc function. We consider the
$K=2,\ldots,10$ as possible values for the number of clusters.

To guide the choice of the number of clusters, the function validation
provides cluster validation measures and graphical tools. Specifically, it
displays the Elbow plot from the computed within-cluster sum of squared values
and returns an object with IGP, the InGroup Proportion (IGP) [Kapp and
Tibshirani, 2006], and metrics a table of computed evaluation measures. This
include the $k$-sample test statistics with corresponding critical value and,
if the true labels are provided, the measures of adjusted Rand index (ARI),
Macro-Precision and Macro-Recall.

Figure 8: Set of plots of the wireless data displayed with respect to the
given labels. Density plots of individual variables (main diagonal); pair-wise
scatter plots (lower triangular); pair-wise correlations, overall and per
group. Figure 9: Set of plots of the wireless data with spherically
transformed data points, displayed with respect to the given labels. Density
plots of individual variables (main diagonal); pair-wise scatter plots (lower
triangular); pair-wise correlations, overall and per group.

Listing 24 illustrates the usage of the validation function on the clustering
results of the wireless data set, together with the obtained results. Figure
10 shows the generated elbow plots.

Listing 24: The validation function for the wireless data set.

⬇

> set.seed(2468)

> res_validation <- validation(res_pk, true_label = labels)

> res_validation$IGP

[[1]]

NULL

[[2]]

[1] 0.9775475 0.9908515

[[3]]

[1] 0.9706215 0.9860558 0.9429038

[[4]]

[1] 0.9526627 0.9733607 0.9662698 0.9880240

[[5]]

[1] 0.9713701 0.7727273 0.9880240 0.9639831 0.9433498

[[6]]

[1] 0.9544469 0.8664122 0.7226891 0.9880240 0.8608247 0.9848812

[[7]]

[1] 0.8041237 0.7868852 0.9705160 0.8763441 0.9102990 0.9713701 0.9422111

[[8]]

[1] 0.7641509 0.9848812 0.8814229 0.9124579 0.9705160 0.8497409 0.8858696

[8] 0.8041237

[[9]]

[1] 0.7901235 0.9705160 0.8378378 0.7596154 0.9078498 0.7211538 0.8041237

[8] 0.8810811 0.9848812

[[10]]

[1] 0.8866397 0.8640000 0.8279570 0.9274611 0.8227848 0.8922156 0.7474747

[8] 0.8014184 0.9729064 0.9512195

> round(res_validation$metrics, 8)

k=2 k=3 k=4 k=5 k=6

k 2.00000000 3.00000000 4.00000000 5.00000000 6.00000000

Test statistic 0.00000024 0.00000176 0.00000551 0.00000829 0.00001086

Critical value 0.00000000 0.00000001 0.00000001 0.00000011 0.00000010

Reject_H0 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000

ARI 0.31099198 0.69525955 0.94030864 0.91408726 0.85015687

Macro_Precision 0.04760350 0.48150246 0.97718706 0.00328542 0.17222153

Macro_Recall 0.02800000 0.48700000 0.97700000 0.00400000 0.25700000

k=7 k=8 k=9 k=10

k 7.00000000 8.00000000 9.00000000 10.00000000

Test statistic 0.00001419 0.00002141 0.00001785 0.00003005

Critical value 0.00000023 0.00000023 0.00000036 0.00000034

Reject_H0 1.00000000 1.00000000 1.00000000 1.00000000

ARI 0.80291565 0.73174521 0.71924474 0.59713445

Macro_Precision 0.28851541 0.25040392 0.10720486 0.19087003

Macro_Recall 0.10500000 0.26250000 0.12350000 0.27550000

Figure 10: Elbow plots of the within-cluster sum of squares, using the
Euclidean distance (left) and cosine similarity (right), versus the number of
clusters, from 2 to 10, for the wireless data.

The elbow plots and the reported metrics suggest $K=4$ as number of clusters.
This is in accordance with the known ground truth. Once the number of clusters
is selected, the function summary_stat in the QuadratiK package can be used to
obtain additional summary information with respect to the clustering results.
In particular, the function provides mean, standard deviation, median, inter-
quantile range, minimum and maximum computed for each variable, overall and by
the assigned membership, as well as a representation of data points. For $d=2$
and $d=3$, observations are displayed on the unit circle and unit sphere,
respectively. If $d>3$, the spherical PCA is applied on the data set, and the
first 3 principal components are used for visualizing data points on the
sphere, after application of $L_{2}$ normalization. Listing 25 shows the usage
of the summary_stat function on the clustering results of the wireless data
set. Figure 11 shows the data points according to the first three spherical
principal components, colored by the true labels and the assigned membership.

Listing 25: The summary_stat function for the wireless data set.

⬇

> summary_clust <- summary_stat(res_pk,4,true_label=labels)

> summary_clust$metrics

[[1]]

Group 1 Group 2 Group 3 Group 4 Overall

mean -0.29470203 -0.23423419 -0.33976246 -0.34377323 -0.30363137

sd 0.01581506 0.05324085 0.01283014 0.01298448 0.05251967

median -0.29625176 -0.24604104 -0.33952107 -0.34242063 -0.31847692

IQR 0.02182651 0.03014941 0.01674475 0.01835030 0.06412536

min -0.34994496 -0.30643954 -0.37455424 -0.39357081 -0.39357081

max -0.23847076 -0.06308050 -0.29339739 -0.31226867 -0.06308050

[[2]]

Group 1 Group 2 Group 3 Group 4 Overall

mean -0.32636022 -0.35918738 -0.30636324 -0.31540450 -0.32655905

sd 0.01842337 0.02052666 0.01392354 0.01608549 0.02635941

median -0.32487137 -0.35809184 -0.30600005 -0.31538238 -0.32221067

IQR 0.02426854 0.02460137 0.01958541 0.02263003 0.03624468

min -0.40026324 -0.45578394 -0.35600342 -0.36631517 -0.45578394

max -0.27775661 -0.30605235 -0.26903743 -0.27586207 -0.26903743

[[3]]

Group 1 Group 2 Group 3 Group 4 Overall

mean -0.31417147 -0.35786158 -0.32905221 -0.28929074 -0.32231793

sd 0.01710971 0.02459888 0.01563476 0.02094437 0.03163402

median -0.31237754 -0.35575396 -0.32915518 -0.29072716 -0.32209343

IQR 0.02422092 0.03437142 0.01957377 0.02851547 0.03973413

min -0.38676345 -0.43623816 -0.38248214 -0.33639128 -0.43623816

max -0.26460827 -0.28237524 -0.28256746 -0.23485570 -0.23485570

[[4]]

Group 1 Group 2 Group 3 Group 4 Overall

mean -0.30079691 -0.24111635 -0.34918727 -0.35013402 -0.31082318

sd 0.01973759 0.04839750 0.01475930 0.01679759 0.05251032

median -0.30107248 -0.24898566 -0.34858664 -0.35004522 -0.32519736

IQR 0.02672633 0.03356955 0.01898212 0.02286452 0.06907237

min -0.35542814 -0.31690470 -0.40129910 -0.41273961 -0.41273961

max -0.23119131 -0.07399736 -0.30779911 -0.30906423 -0.07399736

[[5]]

Group 1 Group 2 Group 3 Group 4 Overall

mean -0.37583631 -0.43319403 -0.38226045 -0.28257085 -0.36807675

sd 0.01814674 0.02787579 0.01770626 0.02110273 0.05820777

median -0.37549563 -0.43009295 -0.37991233 -0.28337199 -0.37738801

IQR 0.02541214 0.04108696 0.02531565 0.02710382 0.06966093

min -0.42808192 -0.52580164 -0.45385132 -0.34353824 -0.52580164

max -0.33333333 -0.36916034 -0.34027255 -0.21464345 -0.21464345

[[6]]

Group 1 Group 2 Group 3 Group 4 Overall

mean -0.48366352 -0.46327257 -0.45125600 -0.49713921 -0.47391075

sd 0.01814601 0.02195263 0.01446455 0.01412842 0.02488735

median -0.48319020 -0.46519225 -0.45081811 -0.49649227 -0.47411000

IQR 0.02503768 0.02770315 0.01833254 0.01895624 0.03825696

min -0.53130759 -0.53938031 -0.51475369 -0.53840040 -0.53938031

max -0.44107352 -0.40313011 -0.41652096 -0.45551482 -0.40313011

[[7]]

Group 1 Group 2 Group 3 Group 4 Overall

mean -0.48984740 -0.46863418 -0.45728346 -0.49701329 -0.47827795

sd 0.01993870 0.02396734 0.01710219 0.01587742 0.02515570

median -0.49013407 -0.46920446 -0.45703017 -0.49711582 -0.47914650

IQR 0.02784124 0.03050579 0.02311965 0.02030659 0.03643335

min -0.54571977 -0.56835661 -0.50728615 -0.53648350 -0.56835661

max -0.43815927 -0.38995633 -0.41326597 -0.44423198 -0.38995633

The clusters identified with $k=4$ achieve high performance metrics in terms
of ARI, Macro Precision and Macro Recall. Figure 11 also shows that the
identified clusters follows the initial labels.

Figure 11: Data points of the wireless data set displayed according to the
first three principal components. Observations are colored with respect to the
final memberships assigned by the clustering algorithm for $K=4$ clusters
(left) and the true label (right).

##  5 Dashboard application

A graphical user interface (GUI) has been developed for the QuadratiK package,
utilizing the streamlit framework in the Python programming language. The
primary goal of the graphical interface is to facilitate user interaction with
the methods without the need for programming in either R or python. The main
aim of our Dashboard application is to address the needs of practitioners who
wish to analyze findings without extensively engaging with the complexities of
programming.

###  5.1 Instantiating the Dashboard application

The code in Listing 26 is used to start the dashboard application on the host
computer.

Listing 26: Python Code for instantiating the dashboard application.

⬇

>>> from QuadratiK.ui import UI

>>> UI().run()

You can now view your Streamlit app in your browser.

Local URL: http://localhost:8510

Network URL: http://192.168.0.181:8510

###  5.2 Main Menu

Figure 12 shows all the functionalities available on the left hand pane in the
graphical interface. Currently the functionalities available are: ‘Normality
Test’, ‘Two Sample Test’, ‘K-Sample Test’, ‘Tuning Parameter Selection’,
‘Uniformity Test’,‘Data generation from PKBD Models’ and ‘Clustering on
Sphere’.

Figure 12: The first page of the dashboard application displaying the
supported functionalities in the left pane and general information on the
software in the middle.

We show the usage of the functionality of the dashboard through the data sets
generated in the previous examples, providing a detailed explanation of the
steps involved.

### Normality test

Once the user clicks on the “Normality Test” functionality on the left pane,
the interface for performing the “Normality Test” is displayed (Figure 13). A
concise description of the method is included in the beginning, followed by an
expandable component to view the actual Python and R code used in the backend
to generate the results. Presenting the actual code enables users of the GUI
to quickly replicate the results, should they wish to integrate this package
in their pre-existing pipelines.

Figure 13: The user interface for the one-sample test in the dashboard showing
the data upload button and other options available to the user.

The user needs to specify the delimiter for the data files and upload their
datasets in order to perform the test. We allow the users to upload data files
up to a maximum of size 20MBs. Once the data files are uploaded, the user is
presented with a number of input fields where the user can choose to change
the default parameters used for performing the tests. Depending on the
parameters chosen, the test is performed and the results are shown to the
user. We also provide the functionality to the user to download the results.
Furthermore, we present Q-Q plots for each data feature, accompanied by
statistical measures such as the mean, standard deviation, median,
interquartile range (IQR), as well as minimum and maximum values. A sample
output consisting the results of the normality test and the summary statistics
for the uploaded dataset is shown in Figure 14.

Figure 14: The results of the normality test and the summary statistics for
the uploaded dataset.

### Two- and $k$-sample

We implement a similar interface for the ’Two Sample Test’. The user needs to
upload two data files corresponding to the two samples to be tested. Next, we
display the statistical measures for the uploaded datasets and also present
the Q-Q plot for each feature, plotting the two sets of quantiles from the
uploaded datasets.

Figure 15: The user interface for the K-sample test in the dashboard.

In the case of ‘K-Sample test’, shown in Figure 15, the user is required to
upload a file that includes samples along with labels indicating the group to
which each sample belongs. Additionally, the user also needs to specify the
column in which the labels are present, the number of iterations to be used
for critical value estimation, the value of the tuning parameter, and the
proportion of subsampling samples to be used.

Figure 16: The results of the K-sample test performed on the uploaded dataset
alongwith the summary statistics.

The results of the test alongwith summary statistics are shown to the user as
shown in Figure 16. We refrain from presenting the Q-Q plots due to the
escalating complexity as the number of groups of samples increases.

### Tuning Parameter Selection

Figure 17: The user interface for selecting the tuning parameter in the
dashboard.

The interface for Tuning Parameter Selection (Figure 17) is used to compute
the kernel bandwidth of the Gaussian kernel for the proposed kernel-based
quadratic distance (KBQD) tests. We employ a similar interface in which the
user starts by indicating the delimiter, followed by uploading the data file.
Subsequently, the user can specify the column in the data file containing the
label, along with the number of iterations for critical value estimation, the
proportion of subsampling samples to be used, and the alternative, or they can
choose to use the default values. In case of determining $h$ for the normality
test, the user needs to upload the data file with a column containing only a
single label for all observations. Figure 18 shows the results for the tuning
parameter selection algorithm.

Figure 18: The result for the tuning parameter selection algorithm showing the
plot for the $h$ vs power for different $\delta$

### Uniformity Test

To conduct the Uniformity Test, users are again required to indicate the
delimiter utilized in the data file and upload the corresponding file. After
uploading the data file, users have the option to either specify the number of
iterations for critical value estimation and the concentration parameter, or
they can choose to use the default values.

### Data generation from PKBD Models

We also include the functionality for users to generate data from the
supported PKBD algorithms - "rejvmf" and "rejacg".

To generate data from the supported PKBD models users can specify the number
of samples, the concentration parameter value, and the list of location
parameters, as displayed in Figure 19.

If users generate 2-dimensional or 3-dimensional data, we also offer a
visualization of the generated data on a circle or an interactive 3D sphere.
However, if the user generates data for more than 3 dimensions, the
visualization functionality is not supported.

Figure 19: The user interface for generating samples from the supported PKBD
models.

### Clustering on Sphere

Furthermore, we also present an interface for users to use the Poisson Kernel-
Based Clustering algorithm. In order to describe the functionalities and
visualizations incorporated for Poisson Kernel-based clustering, we use the
“Crabs” dataset from R package MASS and perform the clustering on it.

Figure 20: The user interface for the Poisson Kernel-based Clustering
displaying the generated results in a collapsed manner and the corresponding
ARI value Figure 21: Results for the K-Sample Test performed using the cluster
memberships identified by the algorithm

Figure 20 shows the home page of the Dashboard application for the Poisson
Kernel-Based Clustering algorithm. We upload the “Crabs” dataset as a single
data file and explicitly indicate the delimiter used in the file.
Additionally, we also specify whether the true labels for the clusters are
present in the data file and, if so, provide the column number for the same.
Once these parameters are provided, the results are generated. Initially, the
results are displayed in a collapsed manner, requiring the user to click on
them to access the complete generated result (Figure 20). Along with the
results for the clustering, we also show the value for ARI, Macro Precision
and Macro Recall, which for this particular dataset is 0.809, 0.949 and 0.949
respectively. We offer supplementary functionalities such as performing the
K-sample test on the clusters identified by the algorithm, as showed in Figure
21.

Finally, Figure 22 displays additional visualization tools, such as the elbow
plot. The elbow plot is provided to assist the user in selecting the optimal
number of clusters for the given data file. Additionally, a plot of the data
on a three-dimensional interactive sphere is also available to the user.

When visualizing the data on the interactive three-dimensional sphere, if the
true labels are provided, we create two adjacent subplots. In the first
subplot, the data points are colored based on the membership identified by the
PKBD algorithm, and in the next subplot, the data points are colored according
to their true class label. In case the true labels are not present we just
show the first subplot. For the purposes of visualizing the data on the
sphere, if the provided data file has more than three dimensions, PCA is
performed, and only the first three dimensions are chosen for creating the
visualization.

Figure 22: Elbow plot for varied number of clusters.

##  6 Summary

In this article, we presented a new, innovative package that implements
goodness of fit testing and clustering methods that depend on kernel-based
quadratic distance technologies. The package, named QuadratiK, incorporates
many functionalities, such as testing goodness of fit for one, two and
k-sample problems, uniformity testing on the $d$-dimensional sphere and an
algorithm for clustering data on the $d$-dimensional sphere. Furthermore,
algorithms for generating samples from Poisson kernel-based densities, with
given concentration parameters, are implemented. Finally, a Dashboard that
allows users who are not interested in delving into programming details to
incorporate the methods into their workflows is also provided.

The implemented methodology has several advantages, not withstanding the fact
that the methods are mathematically sound and the software is easy to use. Its
implementation in both R and Python programming languages aims to attract a
wider group of users.

## Acknowledgments

We acknowledge help provided by Dr. Yang Chen, Associate Director, GI
Statistics, Takeda Pharmaceutical Company. The second author would like to
acknowledge financial support in the form a research award from KALEIDA Health
Foundation and FDA, that funded the work of the first and third authors.

## References

  * Abbas et al. [2023] S. Abbas, B. Brune, and R. Fried.  robnptests–an r package for robust two-sample location and dispersion tests.  _Journal of Open Source Software_ , 8(82):4947, 2023.
  * Agostinelli and Lund [2023] Claudio Agostinelli and Ulric Lund.  _R package circular: Circular Statistics (version 0.5-0)_, 2023.  URL https://CRAN.R-project.org/package=circular.
  * Ajne [1968] B. Ajne.  A simple test for uniformity of a circular distribution.  _Biometrika_ , 55(2):343–354, 1968.
  * Balogoun et al. [2022] Armando Sosthene Kali Balogoun, Guy Martial Nkiet, and Carlos Ogouyandjou.  A $k$-sample test for functional data based on generalized maximum mean discrepancy, 2022.
  * Banerjee and Ghosh [2002] A. Banerjee and J. Ghosh.  Frequency sensitive competitive learning for clustering on high-dimensional hyperspheres.  In _Proceedings of the 2002 International Joint Conference on Neural Networks. IJCNN’02 (Cat. No.02CH37290)_ , volume 2, pages 1590–1595 vol.2, 2002.  doi: 10.1109/IJCNN.2002.1007755.
  * Banerjee et al. [2005] A. Banerjee, I.S. Dhillon, J. Ghosh, S. Sra, and G. Ridgeway.  Clustering on the unit hypersphere using von mises-fisher distributions.  _Journal of Machine Learning Research_ , 6(9), 2005.
  * Bhatt [2017] Rajen Bhatt.  Wireless Indoor Localization.  UCI Machine Learning Repository, 2017.  DOI: https://doi.org/10.24432/C51880.
  * Bingham [1974] C. Bingham.  An antipodally symmetric distribution on the sphere.  _The Annals of Statistics_ , pages 1201–1225, 1974.
  * Biswas et al. [2014] M. Biswas, M. Mukhopadhyay, and A.K. Ghosh.  A distribution-free two-sample run test applicable to high-dimensional data.  _Biometrika_ , 101(4):913–926, 2014.  ISSN 0006-3444.
  * Bogdan et al. [2002] M Bogdan, K Bogdan, and A Futschik.  A data driven smooth test for circular uniformity.  _Annals of the Institute of Statistical Mathematics_ , 54:29–44, 2002.
  * Botts et al. [2013] Carsten Botts, Wolfgang Hörmann, and Josef Leydold.  Transformed density rejection with inflection points.  _Statistics and Computing_ , 23:251–260, 2013.
  * Cao and Lugosi [2005] R. Cao and G. Lugosi.  Goodness-of-fit tests based on the kernel density estimator.  _Scandinavian Journal of Statistics_ , 32(4):599–616, 2005.  doi: https://doi.org/10.1111/j.1467-9469.2005.00471.x.
  * Chen and Friedman [2017] H. Chen and J.H. Friedman.  A new graph-based two-sample test for multivariate and object data.  _Journal of the American Statistical Association_ , 112(517):397–409, 2017.  ISSN 0162-1459.
  * Chen and Markatou [2020] Y. Chen and M. Markatou.  Kernel tests for one, two, and k-sample goodness-of-fit: State of the art and implementation considerations.  _Statistical Modeling in Biomedical Research: Contemporary Topics and Voices in the Field_ , pages 309–337, 2020.
  * Cuesta-Albertos et al. [2009] Juan A Cuesta-Albertos, Antonio Cuevas, and Ricardo Fraiman.  On projection-based tests for directional and compositional data.  _Statistics and Computing_ , 19:367–380, 2009.
  * Dhillon and Modha [2001] I.S. Dhillon and D.S. Modha.  Concept decompositions for large sparse text data using clustering.  _Machine Learning_ , 42:143–175, 2001.
  * Ding et al. [2023] Y. Ding, M. Markatou, and G. Saraceno.  Poisson kernel-based tests for uniformity on the $d$-dimensional sphere.  _Statistica Sinica_ , 2023.  doi: doi:10.5705/ss.202022.0347.
  * Dortet-Bernadet and Wicker [2008] J.-L. Dortet-Bernadet and N. Wicker.  Model-based clustering on the unit sphere with an illustration using gene expression profiles.  _Biostatistics_ , 9(1):66–80, 2008.
  * Fan [1994] Y. Fan.  Testing the goodness of fit of a parametric density function by kernel method.  _Econometric Theory_ , 10(2):316–356, 1994.  doi: 10.1017/S0266466600008434.
  * Faraway et al. [2019] J. Faraway, G. Marsaglia, J. Marsaglia, and A. Baddeley.  goftest: Classical goodness-of-fit tests for univariate distributions.  _R package version_ , pages 1–2, 2019.
  * Fernández-Durán and Gregorio-Domínguez [2016] Juan José Fernández-Durán and María Mercedes Gregorio-Domínguez.  Circnntsr: An r package for the statistical analysis of circular, multivariate circular, and spherical data using nonnegative trigonometric sums.  _Journal of Statistical Software_ , 70(6):1–19, 2016.  doi: 10.18637/jss.v070.i06.
  * Friedman and Rafsky [1979] J.H. Friedman and L.C. Rafsky.  Multivariate generalizations of the wald-wolfowitz and smirnov two-sample tests.  _The Annals of Statistics_ , 7(4):697–717, 1979\.  ISSN 0090-5364.
  * García-Portugués and Verdebout [2018] Eduardo García-Portugués and Thomas Verdebout.  An overview of uniformity tests on the hypersphere.  _arXiv:1804.00286_ , 2018.  URL https://arxiv.org/abs/1804.00286.
  * García-Portugués and Verdebout [2023] Eduardo García-Portugués and Thomas Verdebout.  _sphunif: Uniformity Tests on the Circle, Sphere, and Hypersphere_ , 2023.  URL https://CRAN.R-project.org/package=sphunif.  R package version 1.1.0.
  * Giné [1975] E. Giné.  Invariant tests for uniformity on compact riemannian manifolds based on sobolev norms.  _The Annals of Statistics_ , 3(6):1243–1266, 1975\.
  * Golzy and Markatou [2020] M. Golzy and M. Markatou.  Poisson kernel-based clustering on the sphere: Convergence properties, identifiability, and a method of sampling.  _Journal of Computational and Graphical Statistics_ , 29(4):758–770, 2020.
  * Golzy et al. [2016] M. Golzy, M. Markatou, and A. Shivram.  Algorithms for clustering on the sphere: Advances & applications.  In _Proceedings of the World Congress on Engineering and Computer Science_ , volume 1, pages 420–425, 2016.
  * Gonzalez-Estrada and Villasenor-Alva [2017] E. Gonzalez-Estrada and J.A. Villasenor-Alva.  goft: Tests of fit for some probability distributions.  _R package version_ , 1(4), 2017.
  * Gretton et al. [2012] A. Gretton, K.M. Borgwardt, M.J. Rasch, B¿ Schölkopf, and A. Smola.  A kernel two-sample test.  _The Journal of Machine Learning Research_ , 13(1):723–773, 2012.
  * Hörmann and Leydold [2022] Wolfgang Hörmann and Josef Leydold.  A generalized transformed density rejection algorithm.  In _Advances in Modeling and Simulation: Festschrift for Pierre L’Ecuyer_ , pages 283–300. Springer, 2022.
  * Hornik and Grün [2014] K. Hornik and B. Grün.  movmf: An r package for fitting mixtures of von mises-fisher distributions.  _Journal of Statistical Software_ , 58(10):1 – 31, 2014.  doi: 10.18637/jss.v058.i10.
  * Hornik et al. [2012] K. Hornik, I. Feinerer, M. Kober, and C. Buchta.  Spherical k-means clustering.  _Journal of Statistical Software_ , 50(10):1–22, 2012.  doi: 10.18637/jss.v050.i10.
  * Hothorn et al. [2008] Torsten Hothorn, Kurt Hornik, Mark A. van de Wiel, and Achim Zeileis.  Implementing a class of permutation tests: The coin package.  _Journal of Statistical Software_ , 28(8):1–23, 2008.  doi: 10.18637/jss.v028.i08.
  * Jupp [2008] P. E. Jupp.  Data-driven Sobolev Tests of Uniformity on Compact Riemannian Manifolds.  _The Annals of Statistics_ , 36(3):1246 – 1260, 2008.  doi: 10.1214/009053607000000541.  URL https://doi.org/10.1214/009053607000000541.
  * Kapp and Tibshirani [2006] Amy V. Kapp and Robert Tibshirani.  Are clusters found in one dataset present in another dataset?  _Biostatistics_ , 8(1):9–31, 04 2006.  ISSN 1465-4644.  doi: 10.1093/biostatistics/kxj029.
  * Kato and Jones [2013] Shogo Kato and M.C. Jones.  An extended family of circular distributions related to wrapped cauchy distributions via brownian motion.  _Bernoulli_ , 19(1):154–171, 2013.
  * Kim et al. [2020] Hyunjoong Kim, Han Kyul Kim, and Sungzoon Cho.  Improving spherical k-means for document clustering: Fast initialization, sparse centroid projection, and efficient cluster labeling.  _Expert Systems with Applications_ , 150:113288, 2020.  ISSN 0957-4174.  doi: https://doi.org/10.1016/j.eswa.2020.113288.  URL https://www.sciencedirect.com/science/article/pii/S0957417420301135.
  * Kruskal and Wallis [1952] William H. Kruskal and W. Allen Wallis.  Use of ranks in one-criterion variance analysis.  _Journal of the American Statistical Association_ , 47(260):583–621, 1952.  doi: 10.1080/01621459.1952.10483441.
  * Lindsay et al. [2008] B.G. Lindsay, M. Markatou, S. Ray, K. Yang, and S.-C. Chen.  Quadratic distances on probabilities: A unified foundation.  _The Annals of Statistics_ , 36(2):983–1006, 2008\.  ISSN 0090-5364.
  * Lindsay et al. [2014] B.G. Lindsay, M. Markatou, and S. Ray.  Kernels, degrees of freedom, and power properties of quadratic distance goodness-of-fit tests.  _Journal of the American Statistical Association_ , 109(505):395–410, 2014.
  * Maitra and Ramler [2010] R. Maitra and I.P. Ramler.  A k-ean-directions algorithm for fast clustering of data on the sphere.  _Journal of Computational and Graphical Statistics_ , 19(2):377–396, 2010.  doi: 10.1198/jcgs.2009.08155.
  * Markatou et al. [2023] M. Markatou, G. Saraceno, and Y. Chen.  Two- and $k$-sample tests based on quadratic distances.  _Manuscript_ , (Department of Biostatistics, University at Buffalo), 2023.
  * Markatou et al. [2017] Marianthi Markatou, Yang Chen, Georgios Afendras, and Bruce G. Lindsay.  _Statistical Distances and Their Role in Robustness_ , pages 3–26.  Springer International Publishing, Cham, 2017.  ISBN 978-3-319-69416-0.  doi: 10.1007/978-3-319-69416-0_1.  URL https://doi.org/10.1007/978-3-319-69416-0_1.
  * Markatou et al. [2021] Marianthi Markatou, Dimitrios Karlis, and Yuxin Ding.  Distance-based statistical inference.  _Annual Review of Statistics and Its Application_ , 8(1):301–327, 2021.  doi: 10.1146/annurev-statistics-031219-041228.
  * Miecznikowski et al. [2013] J. Miecznikowski, A. Vexler, and L. Shepherd.  dbemplikegof: An r package for nonparametric likelihood ratio tests for goodness-of-fit and two-sample comparisons based on sample entropy.  _Journal of Statistical Software_ , 54(3):1–19, 2013.  doi: 10.18637/jss.v054.i03.  URL https://www.jstatsoft.org/index.php/jss/article/view/v054i03.
  * Ng et al. [2001] A. Ng, M. Jordan, and Y. Weiss.  On spectral clustering: Analysis and an algorithm.  In T. Dietterich, S. Becker, and Z. Ghahramani, editors, _Advances in Neural Information Processing Systems_ , volume 14. MIT Press, 2001.  URL https://proceedings.neurips.cc/paper_files/paper/2001/file/801272ee79cfde7fa5960571fee36b9b-Paper.pdf.
  * Panda et al. [2019a] Sambit Panda, Satish Palaniappan, Junhao Xiong, Eric W Bridgeford, Ronak Mehta, Cencheng Shen, and Joshua T Vogelstein.  hyppo: A multivariate hypothesis testing python package.  _arXiv preprint arXiv:1907.02088_ , 2019a.
  * Panda et al. [2019b] Sambit Panda, Cencheng Shen, Ronan Perry, Jelle Zorn, Antoine Lutz, Carey E Priebe, and Joshua T Vogelstein.  Nonpar manova via independence testing.  _arXiv preprint arXiv:1910.08883_ , 2019b.
  * Pewsey and García-Portugués [2021] Arthur Pewsey and Eduardo García-Portugués.  Recent advances in directional statistics.  _Test_ , 30(1):1–58, 2021.
  * Politis et al. [2012] Dimitris N Politis, Joseph P Romano, and Michael Wolf.  _Subsampling_.  Springer Science & Business Media, 2012.
  * Ramler [2008] I.P. Ramler.  _Improved Statistical Methods for k-Means Clustering of Noisy and Directional Data_.  Iowa State University, 2008.
  * Rayleigh [1919] L. Rayleigh.  Xxxi. on the problem of random vibrations, and of random flights in one, two, or three dimensions.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 37(220):321–347, 1919.
  * Rizzo and Székely [2010] Maria L. Rizzo and Gábor J. Székely.  DISCO Analysis: A Nonparametric Extension of Analysis of Variance.  _The Annals of Applied Statistics_ , 4(2):1034 – 1055, 2010.  doi: 10.1214/09-AOAS245.  URL https://doi.org/10.1214/09-AOAS245.
  * Rizzo and Szekely [2022] M.L. Rizzo and G.J. Szekely.  Package ‘energy’.  _R_ , 2022.
  * Rizzo and Székely [2016] M.L. Rizzo and G.J. Székely.  Energy distance.  _WIREs Computational Statistics_ , 8(1):27–38, 2016.  doi: https://doi.org/10.1002/wics.1375.
  * Rohra et al. [2017] Jayant G Rohra, Boominathan Perumal, Swathi Jamjala Narayanan, Priya Thakur, and Rajen B Bhatt.  User localization in an indoor environment using fuzzy hybrid of particle swarm optimization & gravitational search algorithm with neural networks.  In _Proceedings of Sixth International Conference on Soft Computing for Problem Solving: SocProS 2016, Volume 1_ , pages 286–295. Springer, 2017.
  * Rosenbaum [2005] Paul R. Rosenbaum.  An Exact Distribution-Free Test Comparing Two Multivariate Distributions based on Adjacency.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 67(4):515–530, 08 2005.  doi: 10.1111/j.1467-9868.2005.00513.x.
  * Sablica et al. [2023] Lukas Sablica, Kurt Hornik, and Josef Leydold.  Efficient sampling from the PKBD distribution.  _Electronic Journal of Statistics_ , 17(2):2180 – 2209, 2023.  doi: 10.1214/23-EJS2149.  URL https://doi.org/10.1214/23-EJS2149.
  * Scholz and Stephens [1987] F.-W. Scholz and M.A. Stephens.  K-sample anderson–darling tests.  _Journal of the American Statistical Association_ , 82:918–924, 1987.
  * Scholz and Zhu [2019] Fritz Scholz and Angie Zhu.  _kSamples: K-Sample Rank Tests and their Combinations_ , 2019.  URL https://CRAN.R-project.org/package=kSamples.  R package version 1.2-9.
  * Sidak et al. [1999] Zbynek Sidak, Pranab K Sen, and Jaroslav Hajek.  _Theory of Rank Tests_.  Elsevier, 1999.
  * Song [2002] K.-S. Song.  Goodness-of-fit tests based on kullback-leibler discrimination information.  _IEEE Transactions on Information Theory_ , 48(5):1103–1117, 2002.  doi: 10.1109/18.995548.
  * Tsagris et al. [2023] Michail Tsagris, Giorgos Athineou, Christos Adam, Anamul Sajib, Eli Amson, and Micah J. Waldstein.  _Directional: A Collection of Functions for Directional Data Analysis_ , 2023.  URL https://CRAN.R-project.org/package=Directional.  R package version 6.2.
  * Vasicek [1976] O. Vasicek.  A test for normality based on sample entropy.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 38(1):54–59, 1976.  doi: https://doi.org/10.1111/j.2517-6161.1976.tb01566.x.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  Scipy 1.0: Fundamental algorithms for scientific computing in python.  _Nature Methods_ , 17:261–272, 2020.  doi: 10.1038/s41592-019-0686-2.
  * Wald and Wolfowitz [1940] A. Wald and J. Wolfowitz.  On a Test Whether Two Samples are from the Same Population.  _The Annals of Mathematical Statistics_ , 11(2):147 – 162, 1940.  doi: 10.1214/aoms/1177731909.
  * Wang [2023] Xinyue Wang.  The improvement and comparison study of distance metrics for machine learning algorithms for indoor wi-fi localization.  _IEEE Access_ , 11:85513–85524, 2023.  doi: 10.1109/ACCESS.2023.3304334.
  * Yao and Wang [2021] Yaqiong Yao and HaiYing Wang.  A review on optimal subsampling methods for massive datasets.  _Journal of Data Science_ , 19(1):151–172, 2021\.  ISSN 1680-743X.  doi: 10.6339/21-JDS999.
  * Zhong [2005] S. Zhong.  Efficient online spherical k-means clustering.  In _Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005._ , volume 5, pages 3180–3185 vol. 5, 2005.  doi: 10.1109/IJCNN.2005.1556436.