# A likelihood ratio test for circular multimodality

Diego Bolón1<EMAIL_ADDRESS>Rosa M. Crujeiras1
<EMAIL_ADDRESS>Alberto Rodríguez-Casal1<EMAIL_ADDRESS>1Galician Center for Mathematical Research and Technology, CITMAga,
Universidade de Santiago de Compostela, Santiago de Compostela, 15782, Spain

###### Abstract

The modes of a statistical population are high frequency points around which
most of the probability mass is accumulated. For the particular case of
circular densities, we address the problem of testing if, given an observed
sample of a random angle, the underlying circular distribution model is
multimodal. Our work is motivated by the analysis of migration patterns of
birds and the methodological proposal follows a novel approach based on
likelihood ratio ideas, combined with critical bandwidths. Theoretical results
support the behaviour of the test, whereas simulation examples show its finite
sample performance.

###### keywords:

circular data, likelihood ratio test, multimodality, animal orientation

Acknowledgements: This work has been supported by the Spanish Ministerio de
Ciencia, Innovación y Universidades (fellowship FPU20/03960), and by project
PID2020-116587GB-I00 supported by Agencia Estatal de Investigación (AEI) from
Ministerio de Ciencia, Innovación y Universidades. The Supercomputing Centre
of Galicia (CESGA) is acknowledged for providing the computational resources
that allowed to run the simulations in Section 4. We would also like to thank
Jose Ameijeras-Alonso for providing the code of the implementation of the
multimodality test by [5].

##  1 Introduction

Figure 1: Left panel: map of Mexico [30] with the location of the observation
area marked by a red circle. Right panel: location of the new wind farms
erected in the observation area (map retrieved from 24).

A mode of a circular random variable is a point where its density function
reach a local maximum. The existence of multiple modes in a distribution is
intimately related with the presence of different clusters or groups within
the population, being modes high frequency points around which most of the
probability is accumulated. To motivate the need of assessing multimodality in
circular data, that is, to identify whether a (cicular) distribution presents
more than a single mode, we introduce the real data example by [7]. These
authors were interested in determining if the presence of new wind farms had
caused a change on the flight routes of migrating raptors. It is a known fact
that wind turbines have an impact on bird mortality (see, for instance, 20).
But, do migrating birds actively change their migration routes in order to
avoid wind farms? To answer this question, [7] recorded flight trajectories in
six consecutive autumn migration seasons in the Isthmus of Tehuantepec,
Mexico. During that time, two new wind farms were built in the observation
area. Figure 1 shows the location of the observation area within Mexico (letf
panel), and the location of the new wind farms in the observation area (right
panel). These wind farms lie in the middle of the natural route of many
migrating birds that fly from North to South America in order to spend the
boreal winter in warmer latitudes. The authors aimed to check if the new
constructions had altered the raptors’ migration routes and, for that purpose,
they compared the directions of the observed trajectories between two
different periods: 2009-2011, where there was only one wind farm in the
observation zone (pre-construction period), and 2012-2014 where the two new
wind farms were erected (post-construction period).

Figure 2 shows the rose diagrams of observed trajectories during the pre-
construction (left panel) and post-construction (right panel) periods, along
with their mean directions (marked with bullets in the circumferences). From
these simple plots, it is quite clear that there are substantial differences
between both distributions, although it may be argued that the output of a
rose diagram can be influenced by the bin arc-width considered and visual
appreciations are always subjective. [7] applied the Watson and Williams test
(see 29 for details on the procedure) and concluded that there were
significant differences between the two samples. It should be noted that the
Watson and Williams test is devised to assess differences in mean directions,
so any other possible discrepancy is ignored and the conclusion is limited to
a change in the mean, but with no further interpretation on other
characteristics. For the pre-construction observation, which seem to follow a
symmetric and highly concentrated distribution (Figure 2, left), the circular
mean is a suitable summary of the data distribution. However, this is not the
case for the post-construction sample, where observations do not appear to
show a unimodal distribution. Hence, the initial analysis fails to provide a
complete explanation of how the directions of the migration trajectories
changed between the two periods.

Figure 2: Rose diagrams of trajectories directions registered by [7]. Left:
pre-construction period ($n=1183$). Right: post-construction period
($n=1988$). Angle $0$ corresponds with north, $\pi/2$ with east, $\pi$ with
south and $3\pi/2$ with west direction. The mean direction of the sample is
represented with a bullet in the circumferences. Circular mean directions were
approximately $2.475$ rads. for the pre-contruction period and $2.057$ rads.
for the post-contruction period.

Circular kernel density estimators (see equation (2) below) are also a useful
tool for graphical representation of circular data. In Figure 3, three
different kernel density estimators are plotted for the directions of the pre-
construction (left panel) and post-construction (right panel) periods. Figure
3 perfectly illustrates the main flaw of these kinds of estimators: their
highly dependence on the value of the bandwidth $h$. This is particularly
evident for the post-construction data (right panel), where the kernel density
estimators have a very different shape for the three values of $h$ considered.
For $h=0.4$, the estimated density function is unimodal and its mode is
located at $\pi/2$. However, for $h=0.25$, the corresponding estimator has two
modes, one at $\pi/2$ and a secondary one at $\pi$. For the smallest bandwidth
considered ($h=0.1$), a new mode arises between $\pi/2$ and $\pi$.

Figure 3: Linear representations of circular kernel density estimators of
trajectories registered by [7]. Left: pre-construction period ($n=1183$).
Right: post-construction period ($n=1988$). Angle $0$ corresponds with north,
$\pi/2$ with east, $\pi$ with south and $3\pi/2$ with west direction. Each
graph shows the kernel density estimator for three different values of the
bandwidth: $h=0.4$ (blue), $h=0.25$ (green) and $h=0.1$ (red).

[23] use circular kernel density estimates as the basis for a more
sophisticated exploratory tool: the CircSiZer, which allows to determine the
existence of significant modes in a data sample. In Figure 4, the CircSiZers
for each period are shown. In a CircSiZer, as explained in [23], a significant
increasing (decreasing) pattern in the estimated density is shown in blue
(red), so a blue-red pattern indicates the presence of a high frequency point,
that is, a circular mode. Note that, for the pre-construction observations
just the single mode observed in the rose diagram is significant. However, for
the post-construction samples, there is a principal mode (with some
displacement with respect to the mode identified for the pre-construction
sample) but a secondary mode is also noticed, although for a limited range of
smoothing parameters. This will show that most of the raptors fly eastwards,
keeping the wind farms to the south, whereas a small subgroup heads south,
avoiding the wind farms by flying west of them (see Figure 1 for a reference
on where the wind farms are located). In view of the CircSiZer results, we may
wonder if both migration patterns are unimodal and there is just a
displacement trying to avoid the wind farms or if the second group shows a
secondary significant mode, showing the existence of a different migration
strategy.

A formal approach to assess the changes, ensuring the existence of a single
mode in the pre-construction sample and two-modes in the post-construction
observations, is possible through the use of multimodality tests.
Specifically, a multimodality test is a statistical procedure that checks
whether the unknown number of modes is lower or greater than some natural
number. More precisely, if $X$ is a random angle with $j$ modes, a
multimodality test assesses the hypotheses:

$\mathcal{H}_{0}:j\leq k\text{ vs }\mathcal{H}_{1}:j>k;$ (1)

where $k$ is a given natural number. In particular, the simplest (yet useful)
case is the unimodality test:

$\mathcal{H}_{0}:j=1\text{ vs }\mathcal{H}_{1}:j>1.$

Figure 4: CircSiZers for the trajectories directions. Left: pre-construction
period. ($n=1183$) Right: post-construction period ($n=1988$). Angle $0$
corresponds with north, $\pi/2$ with east, $\pi$ with south and $3\pi/2$ with
west direction.

The problem of testing for multimodality for real-valued random variables has
been addressed by different authors. Based on a kernel density estimator, [27]
proposed the use of the critical bandwidth as a test statistic (whose
calibration was corrected by 13), whereas [21] and [8] considered excess-mass
ideas to introduce a testing procedure. [4] combined both views in a procedure
that outperforms the existing proposals. For circular random variables, up to
the authors’ knowledge, the only available tests for multimodality are the
ones proposed by [11], which has been shown to present an unsatisfactory
behaviour by [5], extending these authors the ideas of excess-mass combined
with critical bandwidth to assess multimodality in the circular setting.

In this work, we aim to introduce a multimodality test for circular data that
first, is at least competitive with the existing proposal and outperforms it
when this one fails; and second, the methodology is general enough to be
directly adapted to other settings (e.g. directional, toroidal or cylindrical
data). For that purpose, we will consider the likelihood ratio ideas, with a
nonparametric view, using kernel circular density estimators. All these tools
will be introduced in subsequent sections.

The manuscript is organized as follows. In Section 2, some basic ideas for
circular data modelling are presented, jointly with the existing proposals for
multimodality testing. Section 3 introduces the new multimodality test for
circular data, detailing the calculation of the test statistic and the
calibration method. In Section 4 the behavior of the new test in practice is
analyzed with a simulation study, and its performance is compared with its
competitor. The new proposal is applied to the real data example in Section 5.
Section 6 contains the main conclusions of this paper and some discussion on
possible extensions. The paper is acompained with two appendices with
supplementary material. Appendix A contains the proof of the main result of
this work. For reproducibility purposes, data and code availability are
addressed in Appendix B. Appendix C lists all the models considered in the
simulation study.

##  2 Some background on testing multimodality for circular data

Circular observations are points in the circumference
$\mathbb{S}^{1}=\\{x\in\mathbb{R}^{2}:\|x\|=1\\}$. The statistical analysis of
these types of data has received substantial attention in the literature; see,
for instance, [10], [19], [25], [17] and [26]. A common strategy is to
identify each point of $\mathbb{S}^{1}$ with an angle between $0$ and $2\pi$.
This allows to work with circular data as random angles. A random angle $X$ is
a Borel-measurable mapping from a probability space
$(\Omega,\mathcal{A},\mathbb{P})$ to the interval $(0,2\pi]$:

$X:\Omega\rightarrow(0,2\pi].$

The distribution of an absolutely continuous random angle $X$ is fully
determined by its density function $f$ [19, Sec. 3.2]. This function $f$
verifies:

  * •

$f(x)\geq 0$ almost everywhere.

  * •

$f(x)=f(x+2\pi)$ almost everywhere.

  * •

$\int_{\psi}^{\phi}f(x)dx=\mathbb{P}(\psi<X\leq\phi)$, for all
$\phi,\psi\in(0,2\pi]$ such that $\psi\leq\phi$.

Once the circular density function is determined, the concept of mode
naturally arises. A mode of a random angle $X$ is a point in the interval
$(0,2\pi]$ where its density function reaches a local maximum. Analogously, an
antimode of a random angle is a point in the interval $(0,2\pi]$ where its
density function reaches a local minimum.

Some of the most important families of circular models include the von Mises,
cardioid, wrapped Cauchy and wrappped normal distributions. For a complete
review of distribution circular models, see for instance [19, Sec. 3.5] or
[10, Ch. 3]. For technical reasons, the parametrization for the wrapped normal
distribution used in this paper is slightly different from both references
above. Hereinafter, a wrapped normal distribution
$\mathrm{WN}(\mu,\sigma^{2})$, where $\mu\in(0,2\pi]$ and $\sigma^{2}>0$, is
the circular distribution with density function

$f(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\sum_{m=-\infty}^{+\infty}\exp\left(\frac{-(x-\mu+2\pi
m)^{2}}{2\sigma^{2}}\right),\quad x\in(0,2\pi].$

###  2.1 Multimodality tests for circular data

Since the modes are the local maxima of the density function, the
multimodality testing problem introduced in (1) is closely related with
density function estimation. Most multimodality tests ultimately rely on an
estimator of the density function, like the kernel density estimator (see, for
example, 22). Let $X_{1},\ldots,X_{n}$ be an i.i.d. sample of a random angle
$X$. The circular kernel density estimator of $X$ is given by

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

where $h\in(0,+\infty)$ is a positive real number, and $\mathcal{K}_{h}$ is a
family of symmetric, centered at zero, circular densities indexed by the
parameter $h$. The $\mathcal{K}_{h}$ are called kernel functions, and $h$ is
named bandwidth. A common choice of the kernel function is the wrapped normal
density $\mathrm{WN}(0,h^{2})$:

$\mathcal{K}_{h}(x)=\frac{1}{\sqrt{2\pi
h^{2}}}\sum_{m=-\infty}^{+\infty}\exp\left(\frac{-(x+2\pi
m)^{2}}{2h^{2}}\right).$ (3)

The shape of the kernel estimator $\hat{f}_{h}$ introduced in (2) depends
heavily on the value of the bandwidth $h$. Focusing specifically on the
multimodality problem, if the wrapped normal kernel is chosen, the number of
modes of $\hat{f}_{h}$ is a non decreasing function of $h$, as Figure 3
illustrates (see 15 for a proof). This monotonicity property allows to define
the critical bandwidth for circular data [4].

####  2.1.1 Critical bandwidth test

The critical bandwidth for $k$ modes, namely $h_{k}$, is the smallest value of
$h$ such that the number of modes of $\hat{f}_{h}$ is less or equal to $k$.
That is

$h_{k}=\inf\\{h>0\colon\hat{f}_{h}\text{ has }k\text{ or less modes}\\}.$ (4)

Although the critical bandwidth has been used as a test statistic for
multimodality for real-valued random variables (see 27 and 13), multimodality
tests for circular data do not employ the critical bandwidth directly as a
test statistic. [3] justifies this difference since the properties of the
critical bandwidth are not analogue in both contexts. In the real-valued
setting, [13] show that the null distribution of $h_{k}$ depends on the
curvature of the true density $f$ on its critical points. This makes it
difficult to approximate the null distribution unless one assumes that $f$
only has one critical point. Since circular density always have more than one
critical point (even unimodal densities), the null distribution of $h_{k}$ can
not be directly approximated for circular data.

Trying to overcome this issue, [11] propose a multimodality test for circular
data based on the critical bandwidth, but without using it directly as a test
statistic. Instead, they approach the problem as a goodness-of-fit test,
checking how far is the empirical distribution of the sample from the
distribution given by $\hat{f}_{h_{k}}$. The Watson $U^{2}$ statistic [28], a
circular version of the Cramér-von Mises statistic, is chosen for this
purpose.

Specifically, the proposal by [11] employs the following statistic:

$U^{2}_{n}=n^{-1}\int_{-\infty}^{+\infty}\left(F_{n}(x)-\hat{F}_{h_{k}}(x)-\int_{-\infty}^{+\infty}\left[F_{n}(y)-\hat{F}_{h_{k}}(y)\right]dF_{n}(y)\right)^{2}dF_{n}(x);$
(5)

where $F_{n}$ is the empirical distribution function of the sample, $h_{k}$ is
the critical bandwidth for $k$ modes, and
$\hat{F}_{h_{k}}(x)=\int_{0}^{x}\hat{f}_{h_{k}}(y)dy$. The test rejects the
null hypothesis $\mathcal{H}_{0}:j\leq k$ in (1) for large values of
$U^{2}_{n}$. The null distribution of $U^{2}_{n}$ is approximated via smoothed
bootstrap, resampling from $\hat{f}_{h_{k}}$.

The Fisher and Marron’s test has shown a poor calibration in practice (see 3).
A possible explanation of this phenomenon can be found in the different
converge rates of $h_{k}$ and the optimal bandwidth for distribution
estimation. For real-valued data, [18] showed that, if the underling density
of the sample has $k$ modes, then the critical bandwidth $h_{k}$ is of order
$n^{-1/5}$, whereas the optimal bandwidth for distribution estimation is of
order $n^{-1/3}$ (see 2). Then the $U^{2}_{n}$ statistic from (5) uses a
kernel estimation of the distribution function with a too large bandwidth that
may not correctly fit the data. Consequently, $U^{2}_{n}$ may take large
values even under the null hypothesis.

####  2.1.2 Excess mass test

The other main tool for testing multimodality is excess mass. This concept,
introduced by [21], relies on the intuition that around each mode the density
function $f$ has bumps, that is, regions where $f$ takes values over some
threshold. So, the excess mass tries to measure these bumps by cutting them
from a fixed height.

Given $f$ be a circular density function and $\lambda>0$ a positive real
number, the excess mass of $f$ at level $\lambda$ is

$E(\lambda)=\int_{L(\lambda)}f(x)dx-\lambda\mathrm{Vol}\big{(}L(\lambda)\big{)}=\mathbb{P}\big{(}L(\lambda)\big{)}-\lambda\mathrm{Vol}\big{(}L(\lambda)\big{)},$

where $L(\lambda)=\\{x\in(0,2\pi]:f(x)\geq\lambda\\}$ and
$\mathrm{Vol}\big{(}L(\lambda)\big{)}$ is the Lebesgue measure of the set
$L(\lambda)$. $E(\lambda)$ is measuring the area of the region
$\\{(x,y)\in(0,2\pi]\times\mathbb{R}\colon\lambda<y<f(x)\\}$, as Figure 5
illustrates.

Figure 5:  Linear representation of the excess mass of a circular density
function $f$ at level $\lambda$. $E(\lambda)$ is the area of the gray region.
In this particular case, the set $L(\lambda)$ consists of two intervals,
$[a,b]$ and $[c,d]$, each of them linked to a mode of $f$.

If $f$ is $k$-modal, then $f$ has $k$ bumps, each of them linked to a mode.
This means that the set $L(\lambda)$ is a disjoint union of $k+1$ intervals at
most: one for each mode, and one more if one of the modes is located near
$2\pi$ or $0$ (since, in this case, one of the bumps may extend beyond the
interval $(0,2\pi]$). If this last scenario occurs, two of the $k+1$ disjoint
intervals composing $L(\lambda)$ are of the form $(0,a]$ and $[b,2\pi]$. The
$2\pi$-periodicity of $f$ allows to express the union of these two intervals
as a unique interval of the form $[b,a+2\pi]$, and its probability is

$\mathbb{P}\big{(}[b,a+2\pi]\big{)}=\int_{b}^{a+2\pi}f(x)dx.$

Then, if $f$ at least $k$ modes, one can express excess mass as

$E(\lambda)=\sup_{C_{1},\ldots,C_{k}}\left\\{\sum_{c=1}^{k}\mathbb{P}(C_{c})-\lambda\mathrm{Vol}(C_{c})\right\\};$
(6)

where the supremum is taken over all families of $k$ pairwise disjoint
intervals of $(0,4\pi]$, $\\{C_{1},\ldots,C_{k}\\}$. The left side of (6) can
be estimated empirically, replacing the probability of $f$ by the probability
of the empirical distribution. Then,

$E_{n,k}(\lambda)=\sup_{C_{1},\ldots,C_{k}}\left\\{\sum_{c=1}^{k}\mathbb{P}_{n}(C_{c})-\lambda\mathrm{Vol}(C_{c})\right\\}$

is a natural estimator of the excess mass, where $\mathbb{P}_{n}$ is the
probability of the empirical distribution and, as before, the supremum is
taken over all families of $k$ pairwise disjoint intervals of $(0,4\pi]$,
$\\{C_{1},\ldots,C_{k}\\}$.

However, the quantity $E_{n,k}(\lambda)$ will be a consistent estimator of
$E(\lambda)$ if and only if $f$ has at most $k$ modes. Therefore, the
difference $\big{[}E_{n,k+1}(\lambda)-E_{n,k}(\lambda)\big{]}$ should somehow
measure how far the null hypothesis $\mathcal{H}_{0}:j\leq k$ is from being
true. So the quantity

$\Delta_{n,k+1}=\max_{\lambda>0}\big{[}E_{n,k+1}(\lambda)-E_{n,k}(\lambda)\big{]}$

seems to be a good proposal for testing multimodality of circular data. The
null hypothesis will be rejected for large values of $\Delta_{n,k+1}$.

[5] propose a multimodality test for circular data using this test statistic.
Their idea is to approximate the null distribution of $\Delta_{n,k+1}$
resampling from a slightly modified version of $\hat{f}_{h_{k}}$. This
modification is inspired by the asymptotic properties of $\Delta_{n,k+1}$ in
the real-valued case.

Assuming that the results by [8] also hold for circular data, the null
distribution of $\Delta_{n,k+1}$ depends on the curvature of $f$ in its
critical points. Specifically, if $f$ has $k$ modes and
$x_{1},\ldots,x_{2k}\in(0,2\pi]$ are the critical points of $f$, the
asymptotic distribution of $\Delta_{n,k+1}$ only depends on the values

$d_{l}=\frac{|f^{\prime\prime}(x_{l})|}{{f}^{3}(x_{l})},\quad
l\in\\{1,\ldots,2k\\}.$

[5] suggest to estimate these $d_{l}$ values by kernel estimation methods.
Then, they modify the density function $\hat{f}_{h_{k}}$ so its $d_{l}$ values
are equal to their estimates, and finally, resampling is performed from this
modified version of $\hat{f}_{h_{k}}$. The proposed estimators for the $d_{l}$
values are

$\hat{d}_{l}=\frac{|\hat{f}_{h^{\prime}}^{\prime\prime}(\hat{x}_{l})|}{\hat{f}_{h_{j}}^{3}(\hat{x}_{l})},\quad
l=1,\ldots,2k$

where $\hat{f}_{h_{k}}$ and $\hat{f}^{\prime\prime}_{h^{\prime}}$ are the
kernel estimators of $f$ and $f^{\prime\prime}$, $h_{k}$ is the critical
bandwidth for $k$ modes, $h^{\prime}$ is the bandwidth which minimizes the
asymptotic mean integrated squared error of
$\hat{f}_{h^{\prime}}^{\prime\prime}$ assuming that $f$ is a mixture of von
Mises, and $\hat{x}_{1},\ldots,\hat{x}_{2k}$ are the critical points of
$\hat{f}_{h_{k}}$.

The density $\hat{f}_{h_{k}}$ should be slightly modified to obtain a new
density, $g_{0}$, such that its $d_{l}$ values are equal to their estimators.
This new density $g_{0}$ verifies that $g_{0}(x)=\hat{f}_{h_{k}}(x)$ except on
small neighborhoods around the $2k$ critical points. In those parts small
changes are made so that the position of the modes and antimodes of
$\hat{f}_{h_{k}}$ and $g_{0}$ match, and

$\hat{d_{l}}=\frac{|g_{0}^{\prime\prime}(\hat{x}_{l})|}{g_{0}^{3}(\hat{x}_{l})},\quad\quad\hat{f}_{h_{k}}(\hat{x}_{l})=g_{0}(\hat{x}_{l}),$

for $l=1,\ldots,2k$. Finally, [5] approximate the null distribution of
$\Delta_{n,k+1}$ resampling from this new density $g_{0}$. A detailed
explanation about how the density $\hat{f}_{h_{k}}$ is modified can be found
in [5].

##  3 A new circular multimodality test

In this section, a new multimodality test for circular data is introduced.
This new proposal is based on the likelihood ratio test ideas translated to a
non parametric context. First, a non parametric analogue of the likelihood
function is required, which is achieved via the kernel density estimator.

Let $X_{1},\ldots,X_{n}$ be an i.i.d. sample of a random angle $X$ with
density $f$. The pseudo-likelihood function of the sample is the real function

$\mathcal{L}(h)=\prod_{i=1}^{n}\hat{f}_{h}(X_{i}),$ (7)

where $h>0$, and $\hat{f}_{h}$ is the kernel density estimator defined in (2).
The function $\mathcal{L}$ can be interpreted as a parametric likelihood
function where the assumed parametric density family is indexed by $h>0$, that
is, $\\{\hat{f}_{h}\colon h>0\\}$. However, this family is not independent of
the data and it is directly specified by the sample.

The new multimodality test follows the structure of a parametric likelihood
ratio test. So, the two hypotheses in (1) should be translated into a division
of the parameter space $\Theta=(0,+\infty)$ in order to define the test
statistic. This is performed relying on the kernel density estimator again.
Then, the values of $h$ with a kernel density estimator with $k$ or less modes
are linked with the null hypothesis, and the remaining values to the
alternative hypothesis.

This division of the parameter space can be simplified taking into account the
critical bandwidth $h_{k}$ in (4). By definition, $h_{k}$ is the limit between
kernel density estimators with more and less than $k$ modes. If the bandwidth
$h$ is smaller than $h_{k}$, then $\hat{f}_{h}$ has more than $k$ modes, and
if $h$ is larger than $h_{k}$, then the number of modes of $\hat{f}_{h}$ is at
most $k$. This division completely agrees with the separation between the null
and alternative hypotheses in (1). Then, the parameter space $\Theta$ can be
split by the critical bandwidth as follows

$\displaystyle\Theta_{0}$ $\displaystyle=\\{h>0\colon\hat{f}_{h}\text{ has
}k\text{ or less modes}\\}=[h_{k},+\infty);$ $\displaystyle\Theta_{1}$
$\displaystyle=\\{h>0\colon\hat{f}_{h}\text{ has more than }k\text{
modes}\\}=(0,h_{k}).$

Based on these ideas, the test statistic for addressing (1) is given by:

$D_{k}=2\left[\sup_{h>0}\ell(h)-\sup_{h\geq h_{k}}\ell(h)\right],$ (8)

where $\ell(h)=\log\mathcal{L}(h)$. Similarly to parametric likelihood ratio
tests, large values of $D_{k}$ lead to reject the null hypothesis.

However, the $D_{k}$ statistic introduced in (8) is not well defined. It can
be seen that $\lim_{h\rightarrow 0}\mathcal{L}(h)=+\infty$. Therefore
$\mathcal{L}$ is not upper bounded and $\sup_{h>0}\mathcal{L}(h)$ is not a
real number. One way to solve this issue is to redefine the function
$\mathcal{L}$ given by (7) using cross-validation. Let $\hat{f}^{-i}_{h}$ be
the kernel density estimator in (2) without the $i$-th observation:

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

Then, the cross-validation pseudo-likelihood of the sample
$X_{1},\ldots,X_{n}$ is defined as

$\mathcal{L}_{CV}(h)=\prod_{i=1}^{n}\hat{f}^{-i}_{h}(X_{i}),\quad h>0.$ (9)

The new function $\mathcal{L}_{CV}$ is upper bounded, as Proposition 1
guarantees. The proof of this result can be found in Appendix A.

######  Proposition 1.

Let $X_{1},\ldots,X_{n}$ be an i.i.d. sample of an absolutely continuous
random angle $X$. Let $\mathcal{L}_{CV}$ be the cross-validation pseudo-
likelihood function defined in (9), where $\mathcal{K}_{h}$ is the wrapped
normal density $\mathrm{WN}(0,h^{2})$ in (3). Then, with probability $1$, the
function $\mathcal{L}_{CV}(h)$ is upper bounded in $(0,+\infty)$.

Then, the test statistic

$D_{k}=2\left[\max_{h>0}\ell_{CV}(h)-\max_{h\geq h_{k}}\ell_{CV}(h)\right],$
(10)

where $\ell_{CV}(h)=\log\mathcal{L}_{CV}(h)$, is well defined by Proposition
1. As before, large values of $D_{k}$ reject the null hypothesis. It should be
noted that $D_{k}$ will be null whenever $\max_{h>0}\ell_{CV}(h)=\max_{h\geq
h_{k}}\ell_{CV}(h)$. Consequently, $\mathbb{P}(D_{k}=0)>0$ may hold for some
distributions.

Table 1:  Estimates of $\mathbb{P}(D_{1}=0)$ for Models 1–5. Distribution | $\mathbb{P}_{\mathcal{M}}(D_{1}=0)$
---|---
Model 1 | 0.380
Model 2 | 0.244
Model 3 | 0.331
Model 4 | 0.291
Model 5 | 0.037

For test calibration, one might think of deriving the (asymptotic)
distribution of $D_{k}$ under the null hypothesis. However, our simulation
results indicate that the null distribution of $D_{k}$ is not distribution-
free, that is, the test statistic $D_{k}$ has a different distribution for
each circular population under $\mathcal{H}_{0}:j\leq k$. We simulated
$\mathcal{M}=1000$ samples of size $n=1000$ from five unimodal circular
distributions (Models 1 to 5, see Appendix C for a definition of these
distributions), and the test statistic $D_{1}$ was computed for each sample.
Thus, we obtained $\mathcal{M}=1000$ replicates of $D_{1}$,
$D_{1}^{1},\ldots,D_{1}^{1000}$, that were used to approximate the null
distribution of the test statistic. For example, the probability of
$\mathbb{P}(D_{1}=0)$ was estimated by the proportion of replicates equal to
zero, i.e.:

$\mathbb{P}_{\mathcal{M}}(D_{1}=0)=\frac{1}{\mathcal{M}}\sum_{m=1}^{\mathcal{M}}\mathbb{I}(D_{1}^{m}=0).$

Table 1 shows the estimates of $\mathbb{P}(D_{1}=0)$ for the five
distributions considered. They are significantly different for the five
models. Moreover, in Figure 6, the boxplots of the non-null replicates of
$D_{1}$ are plotted. There exist remarkable differences between the
distribution of $D_{1}$ (conditioned on $D_{1}>0$) under the five models.
These differences are especially noticeable on the right tails, which are
extremely important for test calibration since the null hypothesis is rejected
for large values of $D_{1}$.

Therefore, our simulation results raise some skepticism about the practical
utility of determining the null distribution of $D_{k}$. Instead, a fully
data-driven calibration method is proposed.

Figure 6: Boxplots of $1000$ replicates of the test statistic $D_{1}$
(conditioned on $D_{1}>0$) for five unimodal distributions. Each replicate was
computed from a sample of size $n=1000$. See Appendix C for a definition of
the five models.

###  3.1 Calibration procedure

Smoothed bootstrap is employed for test calibration. Following the ideas from
[5], resamples are simulated from the circular kernel density estimator with
the critical bandwidth, namely $\hat{f}_{h_{k}}$. Replicates of $D_{k}$ are
calculated from these resamples, and the p-value of the test is approximated
by the proportion of replicates larger than the observed value of the
statistic.

Finally, as a summary, an outline of the implementation of the new
multimodality test with significance level $\alpha\in(0,1)$ for circular data
is provided.

  1. 1.

With the original sample of random angles $X_{1},\ldots,X_{n}$, calculate the
critical bandwidth

$h_{k}=\min\\{h>0\colon\hat{f}_{h}\text{ has }k\text{ or less modes}\\}$

and the bandwidths which maximize the cross-validation pseudo-likelihood under
the null and alternative hypothesis:

$\ell_{CV}(h_{\mathrm{max}})=\max_{h>0}\\{\ell_{CV}(h)\\};\quad\ell_{CV}(h_{\mathcal{H}_{0}})=\max_{h\geq
h_{k}}\\{\ell_{CV}(h)\\}.$

Compute the test statistic as:

$D_{k}=2\left[\ell_{CV}(h_{\mathrm{max}})-\ell_{CV}(h_{\mathcal{H}_{0}})\right].$

  2. 2.

The resample $X_{1}^{\ast},\ldots,X_{n}^{\ast}$ is simulated from the smoothed
density $\hat{f}_{h_{k}}$. One replicate of the test statistic, denoted by
$D_{k}^{\ast}$, is calculated from the resample.

  3. 3.

Repeat step 2 $B$ times. Then one will obtain $B$ replicates of the test
statistic: $D_{k}^{\ast,1},D_{k}^{\ast,2},\ldots,D_{k}^{\ast,B}$.

  4. 4.

The null hypothesis ($\mathcal{H}_{0}\colon f$ has $k$ or less modes) is
rejected if the approximated p-value

$\frac{1}{B}\sum_{b=1}^{B}\mathbb{I}(D_{k}^{\ast,b}>D_{k})$

is smaller than $\alpha$.

###  3.2 The $h_{\mathrm{max}}$ bandwidth

The pseudo-likelihood in (9) function can be also used for bandwidth selection
for kernel density estimation. Given the circular sample $X_{1},\ldots,X_{n}$,
one could choose as bandwidth the value $h_{\mathrm{max}}$ verifying:

$\mathcal{L}_{CV}(h_{\mathrm{max}})=\max_{h>0}\mathcal{L}_{CV}(h);$

which exists since $\mathcal{L}_{CV}$ is upper bounded.

The main properties of this bandwidth were analyzed by [14], who proved that
$h_{\mathrm{max}}$ is an asymptotically global optimal bandwidth. This means
that the kernel density estimator $\hat{f}_{h_{\mathrm{max}}}$ is an
asymptotically consistent estimator of the underlying density of the
population $f$, if $f$ is assumed to be bounded away from zero. It is expected
that this optimality property of the $h_{\mathrm{max}}$ bandwidth translates
into a good behavior of the new multimodality test.

Furthermore, the conditions needed for asymptotic convergence of
$\hat{f}_{h_{\mathrm{max}}}$ show a possible flaw of the multimodality test
based on pseudo–likelihood. The density function $f$ is required to be bounded
away from zero for guaranteeing the asymptotic convergence of
$\hat{f}_{h_{\mathrm{max}}}$, so the test may have problems dealing with
densities vanishing at a point or in an arc of positive length. The simulation
study in Section 4 includes some densities of this type to check whether or
not this is a real issue of the test in practice.

##  4 Simulation study

A simulation study was performed to check the calibration and power of the new
test, and to compare its behavior in practice with the existing multimodality
tests for circular data. For this purpose, $\mathcal{M}=1000$ samples of size
$n$ (with $n=100$, $n=500$, and $n=1000$) were simulated from various circular
distributions. The new test with significance level $\alpha$ was applied to
every sample, and the proportion of rejections was calculated. Fifteen
different circular distributions have been considered: five with one mode,
five with two modes, and five with three modes. Details of these fifteen
distributions can be found in Appendix C. The three most common significance
levels were chosen for the simulation study: $\alpha=0.01$, $\alpha=0.05$ and
$\alpha=0.1$. In all cases, $B=500$ resamples were employed for approximating
the p-value.

The simulation study addressed the multimodality test introduced in (1) for
two values of $k$: $k=1$ (unimodality vs. multimodality) and $k=2$ (having at
most two modes vs. having more than two modes). We discuss the obtained
results separately for these two cases.

###  4.1 Testing unimodality

We start regarding the behavior in practice of the new test when testing
unimodality vs. multimodality ($k=1)$. As reported in [3], the test proposed
by [11] suffers from poor calibration. Hence, in our simulations, we only
include the multimodality test by [5] for comparison. The obtained results can
be found in Tables 2 and 3. The new proposal is referred as likelihood test in
both tables.

Table 2:  Proportions of rejections for the test based on pseudo-likelihood (10) and the test by [5] with significance level 1%, 5% and 10% calculated from $\mathcal{M}=1000$ samples. $B=500$ resamples were employed for test calibration. Each row represents a different combination of sampling distribution and sample size. All the considered distributions are unimodal. See Appendix C for a definition of the models. |  | Likelihood test | Ameijeiras et al. test
---|---|---|---
Model | Size | 1% | 5% | 10% | 1% | 5% | 10%
| $n=100$ | 0.001 | 0.013 | 0.037 | 0.007 | 0.036 | 0.067
| $n=500$ | 0.005 | 0.032 | 0.060 | 0.006 | 0.029 | 0.073
M1 | $n=1000$ | 0.003 | 0.027 | 0.058 | 0.003 | 0.039 | 0.072
| $n=100$ | 0.005 | 0.036 | 0.080 | 0.012 | 0.049 | 0.113
| $n=500$ | 0.011 | 0.054 | 0.096 | 0.013 | 0.049 | 0.105
M2 | $n=1000$ | 0.011 | 0.056 | 0.110 | 0.009 | 0.052 | 0.109
| $n=100$ | 0.003 | 0.029 | 0.059 | 0.008 | 0.032 | 0.077
| $n=500$ | 0.013 | 0.046 | 0.087 | 0.004 | 0.036 | 0.079
M3 | $n=1000$ | 0.006 | 0.033 | 0.070 | 0.004 | 0.038 | 0.087
| $n=100$ | 0.011 | 0.048 | 0.075 | 0.003 | 0.031 | 0.075
| $n=500$ | 0.009 | 0.039 | 0.081 | 0.005 | 0.023 | 0.069
M4 | $n=1000$ | 0.006 | 0.034 | 0.057 | 0.005 | 0.035 | 0.072
| $n=100$ | 0.011 | 0.063 | 0.133 | 0.011 | 0.039 | 0.087
| $n=500$ | 0.030 | 0.138 | 0.232 | 0.006 | 0.034 | 0.082
M5 | $n=1000$ | 0.052 | 0.155 | 0.263 | 0.007 | 0.038 | 0.085

Table 2 shows the calibration results. For Models 1 to 4, supported in the
whole circumference, the proportions of rejections of the likelihood test are
similar to or smaller than the significance level for every sample size. The
test is notably conservative in Models 1, 3 and 4, with proportions of
rejections much smaller than the significance levels. Moreover, in these three
cases increasing the sample size from $n=500$ to $n=1000$ results in a
reduction of the proportions of rejections. This phenomenon may seem
surprising at first, but it turns out to be a direct consequence of the nature
of the test.

The type I error of likelihood ratio tests with composite null hypothesis is
not equal for every distribution under the null hypothesis. These types of
tests only achieve a type I error equal to the significance level for
scenarios in the boundary between the null and alternative hypothesis. The new
multimodality test is based on likelihood ratio tests, so one expects to
obtain proportions of rejections much smaller than the significance levels for
distributions that are clearly unimodal, such as Models 1 and 4. On the
contrary, the test should reach proportions of rejections similar to the
significance levels for distributions that are close to be multimodal. This
seems to be the case for Model 2, where the test achieves proportions of
rejections similar to the significance levels for large sample sizes ($n=500$
and $n=1000$).

Finally, the new test shows proportions of rejections larger than the
significance levels for Model 5. This confirms that the test does not perform
as well as it should when the density function vanishes in a circular sector
with positive measure (which is the case for this model), as it was suggested
in Section 3.2.

The test by [5] obtains proportions of rejections similar to or smaller than
the significance levels in every considered scenario. The test exhibits a
conservative behavior in Models 1 and 4, but obtains larger proportions of
rejections than the new proposal in both cases. Finally, for Model 5 the
proportions of rejections are smaller than the significance levels for all
sample sizes. This seems to point out that this test has no problem with
densities vanishing in an arc of positive length, unlike the likelihood ratio
test.

Table 3:  Proportions of rejections for the test based on pseudo-likelihood (10) and the test by [5] with significance level 1%, 5% and 10% calculated from $\mathcal{M}=1000$ samples. $B=500$ resamples were employed for test calibration. Each row represents a different combination of sampling distribution and sample size. All the considered distributions are bimodal. See Appendix C for a definition of the models. |  | Likelihood test | Ameijeiras et al. test
---|---|---|---
Model | Size | 1% | 5% | 10% | 1% | 5% | 10%
| $n=100$ | 0.143 | 0.290 | 0.395 | 0.096 | 0.235 | 0.352
| $n=500$ | 0.423 | 0.720 | 0.846 | 0.486 | 0.711 | 0.823
M6 | $n=1000$ | 0.963 | 0.993 | 0.999 | 0.819 | 0.928 | 0.960
| $n=100$ | 0.034 | 0.083 | 0.129 | 0.020 | 0.091 | 0.160
| $n=500$ | 0.118 | 0.290 | 0.384 | 0.070 | 0.217 | 0.314
M7 | $n=1000$ | 0.231 | 0.422 | 0.542 | 0.118 | 0.292 | 0.398
| $n=100$ | 0.021 | 0.072 | 0.140 | 0.011 | 0.069 | 0.142
| $n=500$ | 0.034 | 0.135 | 0.217 | 0.021 | 0.084 | 0.142
M8 | $n=1000$ | 0.052 | 0.184 | 0.278 | 0.020 | 0.076 | 0.158
| $n=100$ | 0.081 | 0.208 | 0.312 | 0.003 | 0.037 | 0.079
| $n=500$ | 0.317 | 0.513 | 0.616 | 0.006 | 0.041 | 0.079
M9 | $n=1000$ | 0.556 | 0.683 | 0.749 | 0.006 | 0.042 | 0.079
| $n=100$ | 0.154 | 0.379 | 0.530 | 0.006 | 0.032 | 0.080
| $n=500$ | 0.820 | 0.903 | 0.936 | 0.015 | 0.097 | 0.168
M10 | $n=1000$ | 0.957 | 0.980 | 0.989 | 0.071 | 0.260 | 0.391

The power behavior results are shown in Table 3. The new test detects the
alternative hypothesis for all scenarios. Furthermore, the results are
coherent with the shape of the sampling distributions. For example, the power
of the test increases as the distance between the two modes increases (Model 6
vs. Model 7). The power also rises with the excess mass of the secondary mode
(Model 9 vs. Model 10). For Model 8 the proportions of rejections are much
smaller, but always above the significance level. Additionally, the
proportions of rejections grows as the sample size gets larger in all the
scenarios.

Comparing the results of both tests, the new test is clearly more powerful
than the proposal by [5], that only achieves proportions of rejections similar
to the ones obtained by the new test under certain scenarios with small sample
size (Models 7 and 8 with $n=100$). For most cases, the likelihood-ratio test
obtains much larger proportions of rejections than its competitor. This is
specially remarkable for distributions with a secondary mode with small excess
mass, such as Models 9 and 10, since the proposal by [5] relies on excess mass
for testing unimodality. Moreover, this test fails to detect multimodality in
Model 9 due to the extremely small excess mass of the secondary mode. For this
distribution, all the proportions of rejections obtained by Ameijeiras’ test
are smaller than the significance levels.

###  4.2 Testing for at most two modes

The results obtained when testing for at most two modes vs. more than two
modes ($k=2$) are collected in Tables 4 and 5. In Section 4.1, it has been
shown that the new test outperforms the proposal by [4] in terms of power.
Therefore, we only include the results of the new test from now on.

Table 4 contains the calibration results. For the unimodal distributions
(Models 1 to 4), the test obtains proportions of rejections that are notably
smaller than the significance levels for all sample sizes. The behavior of the
new test is very different for Model 5. In this case, the test suffers from
the same problem of calibration reported in Section 4.1, and the proportion of
rejections are larger than the signification levels even for the largest
sample size. The new procedure still has problems with densities that are not
bounded away from zero when testing bimodality.

Table 4:  Proportions of rejections for the test based on pseudo-likelihood (10) with significance level 1%, 5% and 10% when testing bimodality. All proportions of rejections were calculated from $\mathcal{M}=1000$ samples. $B=500$ resamples were employed for test calibration. Each row represents a different combination of sampling distribution and sample size. All the distributions considered had 2 or less modes. See Appendix C for a definition of the models. Model | Size | 1% | 5% | 10% | Model | Size | 1% | 5% | 10%
---|---|---|---|---|---|---|---|---|---
| $n=100$ | 0.002 | 0.014 | 0.038 |  | $n=100$ | 0.001 | 0.051 | 0.138
| $n=500$ | 0.005 | 0.025 | 0.041 |  | $n=500$ | 0.004 | 0.013 | 0.052
M1 | $n=1000$ | 0.004 | 0.019 | 0.043 | M6 | $n=1000$ | 0.006 | 0.020 | 0.047
| $n=100$ | 0.004 | 0.025 | 0.056 |  | $n=100$ | 0.002 | 0.016 | 0.050
| $n=500$ | 0.006 | 0.028 | 0.059 |  | $n=500$ | 0.005 | 0.032 | 0.064
M2 | $n=1000$ | 0.008 | 0.027 | 0.059 | M7 | $n=1000$ | 0.004 | 0.037 | 0.071
| $n=100$ | 0.003 | 0.026 | 0.051 |  | $n=100$ | 0.010 | 0.031 | 0.058
| $n=500$ | 0.005 | 0.038 | 0.074 |  | $n=500$ | 0.008 | 0.046 | 0.092
M3 | $n=1000$ | 0.005 | 0.023 | 0.052 | M8 | $n=1000$ | 0.021 | 0.055 | 0.107
| $n=100$ | 0.002 | 0.015 | 0.041 |  | $n=100$ | 0.007 | 0.026 | 0.064
| $n=500$ | 0.005 | 0.018 | 0.041 |  | $n=500$ | 0.013 | 0.056 | 0.104
M4 | $n=1000$ | 0.004 | 0.021 | 0.053 | M9 | $n=1000$ | 0.011 | 0.059 | 0.105
| $n=100$ | 0.009 | 0.039 | 0.091 |  | $n=100$ | 0.009 | 0.046 | 0.085
| $n=500$ | 0.017 | 0.078 | 0.141 |  | $n=500$ | 0.011 | 0.042 | 0.082
M5 | $n=1000$ | 0.034 | 0.104 | 0.190 | M10 | $n=1000$ | 0.006 | 0.034 | 0.065

For the bimodal distributions (Models 6 to 10), the proportions of rejections
are similar to or smaller than the significance levels in most cases. The only
exception is Model 6, for which the obtained proportion of rejections at
$10\%$ is larger than $0.1$ for the smallest sample size ($n=100$). But this
calibration problem is rapidly corrected increasing the sample size, and the
proportions of rejections are always smaller than the signification levels for
the other two sample sizes.

Comparing the results of the unimodal and bimodal distributions, the
proportions of rejections are always smaller for the unimodal models than the
bimodal ones. Furthermore, the only two cases where the proportions of
rejections are similar to the signification levels for the largest sample size
are bimodal (Models 8 and 9).

Table 5 shows the power results. The new test detects satisfactorily the
alternative hypothesis in the five considered scenarios. For all the models,
the test obtains proportions of rejections that are larger than the
significance level for medium ($n=500$) or large ($n=1000$) sample size, and
the test power increases as $n$ grows. Comparing the results of the five
models leads to some surprising conclusions. The results of the Models 13 to
15 show that the power of the new test naturally decreases as the excess mass
of a secondary mode becomes smaller (Model 13 vs. Model 14) or as the two
secondary modes approach each other (Models 13 vs. Model 15, see Appendix C
for the definition of the fifteen models). Nevertheless, the test obtains
smaller proportions of rejections for Model 11 than for Model 12, which, quite
unexpectedly, suggests that the power of the test grows as the three modes
come closer. This strange behavior may be caused by the particular
distribution of Model 12, that concentrates almost all the probability between
$\pi/2$ and $3\pi/2$. As previously seen with Model 5, this sectorial shape
could lead to inexplicably large proportions of rejections and, therefore, be
the cause of this strange phenomenon.

Table 5:  Proportions of rejections for the test based on pseudo-likelihood (10) with significance level 1%, 5% and 10% when testing bimodality. All proportions of rejections were calculated from $\mathcal{M}=1000$ samples. $B=500$ resamples were employed for test calibration. Each row represents a different combination of sampling distribution and sample size. All the distributions considered had $3$ modes. See Appendix C for a definition of the models. Model | Size | 1% | 5% | 10%
---|---|---|---|---
| $n=100$ | 0.048 | 0.048 | 0.052
| $n=500$ | 0.188 | 0.189 | 0.190
M11 | $n=1000$ | 0.435 | 0.441 | 0.443
| $n=100$ | 0.037 | 0.113 | 0.183
| $n=500$ | 0.209 | 0.395 | 0.496
M12 | $n=1000$ | 0.465 | 0.643 | 0.723
| $n=100$ | 0.588 | 0.715 | 0.776
| $n=500$ | 0.999 | 0.999 | 0.999
M13 | $n=1000$ | 1.000 | 1.000 | 1.000
| $n=100$ | 0.341 | 0.530 | 0.650
| $n=500$ | 0.851 | 0.904 | 0.927
M14 | $n=1000$ | 0.959 | 0.977 | 0.983
| $n=100$ | 0.127 | 0.186 | 0.259
| $n=500$ | 0.432 | 0.584 | 0.694
M15 | $n=1000$ | 0.640 | 0.764 | 0.821

##  5 Real data illustration

Once the new test calibration and power is checked and compared to the
proposal by [5], the performance in practice of both tests is illustrated by
analyzing the migratory bird data [7] described in Section 1. First, the post-
construction observations are regarded in order to assess whether there exists
more than one mode. It should be noticed that, in the simulation study
performed in Section 4, Model 8 was chosen to resemble the distribution of the
post-construction data (see Appendix B for details). For the new test, the
p-value is approximated with $B=1000$ bootstrap resamples, and the estimated
p-value is $0$. For the proposal by [5], the p-value is also estimated with
$B=1000$ resamples, and the estimation is $0$ too. Therefore, both tests
reject the unimodality of the data for the most common significance levels.
This provides statistical evidence that migratory raptors follow at least two
different strategies to avoid the two new wind farms developed in the
observation area, as discussed in Section 1.

Now, the tests are applied to the pre-construction data. In this case, the
estimated p-values are $0.294$ for the new test, and $0.805$ for the test by
[5] (again, $B=1000$ bootstrap resamples were used for approximating the two
p-values). Consequently the two tests confirm the unimodality of the pre-
construction data. Since the post-construction observations have been shown to
be multimodal, this also points out that the construction of the two new wind
farms caused a change in the migratory raptor routes in order to avoid them.

##  6 Conclusions and discussion

This paper introduces a new multimodality test for circular data, based on
likelihood ratio ideas adapted to a nonparametric context. In the exhaustive
simulation study presented in Section 4, the new proposal exhibits a
satisfactory performance for most scenarios, showing a good calibration under
the null hypothesis of unimodality (bimodality) and successfully detecting the
alternative of having more than one mode (two modes). Furthermore, simulation
results show that the new test is consistently more powerful than its main
existing competitor, the proposal by [5], which obtains smaller rejection
proportions in all the considered cases. Finally, the new test is then applied
to a real dataset, illustrating the practical relevance of such procedures.

Perspectives for future research are rich and diverse. First, it could be
interesting to derive the asymptotic distribution of $D_{k}$ defined in
equation (10) under the null hypothesis of $\mathcal{H}_{0}:j\leq k$ in (1).
However, the absence of a closed form expression for $D_{k}$ makes it
challenging to run a theoretical analysis. In addiction, as it was noted in
Section 3, our simulations results raise some skepticism about the utility in
practice of determining the null distribution of $D_{k}$. Anyway, the
asymptotic distribution of $D_{k}$ could help to understand the conservative
behavior of our test under the null hypothesis.

The most promising avenue for future research consists of adapting the testing
procedure to other contexts. The pseudo-likelihood approach only requires a
kernel density estimator $\hat{f}_{h}$ for test implementation, which has
natural extensions to a large variety of spaces, such as the sphere and the
hypersphere (check 19, Ch. 12), the torus [9], or the cylinder [12]. Then, an
analogous reasoning can be followed to test multimodality for data on those
spaces too.

For example, consider the case where $f$ is a density function supported on
the $d$-dimensional hypersphere $\mathbb{S}^{d}$. Given
$\bm{X}_{1},\ldots,\bm{X}_{n}\in\mathbb{S}^{d}$ an i.i.d. sample from $f$, the
kernel density estimator of $f$ (with von Mises-Fisher kernel) is defined as

$\hat{f}_{\kappa}(\bm{x})=a_{d}(\kappa)\frac{1}{n}\sum_{i=1}^{n}\exp(\kappa\bm{X}_{i}^{\prime}\bm{x}),$
(11)

where the quantity $\kappa\in(0,+\infty)$ controls the degree of smoothing and

$a_{d}(\kappa)=\log\bigg{[}\bigg{(}\frac{\kappa}{2}\bigg{)}^{1-d/2}\Gamma\bigg{(}\frac{d}{2}\bigg{)}I_{d/2-1}(\kappa)\bigg{]}.$

Based on the kernel density estimator in (11), we can introduce the pseudo-
likelihood function for hyperspherical data as

$\mathcal{L}_{CV}(\kappa)=\prod_{i=1}^{n}\hat{f}^{-i}_{\kappa}(\bm{X}_{i}),\quad
h>0;$

where

$\hat{f}^{-i}_{\kappa}(\bm{x})=a_{d}(\kappa)\frac{1}{n-1}\sum_{m=1,m\neq
i}^{n}\exp(\kappa\bm{X}_{m}^{\prime}\bm{x}).$

Then, one can address the testing problem in (1) for hyperspherical data using
a new version of the test statistic $D_{k}$ in (10) based on this pseudo-
likelihood function. As in the circular case, large values of $D_{k}$ reject
the null hypothesis.

However, this extension is not as straightforward as it seems at first glance.
For example, the new multimodality test requires that the underling density
function of the data is bounded away from zero for correct calibration,
something that has been already confirmed through simulations. This
requirement indicates that this procedure will only exhibit a good behavior in
compact manifolds (such as the sphere or the torus), since in non-compact
manifolds no density can be bounded away from zero. Moreover, a key fact for
the construction of the test for circular data is that the number of modes of
the kernel density estimator $\hat{f}_{h}$ is monotonous with $h$, a property
that is not preserved in multidimensional spaces. Consequently, the division
of the parameter space between the null and alternative hypotheses is not as
simple as in the circular case, making the computation of $D_{k}$ more
challenging. The lack of monotonicity also impedes defining an analogue of the
critical bandwidth in multidimensional spaces. Then, the current calibration
method, which heavily relies on critical bandwidth, must be reformulated.

##  Appendix A Proof of Proposition 1

Let $\mathcal{K}_{h}$ the density function of a wrapped normal WN$(0,h^{2})$
defined by (3). It is easy to check that
$\lim_{h\rightarrow+\infty}\mathcal{K}_{h}(x)=1/2\pi$ for all
$x\in\mathbb{R}$. So:

$\lim_{h\rightarrow+\infty}\mathcal{L}_{CV}(h)=\lim_{h\rightarrow+\infty}\prod_{i=1}^{n}\hat{f}^{-i}_{h}(X_{i})=(2\pi)^{-n}.$

Split the series present in (3) in half, getting one series for $m\geq 0$ and
other for $m<0$. For the first series, one obtains:

$\sum_{m=0}^{+\infty}\exp\left[-\frac{(x+2\pi
m)^{2}}{2h^{2}}\right]\leq\sum_{m=0}^{+\infty}\exp\left(-\frac{x^{2}}{2h^{2}}-\frac{2\pi^{2}m^{2}}{h^{2}}\right)\leq\\\
\leq\exp\left(-\frac{x^{2}}{2h^{2}}\right)\sum_{m=0}^{+\infty}\exp\left(-\frac{2\pi^{2}m}{h^{2}}\right)=\exp\left(-\frac{x^{2}}{2h^{2}}\right)\left[1-\exp\left(-\frac{2\pi^{2}}{h^{2}}\right)\right]^{-1}.$
(12)

The first inequality of (12) is given by $(a+b)^{2}\geq a^{2}+b^{2}$ for all
$a,b\geq 0$. The last inequality is straightforward, since one is adding more
terms to the series (there are more natural numbers than perfect squares).

A similar reasoning shows that

$\sum_{m=1}^{+\infty}\exp\left[-\frac{(x-2\pi
m)^{2}}{2h^{2}}\right]\leq\exp\left[-\frac{(2\pi-x)^{2}}{2h^{2}}\right]\left[1-\exp\left(-\frac{2\pi^{2}}{h^{2}}\right)\right]^{-1}.$
(13)

With (3), (12) and (13), one deduces

$\mathcal{K}_{h}(x)\leq\left[1-\exp\left(-\frac{2\pi^{2}}{h^{2}}\right)\right]^{-1}\left\\{\frac{1}{\sqrt{2\pi
h^{2}}}\exp\left[-\frac{(2\pi-x)^{2}}{2h^{2}}\right]+\frac{1}{\sqrt{2\pi
h^{2}}}\exp\left(-\frac{x^{2}}{2h^{2}}\right)\right\\};$

so $\lim_{h\rightarrow 0}\mathcal{K}_{h}(x)=0$ for every $x\in(0,2\pi)$.

Assume that the sample $X_{1},\ldots,X_{n}$ has no repeated observations, that
is, $X_{i}\neq X_{m}$ for all $i\neq m$. Then,

$\lim_{h\rightarrow 0}\hat{f}^{-i}_{h}(X_{i})=\lim_{h\rightarrow
0}\frac{1}{n-1}\sum_{m=1,m\neq
i}^{n}\mathcal{K}_{h}\left(X_{i}-X_{m}\right)=0,$

for all $i\in\\{1,\ldots,n\\}$, so $\lim_{h\rightarrow
0}\mathcal{L}_{CV}(h)=0$.

Summarizing, if there are no repeated observations in $X_{1},\ldots,X_{n}$, it
has been shown that

$\lim_{h\rightarrow
0}\mathcal{L}_{CV}(h)=0,\quad\lim_{h\rightarrow+\infty}\mathcal{L}_{CV}(h)=(2\pi)^{-n}.$
(14)

$\mathcal{L}_{CV}(h)$ is a continuous function in $(0,+\infty)$. Therefore,
(14) implies that $\mathcal{L}_{CV}(h)$ is bounded in $(0,+\infty)$. Since $X$
is an absolutely continuous angle, this happens with probability 1.

##  Appendix B Data and code

The data analyzed in this paper were collected by [6] and are available in
Figshare (under license CC BY 4.0,
https://creativecommons.org/licenses/by/4.0/) with the identifier:
https://doi.org/10.6084/m9.figshare.3123100.v2.

The code for reproducing the simulation study can be found in
https://github.com/diboro/testingnmodes. Once the manuscript will be revised,
the final code will be delivered through R package NPCirc, devoted to
nonparametric methods for circular data.

##  Appendix C Simulation models

Table 6 contains the distributions of all the models that had been considered
in the simulation study presented in Section 4. In that table,
$\text{vM}(\mu,\kappa)$ denotes the von Mises distribution with mean direction
$\mu$ and concentration $\kappa$ as in Section 3.5.4 of [19], and
$\mathrm{ssvM}(\mu,\kappa,\lambda)$ is the sine-skewed von Mises with skewing
parameter $\lambda\in[-1,1]$ as defined by [1].

Table 6:  The fifteen models considered in the simulation study of Section 4. All models are mixtures of von Mises except Models 4 and 5. Model | Distribution
---|---
Unimodal:
M1 | vM$(\pi,1)$
M2 | $0.2\cdot\text{vM}(2\pi/3,3)+0.6\cdot\text{vM}(\pi,1.4)+0.2\cdot\text{vM}(4\pi/3,3)$
M3 | $0.05\cdot\text{vM}(2\pi/3,7)+0.9\cdot\text{vM}(\pi,1)+0.05\cdot\text{vM}(4\pi/3,7)$
M4 | ssvM$(\pi,1,-0.9)$
M5111Model 5 is a standard beta(3,2) distribution (as parameterized by 16, Ch. 25) rescaled so that its support is the interval $[\pi/2,3\pi/2]$. | $\pi\cdot\text{beta}(3,2)+\pi/2$
Bimodal:
M6 | $0.5\cdot\text{vM}(\pi-1.25,1.5)+0.5\cdot\text{vM}(\pi+1.25,1.5)$
M7 | $0.5\cdot\text{vM}(\pi-1,1.5)+0.5\cdot\text{vM}(\pi+1,1.5)$
M8 | $0.5\cdot\text{vM}(1.5,4)+0.5\cdot\text{vM}(3,2)$
M9 | $0.95\cdot\text{vM}(\pi/2,6)+0.05\cdot\text{vM}(3\pi/2,3)$
M10 | $0.9\cdot\text{vM}(\pi/2,6)+0.1\cdot\text{vM}(3\pi/2,3)$
Trimodal:
M11 | $1/3\cdot\text{vM}(\pi-2,7)+1/3\cdot\text{vM}(\pi,7)+1/3\cdot\text{vM}(\pi+2,7)$
M12 | $1/3\cdot\text{vM}(\pi-1,7)+1/3\cdot\text{vM}(\pi,7)+1/3\cdot\text{vM}(\pi+1,7)$
M13 | $0.2\cdot\text{vM}(\pi/2,6)+0.2\cdot\text{vM}(\pi,6)+0.6\cdot\text{vM}(7\pi/4,8)$
M14 | $0.1\cdot\text{vM}(\pi/2,6)+0.25\cdot\text{vM}(\pi,6)+0.65\cdot\text{vM}(7\pi/4,8)$
M15 | $0.2\cdot\text{vM}(\pi/2,6)+0.2\cdot\text{vM}(6\pi/7,6)+0.6\cdot\text{vM}(7\pi/4,8)$

Figures 7 and 8 illustrate the density functions of all the fifteen models
considered in the simulation study. They are represented on the real line
(Figure 7) and on the circumference (Figure 8).

Model 1.

Model 6.

Model 11.

Model 2.

Model 7.

Model 12.

Model 3.

Model 8.

Model 13.

Model 4.

Model 9.

Model 14.

Model 5.

Model 10.

Model 15.

Figure 7:  Density functions of the models considered for the simulation
study. Left column: unimodal models. Central column: bimodal models. Right
column: trimodal models.

Model 1.

Model 6.

Model 11.

Model 2.

Model 7.

Model 12.

Model 3.

Model 8.

Model 13.

Model 4.

Model 9.

Model 14.

Model 5.

Model 10.

Model 15.

Figure 8:  Density functions of the models considered for the simulation
study. Left column: unimodal models. Central column: bimodal models. Right
column: trimodal models.

## References

  *   * Abe and Pewsey [2009] Abe T, Pewsey A (2009) Sine-skewed circular distributions. Statistical Papers 52(3):683–707. 10.1007/s00362-009-0277-x
  * Altman and Léger [1995] Altman N, Léger C (1995) Bandwidth selection for kernel distribution function estimation. Journal of Statistical Planning and Inference 46(2):195–214. 10.1016/0378-3758(94)00102-2
  * Ameijeiras-Alonso [2017] Ameijeiras-Alonso J (2017) Assessing simplifying hypotheses in density estimation. PhD thesis, Universidade de Santiago de Compostela
  * Ameijeiras-Alonso et al [2018] Ameijeiras-Alonso J, Crujeiras RM, Rodríguez-Casal A (2018) Mode testing, critical bandwidth and excess mass. TEST 28(3):900–919. 10.1007/s11749-018-0611-5
  * Ameijeiras-Alonso et al [2019] Ameijeiras-Alonso J, Benali A, Crujeiras RM, et al (2019) Fire seasonality identification with multimodality tests. The Annals of Applied Statistics 13(4):2120–2139. 10.1214/19-aoas1273
  * Cabrera-Cruz and Villegas-Patraca [2016a] Cabrera-Cruz SA, Villegas-Patraca R (2016a) Data from: Response of migrating raptors to an increasing number of wind farms. 10.6084/M9.FIGSHARE.3123100.V2
  * Cabrera-Cruz and Villegas-Patraca [2016b] Cabrera-Cruz SA, Villegas-Patraca R (2016b) Response of migrating raptors to an increasing number of wind farms. Journal of Applied Ecology 53(6):1667–1675. 10.1111/1365-2664.12673
  * Cheng and Hall [1998] Cheng MY, Hall P (1998) Calibrating the excess mass and dip tests of modality. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 60(3):579–589. 10.1111/1467-9868.00141
  * Di Marzio et al [2011] Di Marzio M, Panzera A, Taylor CC (2011) Kernel density estimation on the torus. Journal of Statistical Planning and Inference 141(6):2156–2173. 10.1016/j.jspi.2011.01.002
  * Fisher [1993] Fisher NI (1993) Statistical Analysis of Circular Data. Cambridge University Press, 10.1017/cbo9780511564345
  * Fisher and Marron [2001] Fisher NI, Marron JS (2001) Mode testing via the excess mass estimate. Biometrika 88(2):499–517. 10.1093/biomet/88.2.499
  * García-Portugués et al [2013] García-Portugués E, Crujeiras RM, González-Manteiga W (2013) Kernel density estimation for directional–linear data. Journal of Multivariate Analysis 121:152–175. 10.1016/j.jmva.2013.06.009
  * Hall and York [2001] Hall P, York M (2001) On the calibration of Silverman’s test for multimodality. Statistica Sinica 11(2):515–536
  * Hall et al [1987] Hall P, Watson G, Cabrera J (1987) Kernel density estimation with spherical data. Biometrika 74(4):751–762. https://doi.org/10.1093/biomet/74.4.751
  * Huckemann et al [2016] Huckemann S, Kim KR, Munk A, et al (2016) The circular SiZer, inferred persistence of shape parameters and application to early stem cell differentiation. Bernoulli 22(4):2113–2142. 10.3150/15-bej722
  * Johnson et al [1995] Johnson NL, Kotz S, Balakrishnan N (1995) Continuous Univariate Distributions, vol 2, 2nd edn. Wiley, New York
  * Ley and Verdebout [2017] Ley C, Verdebout T (2017) Modern Directional Statistics. Chapman and Hall/CRC, New York, 10.1201/9781315119472
  * Mammen et al [1992] Mammen E, Marron JS, Fisher NI (1992) Some asymptotics for multimodality tests based on kernel density estimates. Probability Theory and Related Fields 91(1):115–132. 10.1007/bf01194493
  * Mardia and Jupp [2000] Mardia KV, Jupp PE (2000) Directional Statistics. John Wiley & Sons, 10.1002/9780470316979
  * Maurer et al [2020] Maurer JD, Huso M, Dalthorp D, et al (2020) Comparing methods to estimate the proportion of turbine-induced bird and bat mortality in the search area under a road and pad search protocol. Environmental and Ecological Statistics 27(4):769–801. 10.1007/s10651-020-00466-0
  * Müller and Sawitzki [1991] Müller DW, Sawitzki G (1991) Excess mass estimates and tests for multimodality. Journal of the American Statistical Association 86(415):738–746. 10.1080/01621459.1991.10475103
  * Oliveira et al [2012] Oliveira M, Crujeiras R, Rodríguez-Casal A (2012) A plug-in rule for bandwidth selection in circular density estimation. Computational Statistics & Data Analysis 56(12):3898–3908. 10.1016/j.csda.2012.05.021
  * Oliveira et al [2014] Oliveira M, Crujeiras RM, Rodríguez-Casal A (2014) CircSiZer: an exploratory tool for circular data. Environmental and Ecological Statistics 21(1):143–159. 10.1007/s10651-013-0249-0
  * OpenStreetMap contributors [2024] OpenStreetMap contributors (2024) Planet dump retrieved from https://planet.osm.org . https://www.openstreetmap.org
  * Pewsey et al [2013] Pewsey A, Neuhäuser M, Ruxton GD (2013) Circular statistics in R. Oxford University Press, Oxford, United Kingdom
  * SenGupta and Arnold [2022] SenGupta A, Arnold BC (eds) (2022) Directional Statistics for Innovative Applications: A Bicentennial Tribute to Florence Nightingale. Springer Nature Singapore, 10.1007/978-981-19-1044-9
  * Silverman [1981] Silverman BW (1981) Using kernel density estimates to investigate multimodality. Journal of the Royal Statistical Society: Series B (Methodological) 43(1):97–99. 10.1111/j.2517-6161.1981.tb01155.x
  * Watson [1961] Watson GS (1961) Goodness-of-fit tests on a circle. Biometrika 48(1/2):109–114. 10.2307/2333135
  * Watson and Williams [1956] Watson GS, Williams EJ (1956) On the construction of significance tests on the circle and the sphere. Biometrika 43(3/4):344–352. 10.2307/2332913
  * Wikimedia Commons [2011] Wikimedia Commons (2011) Topographic and bathymetric blank map of Mexico, from User:Addicted04. URL https://commons.wikimedia.org/wiki/File:Mexico_topographic_map-blank.svg, [Edited]