# The flexible Gumbel distribution: A new model for inference about the mode

Qingyang Liu
Department of Statistics
University of South Carolina
Columbia, SC 29201
<EMAIL_ADDRESS>
& Xianzheng Huang
Department of Statistics
University of South Carolina
Columbia, SC 29201
<EMAIL_ADDRESS>
& Haiming Zhou
Daiichi Sankyo Inc.
Basking Ridge, NJ 07920
<EMAIL_ADDRESS>

###### Abstract

A new unimodal distribution family indexed by the mode and three other
parameters is derived from a mixture of a Gumbel distribution for the maximum
and a Gumbel distribution for the minimum. Properties of the proposed
distribution are explored, including model identifiability and flexibility in
capturing heavy-tailed data that exhibit different directions of skewness over
a wide range. Both frequentist and Bayesian methods are developed to infer
parameters in the new distribution. Simulation studies are conducted to
demonstrate satisfactory performance of both methods. By fitting the proposed
model to simulated data and data from an application in hydrology, it is shown
that the proposed flexible distribution is especially suitable for data
containing extreme values in either direction, with the mode being a location
parameter of interest. A regression model concerning the mode of a response
given covariates based on the proposed unimodal distribution can be easily
formulated, which we apply to data from an application in criminology to
reveal interesting data features that are obscured by outliers. Computer
programs for implementing all considered inference methods in the study are
available at https://github.com/rh8liuqy/flexible_Gumbel.

_Keywords_ extreme values $\cdot$ mixture distribution $\cdot$ modal
regression $\cdot$ unimodal distribution

##  1 Introduction

The mean, median, and mode are three most commonly used measure of central
tendency of data. When data contain outliers that cause heavy tails or are
potentially skewed, the mode is a more sensible representation of the central
location of data than the mean or median. The timely review on mode estimation
and its application by Chacón [2020] and references therein provide many
examples in various fields of research where the mode serves as a more
informative representative value of data. Most existing methods developed to
draw inference for the mode are semi-/non-parametric in nature, starting from
early works on direct estimation in the 1960s [Chernoff, 1964, Dalenius, 1965,
Venter, 1967] to more recent works based on kernel density estimation [Chen,
2018] and quantile-based methods [Ota et al., 2019, Zhang et al., 2021]. There
are two main reasons contributing to the long-lasting trend of opting to
semi-/non-parametric methods for mode estimation, despite the fact that
inference procedures proposed along these veins are usually less
straightforward to implement (e.g., involving bandwidth selection), and less
efficient than their parametric counterparts. First, a parametric model
typically imposes stringent constraints on the relationship between the mode
and other location parameters that may not be satisfied in a given
application. Second, very few existing named distribution families that allow
inclusion of both symmetric and asymmetric distributions in the same family
can be parameterized so that it is indexed by the mode as the location
parameter along with other parameters, such as shape or scale parameters. In
this study, we alleviate concerns raised by both reasons that discourage use
of parametric methods for mode estimation by formulating a flexible
distribution indexed by the (unique) mode and parameters controlling the shape
and scale.

When it comes to modeling heavy-tailed data, the Gumbel distribution [Gumbel,
1941] is arguably one of the most widely used models in many disciplines.
Indeed, as a case of the generalized extreme value distribution [Jenkinson,
1955], the Gumbel distribution for the maximum (or minimum) is well-suited for
modeling extremely large (or small) events that produce heavy-tailed data. For
example, it is often used in hydrology to predict extreme rainfall and flood
frequency [Loaiciga and Leipnik, 1999, Koutsoyiannis, 2004, Dawley et al.,
2019]. In econometrics, the Gumbel distribution plays an important role in
modeling extreme movements of stock prices and large changes in interest rates
[Bali, 2003, Pratiwi et al., 2019]. The Gumbel distribution is indexed by the
mode and a scale parameter, and thus is convenient for mode estimation.
However, the Gumbel distribution for the maximum (or minimum) is right-skewed
(or left-skewed) with the skewness fixed at around 1.44 (or $-1.44$), and the
kurtosis fixed at 5.4 across the entire distribution family. Thus it may be
too rigid for scenarios where the direction and extremeness of outliers
presented in data are initially unclear, or when the direction and level of
skewness are unknown beforehand. Constructions of more flexible distributions
that overcome these limitations have been proposed. In particular, Cooray
[2010] applied a logarithmic transformation on a random variable following the
odd Weibull distribution to obtain the so-called generalized Gumbel
distribution that includes the Gumbel distribution as a subfamily. But the
mode of the generalized Gumbel distribution is not a location parameter this
distribution is indexed by, or an explicit function of other model parameters.
Shin et al. [2015] considered mixture distributions with one of the components
being the Gumbel distribution and the other component(s) being Gumbel of the
same skewness direction or a different distribution, such as the gamma
distribution. Besides the same drawback pointed out for the generalized Gumbel
distribution, it is difficult to formulate a unimodal distribution following
their construction of mixtures, and thus their proposed models are unsuitable
when unimodality is a feature required to make inferring the mode meaningful,
such as in a regression setting, as in modal regression [Yao et al., 2012, Yao
and Li, 2013, Chen, 2018].

With heavy-tailed data in mind and the mode as the location parameter of
interest, we construct a new unimodal distribution that does not impose
stringent constraints on how the mode relates to other central tendency
measures, while allowing a range of kurtosis wide enough to capture heavy
tails at either direction, as well as different degrees and directions of
skewness. This new distribution, called the flexible Gumbel (FG) distribution,
is presented in Section 2, where we study properties of the distribution and
discuss identifiability of the model. We present a frequentist method and a
Bayesian method for estimating parameters in the FG distribution in Section 3.
Finite sample performance of these methods are inspected in simulation study
in Section 4, followed by an application of the FG distribution in hydrology
in Section 5. Section 6 demonstrates fitting a modal regression model based on
the FG distribution to data from a criminology study. Section 7 highlights
contributions of the study and outlines future research directions.

##  2 The flexible Gumbel distribution

The probability density function (pdf) of the Gumbel distribution for the
maximum is given by

$\displaystyle f(x;\theta,\sigma)$
$\displaystyle=\frac{1}{\sigma}\exp\left\\{-\frac{x-\theta}{\sigma}-\exp\left(-\frac{x-\theta}{\sigma}\right)\right\\},$
(1)

where $\theta$ is the mode and $\sigma>0$ is a scale parameter. The pdf of the
Gumbel distribution for the minimum with mode $\theta$ and a scale parameter
$\sigma$ is given by

$\displaystyle f(x;\theta,\sigma)$
$\displaystyle=\frac{1}{\sigma}\exp\left\\{\frac{x-\theta}{\sigma}-\exp\left(\frac{x-\theta}{\sigma}\right)\right\\}.$
(2)

We define a unimodal distribution for a random variable $Y$ via a mixture of
the two Gumbel distributions specified by (1) and (2) that share the same mode
$\theta$ while allowing different scale parameters, $\sigma_{1}$ and
$\sigma_{2}$, in the two components. We call the resultant distribution the
flexible Gumbel distribution, FG for short, with the pdf given by

$\displaystyle f(y)=$ $\displaystyle\
w\times\frac{1}{\sigma_{1}}\exp\left\\{-\frac{x-\theta}{\sigma_{1}}-\exp\left(-\frac{x-\theta}{\sigma_{1}}\right)\right\\}+$
(3) $\displaystyle\
(1-w)\times\frac{1}{\sigma_{2}}\exp\left\\{\frac{x-\theta}{\sigma_{2}}-\exp\left(\frac{x-\theta}{\sigma_{2}}\right)\right\\},$

where $w\in[0,1]$ is the mixing proportion parameter. Henceforth, we state
that $Y\sim\mbox{FG}(\theta,\sigma_{1},\sigma_{2},w)$ if $Y$ follows the
distribution specified by the pdf in (3).

For each component distribution of FG, the mean and median are both some
simple shift of the mode, with each shift solely determined by the scale
parameter. Because the two components in (3) share a common mode $\theta$, the
mode of $Y$ is also $\theta$, and thus the FG distribution is convenient to
use when one aims to infer the mode as a central tendency measure, or to
formulate parametric modal regression models [Bourguignon et al., 2020, Zhou
and Huang, 2020, 2022]. One can easily show that the mean of $Y$ is
$E(Y)=w(\theta+\sigma_{1}\gamma)+(1-w)(\theta-\sigma_{2}\gamma)=\theta+\\{w(\sigma_{1}+\sigma_{2})-\sigma_{2}\\}\gamma$,
where $\gamma\approx 0.5772$ is the Euler-Mascheroni constant. Thus the
discrepancy between the mode and the mean of FG depends on three other
parameters that control the scale and shape of the distribution. The median of
$Y$, denoted by $m$, is the solution to the following equation,

$w\exp\left\\{-\exp\left(-\frac{m-\theta}{\sigma_{1}}\right)\right\\}+(1-w)\left[1-\exp\left\\{-\exp\left(\frac{m-\theta}{\sigma_{2}}\right)\right\\}\right]=0.5.$

Even though this equation cannot be solved for $m$ explicitly to reveal the
median in closed form, it is clear that $m-\theta$ also depends on all three
other parameters of FG. In conclusion, the relationships between the three
central tendency measures of FG are more versatile than those under a Gumbel
distribution for the maximum or a Gumbel distribution for the minimum.

The variance of $Y$ is
$V(Y)=\\{w\sigma_{1}^{2}+(1-w)\sigma_{2}^{2}\\}\pi^{2}/6+w(1-w)(\sigma_{1}+\sigma_{2})^{2}\gamma^{2}$,
which does not depend on the mode parameter $\theta$. Obviously, by setting
$w=0$ or 1, $\mbox{FG}(\theta,\sigma_{1},\sigma_{2},w)$ reduces to one of the
Gumbel components. Unlike a Gumbel distribution that only has one direction of
skewness at a fixed level (of $\pm 1.44)$, an FG distribution can be left-
skewed, or right-skewed, or symmetric. More specifically, with the mode fixed
at zero when studying the skewness and kurtosis of FG, one can show that the
third central moment of $Y$ is given by

$w\bar{w}(\sigma_{1}+\sigma_{2})^{2}\gamma\left\\{\gamma^{2}(\bar{w}-w)(\sigma_{1}+\sigma_{2})+0.5\pi^{2}(\sigma_{1}-\sigma_{2})\right\\}+2\zeta(3)\left(w\sigma_{1}^{3}-\bar{w}\sigma_{2}^{3}\right),$
(4)

where $\bar{w}=1-w$, and $\zeta(3)\approx 1.202$ is Apéry’s constant. Although
the direction of skewness is not immediately clear from (4), one may consider
a special case with $w=0.5$ where (4) reduces to
$(\sigma_{1}-\sigma_{2})\\{\gamma\pi^{2}(\sigma_{1}+\sigma_{2})^{2}/8+\zeta(3)(\sigma_{1}^{2}+\sigma_{1}\sigma_{2}+\sigma_{2}^{2})\\}$.
Now one can see that $\mbox{FG}(\theta,\sigma_{1},\sigma_{2},0.5)$ is
symmetric if and only if $\sigma_{1}=\sigma_{2}$, and it is left-skewed (or
right-skewed) when $\sigma_{1}$ is less (or greater) than $\sigma_{2}$. The
kurtosis of $Y$ can also be derived straightforwardly, with a more lengthy
expression than (4) that we omit here, which may not shed much light on its
magnitude except for that it varies as the scale parameters and the mixing
proportion vary, instead of fixing at 5.4 as for a Gumbel distribution. An R
Shiny app depicting the pdf of $\mbox{FG}(\theta,\sigma_{1},\sigma_{2},w)$
with user-specified parameter values is available at
https://qingyang.shinyapps.io/gumbel_mixture/, created and maintained by the
first author. Along with the density function curve, the Shiny app provides
skewness and kurtosis of the depicted FG density. From there one can see that
the skewness can be much lower than $-1.44$ or higher than 1.44, and the
kurtosis can be much higher than 5.4, suggesting that inference based on FG
can be more robust to outliers than when a Gumbel distribution is assumed for
data at hand, without imposing stringent assumption on the skewness of the
underlying distribution.

The flexibility of a mixture distribution usually comes with concerns relating
to identifiability [Teicher, 1961, 1963, Yakowitz and Spragins, 1968]. In
particular, there is the notorious issue of label switching when fitting a
finite mixture model [Redner and Walker, 1984]. Take the family of two-
component normal mixture (NM) distributions as an example, defined by
$\\{\mbox{NM}(\mu_{1},\sigma_{1},\mu_{2},\sigma_{2},w):\,w\mathcal{N}(\mu_{1},\sigma_{1}^{2})+(1-w)\mathcal{N}(\mu_{1},\sigma_{2}^{2}),\mbox{
for $\sigma_{1},\sigma_{2}>0$ and $w\in[0,1]$}\\}$. When fitting a data set
assuming a normal mixture distribution, one cannot distinguish between, for
instance, $\mbox{NM}(1,2,3,4,0.2)$ and $\mbox{NM}(3,4,1,2,0.8)$, since the
likelihood of the data is identical under these two mixture distributions. As
another example, for data from a normal distribution, a two-component normal
mixture with two identical normal components and an arbitrary mixing
proportion $w\in[0,1]$ leads to the same likelihood, and thus $w$ cannot be
identified. Teicher [1963] showed that imposing an lexicographical order for
the normal components resolves the issue of non-identifiability, which also
excludes mixtures with two identical components in the above normal mixture
family. Unlike normal mixtures of which all components are in the same family
of normal distributions, the FG distribution results from mixing two
components from different families, i.e., a Gumbel distribution for the
maximum and a Gumbel distribution for the minimum, with weight $w$ on the
former component. By construction, FG does not have the label-switching issue.
And, according to Teicher [1963, Theorem 1], the so-constructed mixture
distribution is always identifiable even when the true distribution is a (one-
component) Gumbel distribution.

##  3 Statistical inference

###  3.1 Frequentist inference method

Based on a random sample of size $n$ from the FG distirbution,
$\mathbf{y}=\\{y_{i}\\}_{i=1}^{n}$, maximum likelihood estimators (MLE) of all
model parameters in $\mathbf{\Omega}=(\theta,\sigma_{1},\sigma_{2},w)$ can be
obtained via the expectation-maximization (EM) algorithm [Dempster et al.,
1977]. To apply the EM algorithm, we introduce a latent variable $Z$ that
follows Bernoulli($w$) such that the joint likelihood of $(Y,Z)$ is

$f_{Y,Z}(y,z)=\\{wf_{1}(y;\theta,\sigma_{1})\\}^{z}\\{(1-w)f_{2}(y;\theta,\sigma_{2})\\}^{1-z},$

where $f_{1}(y;\theta,\sigma_{1})$ is the pdf in (1) evaluated at $y$ with the
scale parameter $\sigma=\sigma_{1}$, and $f_{2}(y;\theta,\sigma_{2})$ is the
pdf in (2) evaluated at $y$ with the scale parameter $\sigma=\sigma_{2}$. A
random sample of size $n$ from Bernoulli($w$),
$\mathbf{z}=\\{z_{i}\\}_{i=1}^{n}$, is viewed as missing data, and
$\\{(y_{i},z_{i})\\}_{i=1}^{n}$ are viewed as the complete data in the EM
algorithm. The complete-data log-likelihood is then

$\ell(\mathbf{\Omega};\mathbf{y},\mathbf{z})=\sum_{i=1}^{n}\\{z_{i}\log(wf_{1}(y_{i};\theta,\sigma_{1}))+(1-z_{i})\log((1-w)f_{2}(y_{i};\theta,\sigma_{2}))\\}.$
(5)

Starting from an initial estimate of $\mathbf{\Omega}$ (at the zero-th
iteration), denoted by $\mathbf{\Omega}^{(0)}$, one iterates two steps
referred to as the E-step and the M-step until a convergence criterion is met.
In the E-step at the $(t+1)$-th iteration, one computes the conditional
expectation of (5) given $\mathbf{y}$ while assuming the true parameter value
to be
$\mathbf{\Omega}^{(t)}=(\theta^{(t)},\sigma_{1}^{(t)},\sigma_{2}^{(t)},w^{(t)})$,
that is,
$\operatorname{E}_{\boldsymbol{\Omega}^{(t)}}\\{\ell(\mathbf{\Omega};\mathbf{y},\mathbf{z})|\mathbf{y}\\}$.
This conditional expectation can be shown to be

$Q\left(\boldsymbol{\Omega}\left|\mathbf{\Omega}^{(t)}\right.\right)=\sum_{i=1}^{n}\left\\{T_{i}^{(t)}\operatorname{log}(wf_{1}\left(y_{i};\theta,\sigma_{1}\right))+\left(1-T_{i}^{(t)}\right)\log((1-w)f_{2}\left(y_{i};\theta,\sigma_{2}\right))\right\\},$
(6)

where

$T_{i}^{(t)}=\operatorname{E}_{{\boldsymbol{\Omega}}^{(t)}}(Z|Y=y_{i})=\frac{w^{(t)}f_{1}(y_{i};\theta^{(t)},\sigma_{1}^{(t)})}{w^{(t)}f_{1}(y_{i};\theta^{(t)},\sigma_{1}^{(t)})+(1-w^{(t)})f_{2}(y_{i};\theta^{(t)},\sigma_{2}^{(t)})}.$
(7)

In the M-step at the $(t+1)$-th iteration, one maximizes
$Q(\boldsymbol{\Omega}|\boldsymbol{\Omega}^{(t)})$ with respect to
$\mathbf{\Omega}$ to obtain an updated estimate $\mathbf{\Omega}^{(t+1)}$.

Our experience with the above EM algorithm for fitting the FG distribution
suggests that maximizing $Q(\boldsymbol{\Omega}|\boldsymbol{\Omega}^{(t)})$ in
(6) can be numerically challenging. We thus exploit the expectation-
conditional maximization (ECM) algorithm [Meng and Rubin, 1993], which
replaces the M-step with a sequence of simpler conditional maximizations
referred to as the CM-step. More specifically, in the CM-step, the updating
formula for $w$ is simply $w^{(t+1)}=\sum_{i=1}^{n}T_{i}^{(t)}/n.$ There is no
closed-form updating formula for the other three parameters in
$\mathbf{\Omega}$, but they can now be easily updated by most well-accepted
one-dimensional optimization algorithms, such as the Newton-Raphson method. To
ensure convergence at the global maximum, as recommended by Wu [1983], one
should implement the ECM algorithm several rounds with different starting
values $\mathbf{\Omega}^{(0)}$.

After obtaining the MLE of $\mathbf{\Omega}$, denoted by
$\hat{\mathbf{\Omega}}$, we propose to use the sandwich variance estimator
[Boos and Stefanski, 2013, Chapter 7] to estimate the variance-covariance
matrix of $\hat{\mathbf{\Omega}}$. One may also estimate the variance-
covariance of $\hat{\mathbf{\Omega}}$ based on the observed information matrix
as described in Louis [1982] and Oakes [1999]. The benefit of using the
sandwich variance estimator is its robustness to model misspecification.
Finally, the EM and ECM algorithms bear a strong resemblance to data
augmentation [Wei and Tanner, 1990] in the Bayesian framework, which we turn
to next for inferring $\mathbf{\Omega}$.

###  3.2 Bayesian inference method

In the Bayesian framework, we formulate hierarchical models starting with the
FG distribution,

$Y|\theta,\sigma_{1},\sigma_{2},w\sim\text{FG}(\theta,\sigma_{1},\sigma_{2},w),$

followed by independent weakly informative or non-informative priors for
elements in $\mathbf{\Omega}$,

$\displaystyle\theta$ $\displaystyle\sim\mathcal{N}(0,10^{4}),$
$\displaystyle\sigma_{j}$ $\displaystyle\sim\text{inv-Gamma}(1,1),\text{ for
}j=1,2,$ $\displaystyle w$ $\displaystyle\sim\text{Beta}(1,1),$

where inv-Gamma refers to the inverse Gamma distribution. We choose the above
prior for the scale parameters by following the prior selection for variance
parameters suggested in Gelman [2006]

We employ the Metropolis-within-Gibbs sampler [Müller, 1991, 1993] to obtain
an estimate of $\mathbf{\Omega}$ from the posterior distribution of
$\mathbf{\Omega}$ given observed data $\mathbf{y}$. Similar to the EM/ECM
algorithm in Section 3.1, the latent variable $Z$ is also introduced as a
device to carry out data augmentation. And the iterative algorithm presented
next is based on the following two conditional distributions that can be
easily proved,

$\displaystyle
z_{i}|\theta,\sigma_{1},\sigma_{2},w,\mathbf{z}_{-i},\mathbf{y}$
$\displaystyle\sim\text{Bernoulli}\left(\frac{wf_{1}(y_{i};\theta,\sigma_{1})}{wf_{1}(y_{i};\theta,\sigma_{1})+(1-w)f_{2}(y_{i};\theta,\sigma_{2})}\right),$
$\displaystyle w|\theta,\sigma_{1},\sigma_{2},\mathbf{z},\mathbf{y}$
$\displaystyle\sim\text{Beta}\left(1+\sum_{i=1}^{n}z_{i},\,n+1-\sum_{i=1}^{n}z_{i}\right),$

where $\mathbf{z}_{-i}$ results from dropping $z_{i}$ from $\mathbf{z}$, and
the first result above is also from which (7) is deduced.

The Metropolis-within-Gibbs sampler at the $(t+1)$-th iteration involves four
steps outlined below.

  * •

Step 1: For $i=1,\dots,n$, draw $z^{(t+1)}_{i}$ from
$\mbox{Bernoulli}(T_{i}^{(t)})$, where $T_{i}^{(t)}$ is given in (7).

  * •

Step 2: Draw $w^{(t+1)}$ from
$\mbox{Beta}\left(1+\sum_{i=1}^{n}z^{(t+1)}_{i},\,n+1-\sum_{i=1}^{n}z^{(t+1)}_{i}\right).$

  * •

Step 3: Draw $\tilde{\theta}$ from $\mathcal{N}(\theta^{(t)},\tau_{0})$, and
update $\theta^{(t)}$ to $\theta^{(t+1)}$ according to the following decision
rule,

$\theta^{(t+1)}=\begin{cases}\tilde{\theta},\quad\text{with probability
}\displaystyle{q=\min\left\\{\frac{p(\tilde{\theta}|w^{(t+1)},\sigma^{(t)}_{1},\sigma^{(t)}_{2},\mathbf{y})}{p(\theta^{(t)}|w^{(t+1)},\sigma^{(t)}_{1},\sigma^{(t)}_{2},\mathbf{y})},\,1\right\\}},\\\
\theta^{(t)},\quad\text{with probability $1-q$}.\end{cases}$

  * •

Step 4: For $j=1,2$, draw $\tilde{\sigma}_{j}$ from
$\mathcal{N}(\sigma_{j}^{(t)},\tau_{j})$, and update $\sigma_{j}^{(t)}$ to
$\sigma_{j}^{(t+1)}$ according to the following decision rule, for $k\neq j$,

$\sigma_{j}^{(t+1)}=\begin{cases}\tilde{\sigma}_{j},\quad\text{with
probability
}\displaystyle{q=\min\left\\{\frac{p(\tilde{\sigma}_{j}|\theta^{(t+1)},\sigma^{(t)}_{k},w^{(t+1)},\mathbf{y})}{p(\sigma_{j}^{(t)}|\theta^{(t+1)},\sigma^{(t)}_{k},w^{(t+1)},\mathbf{y})},\,1\right\\}},\\\
\sigma_{j}^{(t)},\quad\text{with probability $1-q$}.\end{cases}$

In Steps 3 and 4, $p(\cdot|\cdot)$ refers to a conditional pdf generically,
$\tau_{0}$, $\tau_{1}$, and $\tau_{2}$ are three positive tuning parameters
whose values should be chosen so that the acceptance rate at each step is
around $23\%$ [Gelman et al., 1997]. To draw samples from the joint posterior
distribution, there are numerous ways to design the Markov chain Monte Carlo
(MCMC) sampler. Instead of the Metropolis-within-Gibbs sampler we adopt here,
one may use other existing MCMC software, such as Stan [Stan Development Team,
2021], JAGS [Plummer et al., 2003], and BUGS [Spiegelhalter et al., 1996, Lunn
et al., 2009], two of which we demonstrate in the Appendix. After obtaining
enough high quality samples from the joint posterior distribution
$p(\theta,\sigma_{1},\sigma_{2},w|\mathbf{y})$, Bayesian inference is
straightforward, including point estimation, interval estimation, and
uncertainty assessment.

##  4 Simulation study

Large-sample properties of MLEs and likelihood-based Bayesian inference under
a correct model for data have been well studied. To assess finite-sample
performance of the frequentist method and Bayesian method proposed in Section
3, we carried out simulation study with two specific aims: first, to compare
inference results from the two methods; second, to compare goodness of fit for
data from distributions outside of the FG family when one assumes an FG
distribution and when one assumes a two-component normal mixture distribution
for the data.

In the first experiment, referred to as (E1) in the sequel, we drew a random
sample of size $n\in\\{100,200\\}$ from an FG distribution with $\theta=0$,
$\sigma_{1}=1$, $\sigma_{2}=5$, and $w=0.5$. Based on each simulated data set,
we estimated $\mathbf{\Omega}$ by applying the ECM algorithm and the
Metropolis-within-Gibbs algorithm. The former algorithm produced the MLE of
$\mathbf{\Omega}$, and we used the median of the posterior distribution of
$\mathbf{\Omega}$ at convergence of the latter algorithm as another point
estimate of $\mathbf{\Omega}$. Table 1 presents summary statistics of these
estimates of $\mathbf{\Omega}$ and estimates of the corresponding standard
deviation across 1000 Monte Carlo replicates.

Table 1: Frequentist and Bayesian inference results in experiment (E1) across 1000 Monte Carlo replicates. Here, point.est stands for the average of 1000 point estimates for each parameter from each method, $\widehat{\text{s.d.}}$ stands for the average of the corresponding 1000 estimated standard deviations, and s.d. refers to the empirical standard deviation of the 1000 point estimates from each method. Numbers in parentheses are $100\times$ Monte Carlo standard errors associated with the averages. |  | Frequentist | Bayesian
---|---|---|---
sample size | parameter | point.est | $\widehat{\text{s.d.}}$ | s.d. | point.est | $\widehat{\text{s.d.}}$ | s.d.
| $\theta$ | 0.002 | 0.198 (0.20) | 0.201 | 0.013 | 0.205 (0.15) | 0.203
| $\sigma_{1}$ | 0.979 | 0.204 (0.41) | 0.216 | 1.014 | 0.224 (0.27) | 0.214
| $\sigma_{2}$ | 4.932 | 0.590 (0.56) | 0.613 | 4.813 | 0.666 (0.44) | 0.615
$n=100$ | $w$ | 0.495 | 0.091 (0.09) | 0.090 | 0.484 | 0.090 (0.04) | 0.088
| $\theta$ | 0.008 | 0.136 (0.08) | 0.129 | 0.011 | 0.137 (0.07) | 0.130
| $\sigma_{1}$ | 0.999 | 0.143 (0.21) | 0.144 | 1.013 | 0.144 (0.10) | 0.141
| $\sigma_{2}$ | 4.993 | 0.435 (0.32) | 0.431 | 4.940 | 0.457 (0.20) | 0.434
$n=200$ | $w$ | 0.500 | 0.064 (0.04) | 0.063 | 0.495 | 0.063 (0.02) | 0.062

According to Table 1, all estimates for parameters in $\mathbf{\Omega}$ are
reasonably close to the truth. A closer inspection on the reported empirical
mean of these estimates along with their empirical standard error suggests
that, when $n=100$, the Bayesian method may slightly underestimate
$\sigma_{2}$, the larger of the two scale parameters of FG. We believe that
this is due to the inverse gamma prior imposed on the scale parameters that is
sharply peaked near zero, and thus the posterior median of the larger scale
parameter tends to be pulled downwards when the sample size is not large. As
the sample size increases to $200$, this trend of underestimation appears to
diminish. The empirical means of the standard deviation estimates from both
methods are close to the corresponding empirical standard deviations, which
indicate that the variability of a point estimator is accurately estimated,
whether it is based on the sandwich variance estimator in the frequentist
framework, or based on the posterior sampling in the Bayesian framework. In
summary, the methods proposed in Section 3 under both frameworks provide
reliable inference for $\mathbf{\Omega}$ along with accurate uncertainty
assessment of the point estimators when data arise from an FG distribution.

Among all existing mixture distributions, normal mixtures probably have the
longest history and are most referenced in the literature. In another
experiment, we compared the model fitting of normal mixture with that of FG
when data arise from three heavy-tailed distributions: (E2) Laplace with the
location parameter equal to zero and the scale parameter equal to 2; (E3) a
mixture of two Gumbel distributions for the maximum, with a common mode at
zero, scale parameters in the two components equal to 2 and 6, respectively,
and the mixing proportion equal to 0.5; (E4) a $t$ distribution with degrees
of freedom equal to 5. From each of the three distributions in (E2)–(E4), we
generated a random sample of size $n=200$, following which we fit a two-
component normal mixture model via the EM algorithm implemented using the R
package mixtools, and also fit an FG model via the two algorithms described in
Section 3. This model fitting exercise was repeated for 1000 Monte Carlo
replicates under each of (E2)–(E4).

We used an empirical version of the Kullback-Leibler divergence as the metric
to assess the quality of modeling fitting. We denote the true density function
as $p(\cdot)$, and let $\hat{p}(\cdot)$ be a generic estimated density
resulting from one of the three considered model fitting strategies. Under
each setting in (E2)–(E4), a random sample of size 50000,
$(x_{1},\ldots,x_{50000})$, were generated from the true distribution, and an
empirical version of the Kullback-Leibler divergence from $\hat{p}(\cdot)$ to
$p(\cdot)$ is given by
$D_{\hbox{\mbox{KL}}}=(1/50000)\sum_{i=1}^{50000}\log(p(x_{i})/\hat{p}(x_{i}))$.
Figure 1 shows the boxplots of $D_{\hbox{\mbox{KL}}}$ across 1000 Monte Carlo
replicates corresponding to each model fitting scheme under (E2)–(E4).

Figure 1: Boxplots of the empirical Kullback-Leibler divergence from an
estimated density to the true density under each of the true-model settings in
(E2)–(E4). Under each setting, the three considered model fitting strategies
are, from left to right in the figure, (i) using the ECM algorithm to fit an
FG distribution (FG ECM), (ii) using the Bayesian method to fit an FG
distribution (FG Bayes), and (iii) using the EM algorithm to fit a normal
mixture distribution (Normal Mixture Distribution EM).

Judging from Figure 1, the FG distribution clearly outperform the normal
mixture when fitting data from any of the three heavy-tailed distributions in
(E2)–(E4), and results from the frequentist method are comparable with those
from the Bayesian method for fitting an FG model. When implementing the ECM
algorithm for fitting the FG model and the EM algorithm for fitting the normal
mixture, we set a maximum number of iterations at 1000. Our ECM algorithm
always converged in the simulation, i.e., converged to a stationary point
within 1000 iterations. But the EM algorithm for fitting a normal mixture
often had trouble achieving that, with more difficulty when data come from a
heavier-tailed distribution. More specifically, under (E4), which has the
highest kurtosis (equal to 9) among the three settings, the EM algorithm
failed to converge in 59.9% of all Monte Carlo replicates; under (E2), which
has the second highest kurtosis (equal to 6), it failed to converge in 6.7% of
the replicates. Results associated with the normal mixture from these failing
replicates were not included when producing the boxplots in Figure 1. In
conclusion, the FG distribution is more suitable for symmetric or asymmetric
heavy-tailed data than the normal mixture distribution.

##  5 An application in hydrology

Daily maximum water elevation changes of a waterbody, such as ocean, lake, and
wetland, are of interest in hydrologic research. These changes may be close to
zero in most days, but can be extremely large or small under extreme weather.
From National Water Information System (https://waterdata.usgs.gov/), we
downloaded water elevation data for Lake Murray near Columbia, South Carolina,
United States, recorded from September 18, 2020 to September 18, 2021. The
water elevation change of a given day was calculated by contrasting the
maximum elevation and the minimum elevation on that day, returning a positive
(negative) value if the maximum record of the day comes after (before) the
minimum record on the same day. We fit the FG distribution to the resultant
data with $n=366$ records using the frequentist method and the Bayesian
method, with results presented in Table 2. The two inference methods produced
very similar estimates for most parameters, although small differences were
observed. For example, one would estimate the mode of daily maximum water
elevation change to be $-0.795$ feet based on the frequentist method, but
estimate it to be $-0.486$ feet using the Bayesian method. The discrepancy
between these two mode estimates is minimal considering that the daily maximum
water elevation changes range from $-38$ feet to 49.4 feet within this one-
year period. In fact, taking into account the uncertainty in these point
estimates, we do not interpret any of these differences as statistically
significant because a parameter estimate from one method always falls in the
interval estimate for the same parameter from the other method according to
Table 2. Using parameter estimates in Table 2 in the aforementioned R Shiny
app, we obtained an estimated skewness of $-0.102$ and an estimated kurtosis
of 6.384 based on the frequentist inference results, whereas the Bayes
inference yielded an estimated skeweness of 0.058 and an estimated kurtosis of
6.074. Combining these two sets of results, we concluded that the underlying
distribution of daily maximum water elevation change may be nearly symmetry,
with outliers on both tails that cause tails heavier than that of a Gumbel
distribution.

Table 2: Frequentist and Bayesian inferences about daily maximum water elevation changes of Lake Murray, South Carolina, United States. Besides parameter estimates (under point.est) and the estimated standard deviations of these parameter estimates (under $\widehat{\mbox{s.d.}}$), 95% confidence intervals of the parameters from the frequentist method, and 95% credible intervals from the Bayesian method are also provided (under lower 95 and upper 95). | Frequentist | Bayesian
---|---|---
parameter | point.est | $\widehat{\text{s.d.}}$ | lower 95 | upper 95 | point.est | $\widehat{\text{s.d.}}$ | lower 95 | upper 95
$\theta$ | -0.795 | 0.796 | -2.355 | 0.764 | -0.486 | 0.694 | -1.679 | 0.973
$\sigma_{1}$ | 5.186 | 0.541 | 4.124 | 6.247 | 5.399 | 0.651 | 4.534 | 6.916
$\sigma_{2}$ | 6.237 | 1.735 | 2.836 | 9.638 | 5.734 | 1.031 | 4.390 | 8.029
$w$ | 0.698 | 0.169 | 0.367 | 1.029 | 0.630 | 0.141 | 0.329 | 0.847

Figure 2 presents the estimated density functions from these two methods, in
contrast with the estimated density curve resulting from fitting the data to a
two-component normal mixture, and a kernel density estimate using a Gaussian
kernel with the bandwidth selected according to the method proposed by
Sheather and Jones [1991]. The last estimate is fully nonparametric and served
as a benchmark against which the other three density estimates were assessed
graphically. The kernel density estimate is more flexible at describing
varying tail behaviors, but such flexibility comes at the cost of statistical
efficiency and interpretability. With the wiggly tails evident in Figure 2 for
this estimate, we suspected certain level of overfitting of the kernel density
estimate. This often happens to kernel-based estimation of a function around a
region where data are scarce, with a bandwidth not large enough for the
region. Between the two FG density estimates, the difference is almost
negligible. They both track the kernel density estimate closely over a wide
range of the support around the mode. The mode of the estimated normal mixture
density is close to the other three mode estimates, but the tails are much
lighter than those of the other three estimated densities.

Figure 2: Four density estimates based on daily maximum water elevation
changes in Lake Murray, including the kernel density estimate (solid line),
the estimated FG density from the ECM algorithm (dotted line), the estimated
FG density from the Bayesian method (dashed line), and the estimated normal
mixture density (dash-dotted line).

Besides comparing the three parametric density estimates pictorially, we also
used the Monte-Carlo based one-sample Kolmogorov–Smirnov test to assess the
goodness of fit. The $p$-values from this test are 0.223, 0.312, and 0.106 for
the frequentist FG density estimate, the Bayesian FG density estimate, and the
estimated normal mixture density, respectively. Although none of the
$p$-values are low enough to indicate lack of fit (at significance level 0.05
for example), the $p$-value associated with the normal mixture is much lower
than those for FG. This provides quantitative evidence that an FG distribution
fits the current data better than a normal mixture. It is also worth noting
that the Kolmogorov–Smirnov test is known to have low power to detect
deviations from a posited distribution that occur in the tails [Mason and
Schuenemeyer, 1983]. This may explain the above-0.05 $p$-value for the normal
mixture fit of the data even though the tail of this posited distribution may
be too thin for the current data.

We used Stan to implement the Bayesian inference for the Lake Murray data, and
the code and posterior output are given in the Appendix. The output provided
there indicates that our MCMC chain has converged (see the Rhat statistics).
The JAGS code for fitting the FG distribution is also given in the Appendix.

##  6 An application in criminology

Table 3: Frequentist and Bayesian modal regression models based on the FG distribution fitted to the crime data. Besides parameter estimates (under point.est) and the estimated standard deviations of these parameter estimates (under $\widehat{\mbox{s.d.}}$), 95% confidence intervals of the parameters from the frequentist method, and 95% credible intervals from the Bayesian method are also provided (under lower 95 and upper 95). | Frequentist | Bayesian
---|---|---
parameter | point.est | $\widehat{\text{s.d.}}$ | lower 95 | upper 95 | point.est | $\widehat{\text{s.d.}}$ | lower 95 | upper 95
$\beta_{1}$ | -0.166 | 0.071 | -0.306 | -0.026 | -0.162 | 0.078 | -0.315 | -0.008
$\beta_{2}$ | 0.217 | 0.110 | 0.001 | 0.433 | 0.231 | 0.123 | -0.009 | 0.475
$\beta_{3}$ | 0.067 | 0.013 | 0.042 | 0.093 | 0.067 | 0.014 | 0.039 | 0.095
Table 4: Mean regression model based on the normal distribution fitted to the crime data. Besides parameter estimates (under point.est) and the estimated standard deviations of these parameter estimates (under $\widehat{\mbox{s.d.}}$), 95% confidence intervals of the parameters parameter | point.est | $\widehat{\text{s.d.}}$ | lower 95 | upper 95
---|---|---|---|---
$\beta_{1}$ | 0.467 | 0.161 | 0.142 | 0.792
$\beta_{2}$ | 1.140 | 0.224 | 0.689 | 1.591
$\beta_{3}$ | 0.068 | 0.034 | 0.000 | 0.136

With the location parameter $\theta$ signified in the FG distribution as the
mode, it is straightforward to formulate a modal regression model that
explores the relationship between the response variable and predictors. To
demonstrate the formulation of a modal regression model based on the FG
distribution, we analyze a data set from Agresti et al. [2021] in the area of
criminology. This data set contains the percentage of college education,
poverty percentage, metropolitan rate, and murder rate for the 50 states in
the United States and the District of Columbia from year 2003. The poverty
percentage is the percentage of the residents with income below the poverty
level; the metropolitan rate is defined as the percentage of population living
in the metropolitan area; and the murder rate is the annual number of murders
per $100,000$ people in the population.

We fit the following modal regression model to investigate the association
between the murder rate ($Y$) and the aforementioned demographic variables,

$Y\mid\boldsymbol{\beta},\sigma_{1},\sigma_{2}\sim\operatorname{FG}(\beta_{0}+\beta_{1}\times\text{
college }+\beta_{2}\times\text{ poverty }+\beta_{3}\times\text{ metropolitan
},\sigma_{1},\sigma_{2},w),$

where $\boldsymbol{\beta}=[\beta_{0},\beta_{1},\beta_{2},\beta_{3}]^{\top}$
includes all regression coefficients. For the prior elicitation in Bayesian
inference, we assume that
$\beta_{0},\dots,\beta_{3}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}\mathcal{N}(0,10^{4})$
and use the same priors for $\sigma_{1}$, $\sigma_{2}$ and $w$ as those in
Section 3.2. As a more conventional regression analysis to compare with our
modal regression, we also fit the mean regression model assuming mean-zero
normal model error to the data.

Table 3 shows the inference results from the modal regression model, and Table
4 presents the inference results from the mean regression model. At $5\%$
significance level, both frequentist and Bayesian modal regression analyses
confirm that there exists a negative association between the percentage of
college education and the murder rate, as well as a positive association
between the metropolitan rate and the murder rate. In contrast, according to
the inferred mean regression model, there is a positive association between
the percentage of college education and the murder rate. Such claimed positive
association is intuitively difficult to justify and contradicts with many
published results in criminology [Hjalmarsson and Lochner, 2012, Lochner,
2020].

The scatter plot of the data in Figure 3 can shed some light on why one
reaches to such a drastically different conclusion on a covariate effect when
mean regression is considered in place of modal regression. As shown in Figure
3, the exists an obvious outlier, District of Columbia (DC), in panels of the
first row of the scatter plot matrix for instance. Mean regression reacts to
this one extreme outlier by inflating the covariate effect associated with the
percentage of college education in the inferred mean regression function.
Thanks to the heavy-tailed feature of the FG distribution, modal regression
based on this distribution is robust to outliers, which strives to capture
data features suggested by majority of the data and is not distracted by the
extreme outlier when inferring covariate effects in this application.

Figure 3: Scatter plot matrix of the crime data.

##  7 Discussion

The mode had been an overlooked location parameter in statistical inference
until recently when the statistics community witnessed a revived interest in
modal regression among statisticians [Chen, 2018, Chacón, 2020, Feng et al.,
2020, Xu et al., 2020, Ullah et al., 2021, Wang and Li, 2021, Xiang and Yao,
2022]. Historically, statistical inference for the mode have been mostly
developed under the nonparametric framework for reasons we point out in
Section 1. Existing semiparametric methods for modal regression only introduce
parametric ingredients in the regression function, i.e., the conditional mode
of the response, with the mode-zero error distribution left in a nonparametric
form [Yao and Li, 2013, Liu et al., 2013, Zhang et al., 2013, Yang and Yang,
2014, Zhao et al., 2014, Krief, 2017, Tian et al., 2017, Li and Huang, 2019].
The few recently proposed parametric modal regression models all impose
stringent parametric assumptions on the error distribution [Bourguignon et
al., 2020, Zhou and Huang, 2020, 2022]. Our proposed flexible Gumbel
distribution greatly alleviates concerns contributing to data scientists’
reluctance to adopt a parametric framework when drawing inference for the
mode. This new distribution is a heterogeneous mixture in the sense that the
two components in the mixture belong to different Gumbel distribution
families, which is a feature that shields it from the non-identifiability
issue most traditional mixture distributions face, such as the normal
mixtures. The proposed distribution is indexed by the mode along with shape
and scale parameters, and thus is convenient to use to draw inference for the
mode while remaining flexible. It is also especially suitable for modeling
heavy-tailed data, whether the heaviness in tails is due to extremely large or
extremely small observations, or both. These are virtues of FG that cannot be
achieved by the popular normal mixture and many other existing mixture
distributions.

We develop the numerically efficient and stable ECM algorithm for frequentist
inference for the FG distribution, and a reliable Bayesian inference method
that can be easily implemented using free software, including Stan, JAGS, and
BUGS. Compared with the more widely adopted mean regression framework, the
modal regression model based on FG we entertained in Section 6 shows great
potential in revealing meaningful covariate effects potentially masked by
extreme outliers. With these advances made in this study, we open up new
directions for parametric modal regression and semiparametric modal regression
with a fully parametric yet flexible error distribution, and potentially
nonparametric ingredients incorporated in the regression function.

## Disclosure statement

Computer programs for implementing the FG distribution, related models and
data used in this paper are available at
https://github.com/rh8liuqy/flexible_Gumbel.

## References

  * Chacón [2020] José E. Chacón.  The modal age of statistics.  _International Statistical Review_ , 88(1):122–141, 2020.  doi:10.1111/insr.12340.
  * Chernoff [1964] Herman Chernoff.  Estimation of the mode.  _Annals of the Institute of Statistical Mathematics_ , 16(1):31–41, 1964.
  * Dalenius [1965] Tore Dalenius.  The mode–a neglected statistical parameter.  _Journal of the Royal Statistical Society. Series A (General)_ , 128(1):110–117, 1965.
  * Venter [1967] JH Venter.  On estimation of the mode.  _The Annals of Mathematical Statistics_ , 38(5):1446–1455, 1967.
  * Chen [2018] Yen-Chi Chen.  Modal regression using kernel density estimation: A review.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 10(4):e1431, 2018.
  * Ota et al. [2019] Hirofumi Ota, Kengo Kato, and Satoshi Hara.  Quantile regression approach to conditional mode estimation.  _Electronic Journal of Statistics_ , 13(2):3120–3160, 2019.
  * Zhang et al. [2021] Tao Zhang, Kengo Kato, and David Ruppert.  Bootstrap inference for quantile-based modal regression.  _Journal of the American Statistical Association_ , pages 1–13, 2021\.
  * Gumbel [1941] E. J. Gumbel.  The return period of flood flows.  _The Annals of Mathematical Statistics_ , 12(2):163–190, June 1941.  doi:10.1214/aoms/1177731747.
  * Jenkinson [1955] Arthur F Jenkinson.  The frequency distribution of the annual maximum (or minimum) values of meteorological elements.  _Quarterly Journal of the Royal Meteorological Society_ , 81(348):158–171, 1955.
  * Loaiciga and Leipnik [1999] H. A. Loaiciga and R. B. Leipnik.  Analysis of extreme hydrologic events with gumbel distributions: marginal and additive cases.  _Stochastic Environmental Research and Risk Assessment (SERRA)_ , 13(4):251–259, August 1999.  doi:10.1007/s004770050042.
  * Koutsoyiannis [2004] Demetris Koutsoyiannis.  Statistics of extremes and estimation of extreme rainfall: I. theoretical investigation / statistiques de valeurs extrêmes et estimation de précipitations extrêmes: I. recherche théorique.  _Hydrological Sciences Journal_ , 49(4), August 2004.  doi:10.1623/hysj.49.4.575.54430.
  * Dawley et al. [2019] Shawn Dawley, Yong Zhang, Xiaoting Liu, Peng Jiang, Geoffrey Tick, HongGuang Sun, Chunmiao Zheng, and Li Chen.  Statistical analysis of extreme events in precipitation, stream discharge, and groundwater head fluctuation: distribution, memory, and correlation.  _Water_ , 11(4):707, April 2019.  doi:10.3390/w11040707.
  * Bali [2003] Turan G. Bali.  An extreme value approach to estimating volatility and value at risk.  _The Journal of Business_ , 76(1):83–108, January 2003.  doi:10.1086/344669.
  * Pratiwi et al. [2019] N Pratiwi, C Iswahyudi, and R I Safitri.  Generalized extreme value distribution for value at risk analysis on gold price.  _Journal of Physics: Conference Series_ , 1217(1):012090, May 2019.  doi:10.1088/1742-6596/1217/1/012090.
  * Cooray [2010] Kahadawala Cooray.  Generalized Gumbel distribution.  _Journal of Applied Statistics_ , 37(1):171–179, 2010.
  * Shin et al. [2015] Ju-Young Shin, Taesam Lee, and Taha B. M. J. Ouarda.  Heterogeneous mixture distributions for modeling multisource extreme rainfalls.  _Journal of Hydrometeorology_ , 16(6):2639–2657, November 2015.  doi:10.1175/jhm-d-14-0130.1.
  * Yao et al. [2012] Weixin Yao, Bruce G Lindsay, and Runze Li.  Local modal regression.  _Journal of nonparametric statistics_ , 24(3):647–663, 2012.
  * Yao and Li [2013] Weixin Yao and Longhai Li.  A new regression model: Modal linear regression.  _Scandinavian Journal of Statistics_ , 41(3):656–671, 2013.  doi:10.1111/sjos.12054.
  * Bourguignon et al. [2020] Marcelo Bourguignon, Jeremias Leão, and Diego I Gallardo.  Parametric modal regression with varying precision.  _Biometrical Journal_ , 62(1):202–220, 2020.
  * Zhou and Huang [2020] Haiming Zhou and Xianzheng Huang.  Parametric mode regression for bounded responses.  _Biometrical Journal_ , 62(7):1791–1809, 2020\.  doi:10.1002/bimj.202000039.
  * Zhou and Huang [2022] Haiming Zhou and Xianzheng Huang.  Bayesian beta regression for bounded responses with unknown supports.  _Computational Statistics & Data Analysis_, 167:107345, 2022.
  * Teicher [1961] Henry Teicher.  Identifiability of mixtures.  _The Annals of Mathematical Statistics_ , 32(1):244–248, 1961.
  * Teicher [1963] Henry Teicher.  Identifiability of finite mixtures.  _The Annals of Mathematical Statistics_ , 34(4):1265–1269, December 1963.  doi:10.1214/aoms/1177703862.
  * Yakowitz and Spragins [1968] Sidney J Yakowitz and John D Spragins.  On the identifiability of finite mixtures.  _The Annals of Mathematical Statistics_ , 39(1):209–214, 1968.
  * Redner and Walker [1984] Richard A Redner and Homer F Walker.  Mixture densities, maximum likelihood and the em algorithm.  _SIAM review_ , 26(2):195–239, 1984.
  * Dempster et al. [1977] A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, September 1977.  doi:10.1111/j.2517-6161.1977.tb01600.x.
  * Meng and Rubin [1993] Xiao-Li Meng and Donald B. Rubin.  Maximum likelihood estimation via the ECM algorithm: A general framework.  _Biometrika_ , 80(2):267–278, 1993.  doi:10.1093/biomet/80.2.267.
  * Wu [1983] C. F. Jeff Wu.  On the convergence properties of the EM algorithm.  _The Annals of Statistics_ , 11(1), March 1983.  doi:10.1214/aos/1176346060.
  * Boos and Stefanski [2013] Dennis D Boos and L. A Stefanski.  _Essential Statistical Inference_.  Springer New York, 2013.  doi:10.1007/978-1-4614-4818-1.
  * Louis [1982] Thomas A. Louis.  Finding the observed information matrix when using the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 44(2):226–233, January 1982.  doi:10.1111/j.2517-6161.1982.tb01203.x.
  * Oakes [1999] D. Oakes.  Direct calculation of the information matrix via the EM.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 61(2):479–482, April 1999.  doi:10.1111/1467-9868.00188.
  * Wei and Tanner [1990] Greg CG Wei and Martin A Tanner.  A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms.  _Journal of the American Statistical Association_ , 85(411):699–704, 1990.
  * Gelman [2006] Andrew Gelman.  Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).  _Bayesian Analysis_ , 1(3):515–534, September 2006.  doi:10.1214/06-ba117a.
  * Müller [1991] Peter Müller.  A generic approach to posterior integration and Gibbs sampling.  _Technical report, Purdue University, West Lafayette, Indiana_ , 1991\.
  * Müller [1993] Peter Müller.  Alternatives to the Gibbs sampling scheme.  _Technical report, Institue of Statistics and Decison Sciences, Duke University_ , 1993.
  * Gelman et al. [1997] A. Gelman, W. R. Gilks, and G. O. Roberts.  Weak convergence and optimal scaling of random walk metropolis algorithms.  _The Annals of Applied Probability_ , 7(1), February 1997\.  doi:10.1214/aoap/1034625254.
  * Stan Development Team [2021] Stan Development Team.  RStan: the R interface to Stan, 2021.  URL https://mc-stan.org/.  R package version 2.21.3.
  * Plummer et al. [2003] Martyn Plummer et al.  JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling.  In _Proceedings of the 3rd international workshop on distributed statistical computing_ , volume 124, pages 1–10. Vienna, Austria., 2003.
  * Spiegelhalter et al. [1996] David Spiegelhalter, Andrew Thomas, Nicky Best, and Wally Gilks.  BUGS 0.5: Bayesian inference using Gibbs sampling manual (version ii).  _MRC Biostatistics Unit, Institute of Public Health, Cambridge, UK_ , pages 1–59, 1996.
  * Lunn et al. [2009] David Lunn, David Spiegelhalter, Andrew Thomas, and Nicky Best.  The BUGS project: Evolution, critique and future directions.  _Statistics in medicine_ , 28(25):3049–3067, 2009\.
  * Sheather and Jones [1991] S.J. Sheather and Chris Jones.  A reliable data-based bandwidth selection method for kernel density estimation.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 53(3):683–690, 1991.
  * Mason and Schuenemeyer [1983] David M Mason and John H Schuenemeyer.  A modified Kolmogorov-Smirnov test sensitive to tail alternatives.  _The annals of Statistics_ , pages 933–946, 1983.
  * Agresti et al. [2021] Alan Agresti, Christine Franklin, and Bernhard Klingenberg.  _Statistics: The Art and Science of Learning from Data_.  Pearson Education, 5 edition, 2021.
  * Hjalmarsson and Lochner [2012] Randi Hjalmarsson and Lance Lochner.  The impact of education on crime: International evidence.  _CESifo DICE Report_ , 10(2):49–55, 2012.
  * Lochner [2020] Lance Lochner.  Education and crime.  In _The Economics of Education_ , pages 109–117. Elsevier, 2020.
  * Feng et al. [2020] Yunlong Feng, Jun Fan, and Johan Suykens.  A statistical learning approach to modal regression.  _Journal of Machine Learning Research_ , 21(2):1–35, 2020.
  * Xu et al. [2020] Jiamiao Xu, Fangzhao Wang, Qinmu Peng, Xinge You, Shuo Wang, Xiao-Yuan Jing, and CL Philip Chen.  Modal-regression-based structured low-rank matrix recovery for multiview learning.  _IEEE Transactions on Neural Networks and Learning Systems_ , 32(3):1204–1216, 2020.
  * Ullah et al. [2021] Aman Ullah, Tao Wang, and Weixin Yao.  Modal regression for fixed effects panel data.  _Empirical Economics_ , 60(1):261–308, 2021.
  * Wang and Li [2021] Kangning Wang and Shaomin Li.  Robust distributed modal regression for massive data.  _Computational Statistics & Data Analysis_, 160:107225, 2021.
  * Xiang and Yao [2022] Sijia Xiang and Weixin Yao.  Nonparametric statistical learning based on modal regression.  _Journal of Computational and Applied Mathematics_ , 409:114130, 2022.
  * Liu et al. [2013] Jicai Liu, Riquan Zhang, Weihua Zhao, and Yazhao Lv.  A robust and efficient estimation method for single index models.  _Journal of Multivariate Analysis_ , 122:226–238, 2013\.
  * Zhang et al. [2013] Riquan Zhang, Weihua Zhao, and Jicai Liu.  Robust estimation and variable selection for semiparametric partially linear varying coefficient model based on modal regression.  _Journal of Nonparametric Statistics_ , 25(2):523–544, 2013.
  * Yang and Yang [2014] Hu Yang and Jing Yang.  A robust and efficient estimation and variable selection method for partially linear single-index models.  _Journal of Multivariate Analysis_ , 129:227–242, 2014\.
  * Zhao et al. [2014] Weihua Zhao, Riquan Zhang, Jicai Liu, and Yazhao Lv.  Robust and efficient variable selection for semiparametric partially linear varying coefficient model based on modal regression.  _Annals of the Institute of Statistical Mathematics_ , 66(1):165–191, 2014.
  * Krief [2017] Jerome M Krief.  Semi-linear mode regression.  _The Econometrics Journal_ , 20(2):149–167, 2017\.
  * Tian et al. [2017] Maoxi Tian, Jian He, and Keming Yu.  Fitting truncated mode regression model by simulated annealing.  In _Computational Optimization in Engineering-Paradigms and Applications_. Intech Open, 2017.
  * Li and Huang [2019] Xiang Li and Xianzheng Huang.  Linear mode regression with covariate measurement error.  _Canadian Journal of Statistics_ , 47(2):262–280, 2019.