# Modelling wind speed with a univariate probability distribution depending on
two baseline functions

Fábio V. J. Silveira<EMAIL_ADDRESS>Frank Gomes-Silva  Cícero C.
R. Brito  Jader S. Jale  Felipe R. S. Gusmão  Sílvio F. A. Xavier-Júnior  João
S. Rocha  Federal Institute of Education, Science and Technology of Paraíba,
João Pessoa, PB, Brazil  Department of Statistics and Informatics, Federal
Rural University of Pernambuco, Recife, PE, Brazil  Federal Institute of
Education, Science and Technology of Pernambuco, Recife, PE, Brazil
Department of Statistics, Paraíba State University, Campina Grande, PB, Brazil

###### Abstract

Characterizing the wind speed distribution properly is essential for the
satisfactory production of potential energy in wind farms, being the mixture
models usually employed in the description of such data. However, some mixture
models commonly have the undesirable property of non-identifiability. In this
work, we present an alternative distribution which is able to fit the wind
speed data adequately. The new model, called Normal-Weibull-Weibull, is
identifiable and its cumulative distribution function is written as a
composition of two baseline functions. We discuss structural properties of the
class that generates the proposed model, such as the linear representation of
the probability density function, moments and moment generating function. We
perform a Monte Carlo simulation study to investigate the behavior of the
maximum likelihood estimates of the parameters. Finally, we present
applications of the new distribution for modelling wind speed data measured in
five different cities of the Northeastern Region of Brazil.

###### keywords:

Goodness-of-fit; Identifiability; L-BFGS-B algorithm; Maximum likelihood;
Monte Carlo simulation

##  1 Introduction and proposed class

The concern about the emission of greenhouse gases and environmental
contamination from conventional energy generation procedures like coal and oil
power plants encourages research on alternative resources. A smaller impact on
the environment is an advantage of cleaner and sustainable energy production
techniques, such as solar, geothermal and wind, over the combustion of fossil
fuels.

The installed capacity of wind power in Brazil increased from 29 MW in 2005 to
roughly 16,000 MW (9% of the total capacity of electricity generation) in 2020
[1]. The suitable choice of the wind turbine must match with the wind behavior
at the site of installation. Perkin et al. [2] mention that inadequate turbine
selection results in a financially sub-optimal investment. Thus, setting the
probability distribution appropriately to model the wind speed is essential.
Eltamaly [3] used the two-parameter Weibull distribution in a new computer
program to perform the calculations required to precisely design the wind
energy system and to seek the compatibility between sites and turbines.

Ilhan and Kantar [4] declare that despite the wide acceptance of Weibull
distribution [5, 6, 7, 8], it may sometimes be poor to model all wind speed
data available in nature. Hereupon, they propose using two possible models,
namely, the skewed generalized error distribution [9] and the skewed $t$
distribution [10]. Some other distributions used for wind speed and power
modelling are Rayleigh [5], gamma [11], normal [6], generalized extreme value
[8] and Birnbaum-Saunders [12]. Additionally, applications of nonparametric
methods to wind speed modelling are also found in the literature [13, 14, 15].

Oftentimes, one requires more flexibility from the probability density
function (pdf), as in case of bimodality [16] or calm winds regime [17]. In
general, finite mixture models are more flexible than the typical single ones.
Akdag et al. [18] compared the usual biparametric Weibull and the two-
component mixture Weibull distribution in a study focused on wind regimes
presenting nearly zero percentage of null speeds; they concluded that the
mixture is more suitable to describe such wind conditions. Carta and Ramírez
[19] used three different methods to estimate the parameters of the two-
component mixture Weibull, namely, the method of moments, maximum likelihood
and least squares; they verified that there is no significant difference among
them.

The mixture density can be written as:

$\displaystyle f(x;\bm{\psi})=\sum_{i=1}^{d}w_{i}f_{i}(x;\bm{\theta}_{i})$ (1)

where the vector $\bm{\psi}=(w_{1},\ldots,w_{d-1},\bm{\eta}^{\top})^{\top}$
contains the unknown parameters of the mixture model and the vector
$\bm{\eta}$ contains all the distinct parameters in
$\bm{\theta}_{1},\ldots,\bm{\theta}_{d}$. Since $w_{1},\ldots,w_{d}$ are
positive and sum up to one, the presence of $w_{d}$ in $\bm{\psi}$ is
unnecessary.

The general definition of identifiability states that a family of densities
$\\{f(x;\bm{\psi}):\bm{\psi}\in\bm{\Psi}\\}$ is identifiable if:

$\displaystyle
f(x;\bm{\psi})=f(x;\bm{\psi}^{\star})\Leftrightarrow\bm{\psi}=\bm{\psi}^{\star}\,.$
(2)

It is not seldom that (2) fails when two or more component densities in (1)
belong to the same parametric family. Such is the case of the mixture of
normal densities. Consider $d=2$, $f_{1}$ and $f_{2}$ are normal densities,
$w_{1}=0.5$ in (1) and let
$\bm{\eta}_{1}=(\mu_{1},\sigma_{1},\mu_{2},\sigma_{2})^{\top}$ and
$\bm{\eta}_{2}=(\mu_{2},\sigma_{2},\mu_{1},\sigma_{1})^{\top}$, where
$\mu_{1}\neq\mu_{2}$, $\sigma_{1}\neq\sigma_{2}$. We have
$\bm{\psi}=(w_{1},\bm{\eta}_{1}^{\top})^{\top}\neq(w_{1},\bm{\eta}_{2}^{\top})^{\top}=\bm{\psi}^{\star}\Rightarrow
f(x;\bm{\psi})=f(x;\bm{\psi}^{\star})$. That is, (1) may be invariant under
certain permutations of the elements in the parametric vector. McLachlan and
Peel [20] mention an alternative definition of identifiability for mixture
models, such that the mixture of $d$ normal densities would be identifiable
under specific conditions. Nonetheless, they remark that it does not overcome
the complications due to the interchanging of component labels.

Models like the mixture of normal or Weibull densities are quite flexible
tools, although the parametric estimation is only credible if the distribution
under study is identifiable. We present in this paper a class, whose submodels
may be feasible alternatives to mixtures of two components from the same
parametric family. The class is derived using the method of generating classes
of probability distributions of Brito et al. [21]. Its cumulative distribution
function (cdf) is formulated as a composition of two baselines and under
certain conditions, it satisfies (2), even if both baselines belong to the
same parametric family.

The cdf of the general class is given by:

$F(x)=\zeta(x)\sum_{j=1}^{n}\int_{L_{j}(x)}^{U_{j}(x)}\mathrm{d}H(t)-\nu(x)\sum_{j=1}^{n}\int_{M_{j}(x)}^{V_{j}(x)}\mathrm{d}H(t)$
(3)

where $H$ is a cdf, $n\in\mathbb{N}$, $\zeta,\nu:\mathbb{R}\mapsto\mathbb{R}$
and $L_{j},U_{j},M_{j},V_{j}:\mathbb{R}\mapsto\mathbb{R}\cup\\{\pm\infty\\}$
are special functions that will be discussed in the next section.

###  1.1 The Normal-$(G_{1},G_{2})$ class and some structural properties

The method established by Brito et al. [21] states that if
$H,\zeta,\nu:\mathbb{R}\mapsto\mathbb{R}$ and
$L_{j},U_{j},M_{j},V_{j}:\mathbb{R}\mapsto\mathbb{R}\cup\\{\pm\infty\\}$ for
$j=1,2,3,\ldots,n$ are monotonic and right continuous functions such that:

  * (c1)

$H$ is a cdf and $\zeta$ and $\nu$ are non-negative;

  * (c2)

$\zeta(x)$, $U_{j}(x)$ and $M_{j}(x)$ are non-decreasing and $\nu(x)$,
$V_{j}(x)$, $L_{j}(x)$ are non-increasing $\forall j=1,2,3,\ldots,n$;

  * (c3)

If $\displaystyle\lim_{x\to-\infty}\zeta(x)\neq\lim_{x\to-\infty}\nu(x)$, then
$\displaystyle\lim_{x\to-\infty}\zeta(x)=0$; or
$\displaystyle\lim_{x\to-\infty}U_{j}(x)=\lim_{x\to-\infty}L_{j}(x)\,\forall
j=1,2,3,\ldots,n$, and $\displaystyle\lim_{x\to-\infty}\nu(x)=0$; or
$\displaystyle\lim_{x\to-\infty}M_{j}(x)=\lim_{x\to-\infty}V_{j}(x)\,\forall
j=1,2,3,\ldots,n$;

  * (c4)

If $\displaystyle\lim_{x\to-\infty}\zeta(x)=\lim_{x\to-\infty}\nu(x)\neq 0$,
then $\displaystyle\lim_{x\to-\infty}U_{j}(x)=\lim_{x\to-\infty}V_{j}(x)$ and
$\displaystyle\lim_{x\to-\infty}M_{j}(x)=\lim_{x\to-\infty}L_{j}(x)\,\forall
j=1,2,3,\ldots,n$;

  * (c5)

$\displaystyle\lim_{x\to-\infty}L_{j}(x)\leq\lim_{x\to-\infty}U_{j}(x)$ and if
$\displaystyle\lim_{x\to-\infty}\nu(x)\neq 0$, then
$\displaystyle\lim_{x\to+\infty}M_{j}(x)\leq\lim_{x\to+\infty}V_{j}(x)\,\forall
j=1,2,3,\ldots,n$;

  * (c6)

$\displaystyle\lim_{x\to+\infty}U_{n}(x)\geq\sup\\{x\in\mathbb{R}:H(x)<1\\}$
and
$\displaystyle\lim_{x\to+\infty}L_{1}(x)\leq\inf\\{x\in\mathbb{R}:H(x)>0\\}$;

  * (c7)

$\displaystyle\lim_{x\to+\infty}\zeta(x)=1$;

  * (c8)

$\displaystyle\lim_{x\to+\infty}\nu(x)=0$ or
$\displaystyle\lim_{x\to+\infty}M_{j}(x)=\lim_{x\to+\infty}V_{j}(x)\,\forall
j=1,2,3,\ldots,n$ and $n\geq 1$;

  * (c9)

$\displaystyle\lim_{x\to+\infty}U_{j}(x)=\lim_{x\to+\infty}L_{j+1}(x)\,\forall
j=1,2,3,\ldots,n-1$ and $n\geq 2$;

  * (c10)

$H$ is a cdf without points of discontinuity or all functions $L_{j}(x)$ and
$V_{j}(x)$ are constant at the right of the vicinity of points whose image are
points of discontinuity of $H$, being also continuous in that points.
Moreover, $H$ does not have any point of discontinuity in the set
$\displaystyle\left\\{\displaystyle\lim_{x\to\pm\infty}L_{j}(x),\lim_{x\to\pm\infty}U_{j}(x),\lim_{x\to\pm\infty}M_{j}(x),\lim_{x\to\pm\infty}V_{j}(x)\right\\}$
for some $j=1,2,3,\ldots,n$;

then Equation (3) is a cdf.

Let $n=1$, $H(t)=\mathrm{\Phi}(t)$, namely, the standard normal cdf,
$\zeta(x)=1$, $\nu(x)=0$, $U_{1}(x)=G_{1}(x)/[1-G_{1}(x)]$ and
$L_{1}(x)=\log[1-G_{2}(x)]$, where $G_{1}(x)$ and $G_{2}(x)$ are cdfs. The
function in Equation (3) turns into:

$F_{G_{1},G_{2}}(x)=\int^{\frac{G_{1}(x)}{1-G_{1}(x)}}_{\log[1-G_{2}(x)]}\mathrm{d}\mathrm{\Phi}(t).$
(4)

We took $U_{1}$ and $-L_{1}$ from the table of differentiable and
monotonically non-decreasing functions presented in the well-known paper of
Alzaatreh et al. [22], whose method was used to create generalized
distributions of the T-X family. We have intentionally picked the two simplest
functions from the cited table; alternative (and more complicated) choices for
$U_{1}$ and $L_{1}$ would naturally give rise to different classes. Defining
$M_{1}(x)$ and $V_{1}(x)$ is not relevant, since $\nu(x)=0$. Also, for obvious
reasons, there is no need to verify (c4), (c5) and (c9). The conditions (c1),
(c7), (c8) and (c10) are straightforward. As $U_{1}(x)$ and $\zeta(x)$ are
non-decreasing and $L_{1}(x)$ is non-increasing, (c2) is true. It is easy to
verify that
$\displaystyle\lim_{x\rightarrow-\infty}U_{1}(x)=0=\lim_{x\rightarrow-\infty}L_{1}(x)$;
and since $\displaystyle\lim_{x\rightarrow-\infty}\nu(x)=0$, (c3) is
satisfied. The condition (c6) is also true because
$\displaystyle\lim_{x\to+\infty}U_{1}(x)=+\infty=\sup\\{x\in\mathbb{R}:\mathrm{\Phi}(x)<1\\}$
and
$\displaystyle\lim_{x\to+\infty}L_{1}(x)=-\infty=\inf\\{x\in\mathbb{R}:\mathrm{\Phi}(x)>0\\}$.

Thereby, in agreement with the method exposed above, Equation (4) is a cdf. As
already mentioned, it can be viewed as a composite function of two baselines.
Henceforth, let it be denoted by Normal-$(G_{1},G_{2})$ class of probability
distributions.

Since $\phi(t)=\frac{1}{\sqrt{2\pi}}e^{-t^{2}/2}$, and
$\mathrm{\Phi}(x)=\int_{-\infty}^{x}\phi(t)\mathrm{d}t$, one can write
Equation (4) as follows:

$F_{G_{1},G_{2}}(x)=\mathrm{\Phi}\left(\frac{G_{1}(x)}{1-G_{1}(x)}\right)-\mathrm{\Phi}\left(\log[1-G_{2}(x)]\right)\,.$
(5)

In case of continuous $G_{1}(x)$ and $G_{2}(x)$, one can take the derivative
of Equation (5) with respect to $x$ to obtain the following pdf:

$f_{G_{1},G_{2}}(x)=\phi\left(\frac{G_{1}(x)}{1-G_{1}(x)}\right)\frac{g_{1}(x)}{[1-G_{1}(x)]^{2}}+\phi\left(\log[1-G_{2}(x)]\right)\frac{g_{2}(x)}{1-G_{2}(x)},$
(6)

where $g_{i}(x)$ is the pdf of the random variable whose cdf is $G_{i}(x)$,
for $i\in\\{1,2\\}$.

At this point, we need to define properly the support of the distributions
that emerge from the new class. Submodels of classes that may be written as a
composite function of one single baseline usually have the same support of the
baseline. However, characterizing the support of a submodel from (4) is not so
straightforward, especially if the two baselines have different supports. As
previously mentioned, given that
$U_{1}(G_{1}(x),G_{2}(x))=G_{1}(x)/[1-G_{1}(x)]$,
$L_{1}(G_{1}(x),G_{2}(x))=\log[1-G_{2}(x)]$ and $S_{H}=\mathbb{R}$, namely,
the support of $H(t)$ is $\mathbb{R}$, we have that:

  * (a)

$S_{H}$ is a convex set;

  * (b)

$U_{1}(1,1)=U_{1}(G_{1}(+\infty),G_{2}(+\infty))=+\infty=\sup\\{x\in\mathbb{R}:\mathrm{\Phi}(x)<1\\}$,
$L_{1}(1,1)=L_{1}(G_{1}(+\infty),G_{2}(+\infty))=-\infty=\inf\\{x\in\mathbb{R}:\mathrm{\Phi}(x)>0\\}$,
$U_{1}(G_{1}(x),G_{2}(x))$ and $L_{1}(G_{1}(x),G_{2}(x))$ are monotonic
functions.

According to the Theorem (T4) in [21], (a) and (b) entail that the support of
a distribution from (4) is the union of the supports of $G_{1}$ and $G_{2}$.

In the following lines, we demonstrate that, under specific conditions, the
distributions generated by (5) enjoy the attractive property of
identifiability. It is important because it assures the uniqueness of the
estimates of the parameters.

######  Theorem 1.

Let $G_{1}(x|\bm{\theta}_{1})$ and $G_{2}(x|\bm{\theta}_{2})$ be the baseline
cdfs of the normal-$(G_{1},G_{2})$ cdf $F_{G_{1},G_{2}}(x|\bm{\theta})$ (5),
$\bm{\theta}_{1}=(\theta_{1},\ldots,\theta_{r})\in\bm{\Theta}_{1}$,
$\bm{\theta}_{2}=(\theta_{r+1},\ldots,\theta_{r+m})\in\bm{\Theta}_{2}$ and
$\bm{\theta}=(\theta_{1},\ldots,\theta_{r},\theta_{r+1},\ldots,\theta_{r+m})\in\bm{\Theta}$,
where $\bm{\Theta}_{1}$, $\bm{\Theta}_{2}$ and $\bm{\Theta}$ are the
parametric spaces associated with $G_{1}$, $G_{2}$ and $F_{G_{1},G_{2}}$
respectively. If $G_{1}$ and $G_{2}$ are identifiable, then $F_{G_{1},G_{2}}$
is identifiable.

###### Proof.

Assume that
$\mathrm{\Phi}\left(\frac{G_{1}(x|\bm{\theta}_{1})}{1-G_{1}(x|\bm{\theta}_{1})}\right)=\mathrm{\Phi}\left(\frac{G_{1}(x|\bm{\theta}_{1}^{\star})}{1-G_{1}(x|\bm{\theta}_{1}^{\star})}\right)$,
where $\\{\bm{\theta}_{1},\bm{\theta}_{1}^{\star}\\}\subset\bm{\Theta}_{1}$
and $\bm{\theta}_{1}\neq\bm{\theta}_{1}^{\star}$. Since $\mathrm{\Phi}$ is
injective,
$\frac{G_{1}(x|\bm{\theta}_{1})}{1-G_{1}(x|\bm{\theta}_{1})}=\frac{G_{1}(x|\bm{\theta}_{1}^{\star})}{1-G_{1}(x|\bm{\theta}_{1}^{\star})}\Rightarrow
G_{1}(x|\bm{\theta}_{1})=G_{1}(x|\bm{\theta}_{1}^{\star})$; it is a
contradiction, because it denies the identifiability of $G_{1}$. Therefore, if
$\bm{\theta}_{1}\neq\bm{\theta}_{1}^{\star}$ then
$\mathrm{\Phi}\left(\frac{G_{1}(x|\bm{\theta}_{1})}{1-G_{1}(x|\bm{\theta}_{1})}\right)\neq\mathrm{\Phi}\left(\frac{G_{1}(x|\bm{\theta}_{1}^{\star})}{1-G_{1}(x|\bm{\theta}_{1}^{\star})}\right)$.
Analogously, it is easy to verify that for
$\\{\bm{\theta}_{2},\bm{\theta}_{2}^{\star}\\}\subset\bm{\Theta}_{2}$, if
$\bm{\theta}_{2}\neq\bm{\theta}_{2}^{\star}$ then
$\mathrm{\Phi}\left(\log[1-G_{2}(x|\bm{\theta}_{2})]\right)\neq\mathrm{\Phi}\left(\log[1-G_{2}(x|\bm{\theta}_{2}^{\star})]\right)$.

Now consider $\\{\bm{\theta},\bm{\theta}^{\star}\\}\subset\bm{\Theta}$ such
that $\bm{\theta}\neq\bm{\theta}^{\star}$ and assume that
$F_{G_{1},G_{2}}(x|\bm{\theta})=F_{G_{1},G_{2}}(x|\bm{\theta}^{\star})$. If
$\bm{\theta}_{1}=\bm{\theta}_{1}^{\star}$ and
$\bm{\theta}_{2}\neq\bm{\theta}_{2}^{\star}$, then we can infer from (5) that
$G_{2}(x|\bm{\theta}_{2})=G_{2}(x|\bm{\theta}_{2}^{\star})$, namely, an
absurd. Likewise, if $\bm{\theta}_{1}\neq\bm{\theta}_{1}^{\star}$ and
$\bm{\theta}_{2}=\bm{\theta}_{2}^{\star}$, we get to similar contradiction. If
$\bm{\theta}_{1}\neq\bm{\theta}_{1}^{\star}$ and
$\bm{\theta}_{2}\neq\bm{\theta}_{2}^{\star}$, then the assumption fails since
$F_{G_{1},G_{2}}(x|\bm{\theta})\neq F_{G_{1},G_{2}}(x|\bm{\theta}^{\star})$
for almost all values of $x$ within the support. Therefore, $F_{G_{1},G_{2}}$
is identifiable.

∎

###  1.2 Series representation

The normal cdf can be written in terms of the error function erf as follows:

$\mathrm{\Phi}(z)=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{z}{\sqrt{2}}\right)\right]\;,$
(7)

where $\mathrm{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}\mathrm{d}t$.
Since $\mathrm{erf}(z/\sqrt{2})$ may be linearly represented by:

$\displaystyle\mathrm{erf}\left(\frac{z}{\sqrt{2}}\right)$
$\displaystyle=\frac{2}{\sqrt{\pi}}\sum_{n=0}^{\infty}\frac{(-1)^{n}\cdot(z/\sqrt{2})^{2n+1}}{n!(2n+1)}$
$\displaystyle=\sqrt{\frac{2}{\pi}}\cdot\sum_{n=0}^{\infty}\left(-\frac{1}{2}\right)^{n}\frac{z^{2n+1}}{n!(2n+1)}\;,$
(8)

replacing Equation (8) in Equation (7), we get to:

$\mathrm{\Phi}(z)=\frac{1}{2}+\frac{1}{\sqrt{2\pi}}\sum_{n=0}^{\infty}\left(-\frac{1}{2}\right)^{n}\frac{z^{2n+1}}{n!(2n+1)}\;.$
(9)

Now using the result of Equation (9) in Equation (5), we have:

$\displaystyle
F_{G_{1},G_{2}}(x)=\sum_{n=0}^{\infty}\frac{(-1/2)^{n}}{n!(2n+1)\sqrt{2\pi}}\left[\underbrace{\left(\frac{G_{1}(x)}{1-G_{1}(x)}\right)^{2n+1}}_{\mathrm{A1}}-\underbrace{\left(\log[1-G_{2}(x)]\right)^{2n+1}}_{\mathrm{A2}}\right].$
(10)

A well-known result on power series raised to a positive integer $N$ states
that:

$\left(\sum_{k=0}^{\infty}a_{k}y^{k}\right)^{N}=\sum_{k=0}^{\infty}c_{k}y^{k}\;,$
(11)

where $c_{0}=a_{0}^{N}$,
$c_{k}=\frac{1}{ka_{0}}\sum_{s=1}^{k}(sN-k+s)a_{s}c_{k-s}$ for $k\geq 1$ and
$N\in\mathbb{N}$. Setting $N=2n+1$, $y=G_{1}(x)$ and $a_{k}=1$ for all $k\geq
0$, we can use the result in Equation (11) to rewrite A1 in Equation (10):

$\displaystyle\mathrm{A1}$
$\displaystyle=G_{1}(x)^{2n+1}\left(\frac{1}{1-G_{1}(x)}\right)^{2n+1}=G_{1}(x)^{2n+1}\left(\sum_{k=0}^{\infty}G_{1}(x)^{k}\right)^{2n+1}$
$\displaystyle=G_{1}(x)^{2n+1}\sum_{k=0}^{\infty}c_{1,k}\cdot
G_{1}(x)^{k}=\sum_{k=0}^{\infty}c_{1,k}\cdot G_{1}(x)^{k+2n+1}\;,$ (12)

such that $c_{1,0}=1$ and
$c_{1,k}=\frac{1}{k}\sum_{s=1}^{k}(2s[n+1]-k)c_{1,k-s}$ for $k\geq 1$.
Equation (11) also allows us to rewrite A2 in Equation 10 as follows:

$\displaystyle\mathrm{A2}$
$\displaystyle=\left(-\sum_{m=1}^{\infty}\frac{G_{2}(x)^{m}}{m}\right)^{2n+1}=-\left(\sum_{k=0}^{\infty}\frac{G_{2}(x)^{k+1}}{k+1}\right)^{2n+1}$
$\displaystyle=-G_{2}(x)^{2n+1}\left(\sum_{k=0}^{\infty}\frac{G_{2}(x)^{k}}{k+1}\right)^{2n+1}=-G_{2}(x)^{2n+1}\sum_{k=0}^{\infty}c_{2,k}\cdot
G_{2}(x)^{k}$ $\displaystyle=-\sum_{k=0}^{\infty}c_{2,k}\cdot
G_{2}(x)^{k+2n+1}$ (13)

where $c_{2,0}=1$ and
$c_{2,k}=\frac{1}{k}\sum_{s=1}^{k}\frac{2s(n+1)-k}{s+1}c_{2,k-s}$ for $k\geq
1$. Now inserting (12) and (13) in (10), we have:

$\displaystyle
F_{G_{1},G_{2}}(x)=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}\cdot
G_{i}(x)^{k+2n+1}$ (14)

where $c_{1,n,k}=\frac{(-1/2)^{n}}{n!(2n+1)\sqrt{2\pi}}c_{1,k}$ and
$c_{2,n,k}=\frac{(-1/2)^{n}}{n!(2n+1)\sqrt{2\pi}}c_{2,k}$. Using Fubini’s
theorem on differentiation we can write the derivative of (14) as follows:

$f_{G_{1},G_{2}}(x)=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}\cdot
g_{i,k+2n+1}(x)$ (15)

where $g_{i,k+2n+1}(x)=(k+2n+1)g_{i}(x)G_{i}(x)^{k+2n}$ is the pdf of a random
variable from the exponentiated family [23]. Thus, we can say that (15) is the
Normal-$(G_{1},G_{2})$ pdf (6) expressed as a linear combination of pdfs of
exponentiated distributions.

###  1.3 Raw moments, incomplete moments and moment generating function

Given that $X$ is a random variable following a distribution from the
normal-$(G_{1},G_{2})$ class, we can use (15) to write the $r$-th raw moment
of $X$ as follows:

$\displaystyle E(X^{r})$
$\displaystyle=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}\int_{-\infty}^{\infty}x^{r}g_{i,k+2n+1}(x)\mathrm{d}x$
(16)
$\displaystyle=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}E(Y_{i,k+2n+1}^{r})$
(17)

where $Y_{i,k+2n+1}$ follows the exponentiated distribution whose pdf is
$g_{i,k+2n+1}$.

Let $Q_{i}$ be the quantile function of the baseline $G_{i}$. Replacing $x$ in
(16) by $Q_{i}\left(v^{{1}/{k+2n+1}}\right)$ we can also represent (17) as:

$E(X^{r})=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}\int_{0}^{1}\left[Q_{i}\left(v^{{1}/{k+2n+1}}\right)\right]^{r}\mathrm{d}v\,.$

Similarly, one can write the $r$-th incomplete moment of $X$ as follows:

$\displaystyle m_{r}(z)$
$\displaystyle=\int_{-\infty}^{z}x^{r}f_{G_{1},G_{2}}(x)\mathrm{d}x=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}m_{r}^{\star}(z)$
$\displaystyle=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}\int_{0}^{\left[G_{i}(z)\right]^{k+2n+1}}\left[Q_{i}\left(v^{{1}/{k+2n+1}}\right)\right]^{r}\mathrm{d}v$

where $m_{r}^{\star}(z)$ is the $r$-th incomplete moment of $Y_{i,k+2n+1}$
mentioned above.

The moment generating function (mgf) of $X$ is denoted by
$M_{X}(t)=E\left(e^{tX}\right)$. It can be determined from (15) as:

$\displaystyle M_{X}(t)$
$\displaystyle=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}\int_{-\infty}^{\infty}e^{tx}g_{i,k+2n+1}(x)\mathrm{d}x$
$\displaystyle=\sum_{i=1}^{2}\sum_{n,k=0}^{\infty}c_{i,n,k}M_{Y_{k+2n+1}}(t)\,,$

where $M_{Y_{k+2n+1}}(t)$ is the mgf of $Y_{i,k+2n+1}$.

Other meaningful quantities, as the characteristic function, the probability-
weighted moments, the Rényi entropy and the order statistics can be derived
likewise by using (15).

###  1.4 Estimation and Inference

Let $\bm{X}=(x_{1},\ldots,x_{n})$ be a complete random sample of size $n$ from
the random variable $X\sim$ normal-$(G_{1},G_{2})$. Given that
$\bm{\theta}_{1}=(\theta_{1},\ldots,\theta_{r})^{\top}$ is the $r\times 1$
parametric vector associated with $G_{1}(x)=G_{1}(x|\bm{\theta}_{1})$,
$\bm{\theta}_{2}=(\theta_{r+1},\ldots,\theta_{r+m})^{\top}$ is the $m\times 1$
parametric vector associated with $G_{2}(x)=G_{2}(x|\bm{\theta}_{2})$ and
$f_{G_{1},G_{2}}(x)=f_{G_{1},G_{2}}(x|\bm{\theta})$ where
$\bm{\theta}=(\theta_{1},\ldots,\theta_{r},\theta_{r+1},\ldots,\theta_{r+m})^{\top}$,
we can write the log-likelihood function of $X$ as follows:

$\ell(\bm{\theta}|\bm{X})=\sum_{i=1}^{n}\log\left\\{\phi\left(\frac{G_{1}(x_{i})}{1-G_{1}(x_{i})}\right)\frac{g_{1}(x_{i})}{[1-G_{1}(x_{i})]^{2}}+\phi\left(\log[1-G_{2}(x_{i})]\right)\frac{g_{2}(x_{i})}{1-G_{2}(x_{i})}\right\\}\,.$

The solution of the system of equations $U(\bm{\theta}|\bm{X})=\bm{0}_{r+m}$
provides the maximum likelihood estimates (MLEs) for $\bm{\theta}$, where
$\bm{0}_{r+m}$ is an $(r+m)\times 1$ vector of zeros and
$U(\bm{\theta}|\bm{X})=\nabla_{\bm{\theta}}\ell(\bm{\theta}|\bm{X})$ is the
score vector. The elements of $U(\bm{\theta}|\bm{X})=(u_{j})_{1\leq j\leq
r+m}$ are:

$u_{j}=\sum_{i=1}^{n}\frac{1}{f_{G_{1},G_{2}}(x_{i})}\phi\left(\frac{G_{1}(x_{i})}{1-G_{1}(x_{i})}\right)\frac{1}{(1-G_{1}(x_{i}))^{2}}\left[\frac{\partial}{\partial\theta_{j}}{g_{1}}(x_{i})\right.\\\
\left.+\frac{g_{1}(x_{i})}{1-G_{1}(x_{i})}\left(2-\frac{G_{1}(x_{i})}{[1-G_{1}(x_{i})]^{2}}\right)\frac{\partial}{\partial\theta_{j}}G_{1}(x_{i})\right],\;\mathrm{for}\;1\leq
j\leq r$

and

$u_{j}=\sum_{i=1}^{n}\frac{1}{f_{G_{1},G_{2}}(x_{i})}\frac{\phi\left(\log[1-G_{2}(x_{i})]\right)}{1-G_{2}(x_{i})}\left[\frac{\partial}{\partial\theta_{j}}g_{2}(x_{i})\right.\\\
\left.+\left(1+\log[1-G_{2}(x_{i})]\right)\frac{g_{2}(x_{i})}{1-G_{2}(x_{i})}\frac{\partial}{\partial\theta_{j}}G_{2}(x_{i})\right],\;\mathrm{for}\;r<j\leq
r+m\,.$

For testing hypotheses and constructing confidence intervals for
$\bm{\theta}$, the information matrix $J(\bm{\theta}|\bm{X})$ is needed. The
expectation of $J(\bm{\theta}|\bm{X})$, denoted by
$\mathcal{I}_{\bm{\theta}}$, is the expected Fisher information matrix. Given
that certain conditions of regularity are fulfilled, the quantity
$\sqrt{n}(\widehat{\bm{\theta}}-\bm{\theta})$ follows approximately a
multivariate normal distribution
$N_{r+m}(\mathbf{0}_{r+m},\mathcal{I}_{\bm{\theta}}^{-1})$. A brings the
expression for $J(\bm{\theta}|\bm{X})$.

##  2 The proposed model

The Weibull cdf is given by $G_{W}(x|k,\lambda)=1-e^{-(x/\lambda)^{k}}$, for
$x\geq 0$, $k>0$ and $\lambda>0$. Replacing $G_{1}$ and $G_{2}$ in (5) by
$G_{W}(x|k_{1},\lambda_{1})$ and $G_{W}(x|k_{2},\lambda_{2})$ respectively, we
get to the cdf of the Normal-Weibull-Weibull distribution (NWW, for short):

$\displaystyle
F_{NWW}(x|\bm{\theta})=\mathrm{\Phi}\left(e^{(x/\lambda_{1})^{k_{1}}}-1\right)-\mathrm{\Phi}\left(-\left(\frac{x}{\lambda_{2}}\right)^{k_{2}}\right)\,,$

where $\bm{\theta}=(k_{1},\lambda_{1},k_{2},\lambda_{2})^{\top}$. The
corresponding pdf can be obtained using (6):

$f_{NWW}(x|\bm{\theta})=\phi\left(e^{(x/\lambda_{1})^{k_{1}}}-1\right)\frac{k_{1}}{\lambda_{1}}\left(\frac{x}{\lambda_{1}}\right)^{k_{1}-1}e^{(x/\lambda_{1})^{k_{1}}}+\phi\left(-\left(\frac{x}{\lambda_{2}}\right)^{k_{2}}\right)\frac{k_{2}}{\lambda_{2}}\left(\frac{x}{\lambda_{2}}\right)^{k_{2}-1}.$
(18)

Figure 1 displays some plots of the NWW pdf for different values of the
parameters. The distribution is able to fit unimodal right-skewed data (top-
left chart) and also left-skewed data (top-right chart). Notice the different
shapes of the bimodal curves in the remaining charts. For instance, in the
bottom-left chart, the vertical distance between the modes and the local
minimum in the purple curve is much greater than in the green one. We may also
notice that $\lambda_{1}$ and $\lambda_{2}$ somehow behave like shape
parameters, as in the original Weibull baselines, controlling the shape of the
“bells” (compare purple and gray curves to see the effect of varying
$\lambda_{1}$; same for blue and red curves concerning $\lambda_{2}$).

Figure 1: Plots of the Normal-Weibull-Weibull pdf.

##  3 Simulation

Performing Monte Carlo simulation studies is considerably relevant whenever
one wants to test and confirm assumptions on new statistical methods. In this
work, we want to investigate the behavior of the estimates of the parameters
of the NWW distribution under the method of maximum likelihood. For this
purpose, we used the software R version 3.4.4 [24].

Initially, we employed the Von Neumann’s acceptance-rejection method [25] to
generate pseudo-random samples from the NWW distribution; this simple method
requires only the corresponding pdf $y=f(x)$, a minorant and a majorant for
$x$ and a majorant for $y$. The procedure was replicated 10,000 times and at
each replication, four different sample sizes were considered, namely, $n=50$,
$100$, $200$ and $500$. We examined scenarios with four different values of
the parametric vector
$\bm{\theta}=(k_{1},\lambda_{1},k_{2},\lambda_{2})^{\top}$, which are
presented from the second to fifth columns of Tables 1 and 2.

For each scenario, we calculated the bias and the mean squared error (MSE) as
follows:

${\mathrm{Bias}}_{i}=\frac{1}{10000}\sum_{j=1}^{10000}\left({\widehat{\theta}}_{ij}-\theta_{i}\right)\;,\qquad{\mathrm{MSE}_{i}}=\frac{1}{10000}\sum_{j=1}^{10000}\left({\widehat{\theta}}_{ij}-\theta_{i}\right)^{2}$

where $\theta_{i}$ is the $i$-th element of $\bm{\theta}$ and
${\widehat{\theta}}_{ij}$ is the estimate for $\theta_{i}$ at the $j$-th
replication.

The global maximum of the log-likelihood function was found by using the
L-BFGS-B algorithm. It is based on the gradient projection and uses a limited
memory BFGS matrix to approximate the Hessian of the objective function [26].

Besides presenting small values, the desired behavior for both bias and MSE is
to decrease inasmuch as the sample size increases. According to Tables 1 and
2, the values of bias and MSE for all the estimated parameters are small and
the greater the sample size, the smaller the values. Thus, the results
presented in this section indicate that the MLEs of the parameters of the NWW
distribution are well-behaved when calculated using the L-BFGS-B algorithm.

Table 1: Bias of the estimates under the maximum likelihood method for the NWW model. | Actual value | Bias
---|---|---
$n$ | $k_{1}$ | $\lambda_{1}$ | $k_{2}$ | $\lambda_{2}$ | $\widehat{k}_{1}$ | $\widehat{\lambda}_{1}$ | $\widehat{k}_{2}$ | $\widehat{\lambda}_{2}$
50 | 1.3 | 2 | 1.5 | 1.8 | 0.54169 | 0.30227 | 0.42743 | 0.29174
| 3 | 1.5 | 2.8 | 2.5 | 0.55544 | 0.23582 | 0.63534 | 0.2146
| 2 | 2.2 | 6.5 | 4.1 | 0.30213 | 0.1419 | 1.42262 | 0.12017
| 1.4 | 1.6 | 4.8 | 5.1 | 0.18208 | 0.10899 | 0.86019 | 0.13781
100 | 1.3 | 2 | 1.5 | 1.8 | 0.34967 | 0.24391 | 0.31488 | 0.22752
| 3 | 1.5 | 2.8 | 2.5 | 0.43866 | 0.1199 | 0.41435 | 0.13804
| 2 | 2.2 | 6.5 | 4.1 | 0.24521 | 0.11322 | 0.89279 | 0.08923
| 1.4 | 1.6 | 4.8 | 5.1 | 0.13014 | 0.07167 | 0.54369 | 0.09648
200 | 1.3 | 2 | 1.5 | 1.8 | 0.22345 | 0.18795 | 0.24873 | 0.17187
| 3 | 1.5 | 2.8 | 2.5 | 0.28333 | 0.04346 | 0.26341 | 0.07261
| 2 | 2.2 | 6.5 | 4.1 | 0.161 | 0.06349 | 0.57998 | 0.05034
| 1.4 | 1.6 | 4.8 | 5.1 | 0.09163 | 0.04958 | 0.37906 | 0.069
500 | 1.3 | 2 | 1.5 | 1.8 | 0.13429 | 0.11752 | 0.17774 | 0.11164
| 3 | 1.5 | 2.8 | 2.5 | 0.16708 | 0.01806 | 0.17836 | 0.04992
| 2 | 2.2 | 6.5 | 4.1 | 0.09747 | 0.03622 | 0.34572 | 0.0277
| 1.4 | 1.6 | 4.8 | 5.1 | 0.06333 | 0.03001 | 0.22205 | 0.04259
Table 2: MSE of the estimates under the maximum likelihood method for the NWW model. | Actual value | MSE
---|---|---
$n$ | $k_{1}$ | $\lambda_{1}$ | $k_{2}$ | $\lambda_{2}$ | $\widehat{k}_{1}$ | $\widehat{\lambda}_{1}$ | $\widehat{k}_{2}$ | $\widehat{\lambda}_{2}$
50 | 1.3 | 2 | 1.5 | 1.8 | 0.5477 | 0.13969 | 0.37049 | 0.13065
| 3 | 1.5 | 2.8 | 2.5 | 0.43676 | 0.2146 | 0.88316 | 0.11482
| 2 | 2.2 | 6.5 | 4.1 | 0.1976 | 0.07682 | 3.60838 | 0.05433
| 1.4 | 1.6 | 4.8 | 5.1 | 0.05343 | 0.02546 | 1.32698 | 0.03216
100 | 1.3 | 2 | 1.5 | 1.8 | 0.23088 | 0.09363 | 0.17553 | 0.0809
| 3 | 1.5 | 2.8 | 2.5 | 0.27616 | 0.09641 | 0.30237 | 0.0571
| 2 | 2.2 | 6.5 | 4.1 | 0.16973 | 0.07868 | 1.39857 | 0.06496
| 1.4 | 1.6 | 4.8 | 5.1 | 0.02719 | 0.00888 | 0.49768 | 0.01486
200 | 1.3 | 2 | 1.5 | 1.8 | 0.0914 | 0.0588 | 0.10077 | 0.04583
| 3 | 1.5 | 2.8 | 2.5 | 0.12419 | 0.01817 | 0.11285 | 0.01472
| 2 | 2.2 | 6.5 | 4.1 | 0.05914 | 0.02211 | 0.57197 | 0.01905
| 1.4 | 1.6 | 4.8 | 5.1 | 0.01352 | 0.00396 | 0.23394 | 0.00751
500 | 1.3 | 2 | 1.5 | 1.8 | 0.03154 | 0.02332 | 0.05049 | 0.01887
| 3 | 1.5 | 2.8 | 2.5 | 0.04441 | 0.00052 | 0.04762 | 0.00383
| 2 | 2.2 | 6.5 | 4.1 | 0.01723 | 0.00406 | 0.19449 | 0.00315
| 1.4 | 1.6 | 4.8 | 5.1 | 0.00633 | 0.00143 | 0.07816 | 0.00286

##  4 Results and Discussion

The hourly wind speed data measured at 10 m above ground level were collected
by the National Institute of Meteorology of Brazil (INMET). The anemometers
used for measuring the wind speed (in m/s) are installed in stations located
in five cities spread in four states of the Brazilian Northeastern Region, as
illustrated in Figure 2 (blue dots indicate the geographical position of the
stations). Esperantina (denoted by Station 1) is located in the north part of
the State of Piauí. Jaguaruana (denoted by Station 2) is located in the
mesoregion of Jaguaribe in the State of Ceará. Cabaceiras (denoted by Station
3) and Monteiro (denoted by Station 4) are located in the mesoregion of
Borborema, State of Paraíba. Arapiraca (denoted by Station 5) is located in
the mesoregion of Agreste, State of Alagoas. Table 3 brings further details
about the stations and the years of wind data available.

Figure 2: Northeastern Region of Brazil and geographical position of the stations. Table 3: Details of the regions where the wind speed was measured. | Station | Latitude | Longitude | Altitude (m) | Period
---|---|---|---|---|---
1 | Esperantina |  $3^{\circ}54^{\prime}07^{\prime\prime}$S |  $42^{\circ}14^{\prime}02^{\prime\prime}$W | 59 | 2007–2018
2 | Jaguaruana |  $4^{\circ}50^{\prime}02^{\prime\prime}$S |  $37^{\circ}46^{\prime}51^{\prime\prime}$W | 20 | 2007–2018
3 | Cabaceiras |  $7^{\circ}29^{\prime}20^{\prime\prime}$S |  $36^{\circ}17^{\prime}13^{\prime\prime}$W | 382 | 2008–2018
4 | Monteiro |  $7^{\circ}53^{\prime}20^{\prime\prime}$S |  $37^{\circ}07^{\prime}12^{\prime\prime}$W | 599 | 2007–2018
5 | Arapiraca |  $9^{\circ}45^{\prime}07^{\prime\prime}$S |  $36^{\circ}39^{\prime}39^{\prime\prime}$W | 264 | 2008–2018

As we can see in Table 4, Station 4 has the highest value of the mean among
the stations in the study, whereas the highest value of variance belongs to
Station 2. Except for Station 3, whose skewness is negative, the remaining
stations have different degrees of positive skewness. On the other hand,
Station 1 has the only positive value of kurtosis and Station 5 has the lowest
one. Thus, the descriptive statistics indicate that the statistical
characteristics of the wind speed differ from station to station.

Table 4: Descriptive statistics. St. | $n$ | mean | median | min | max | variance | skewness | kurtosis
---|---|---|---|---|---|---|---|---
1 | 72637 | 1.56472 | 1.5 | 0.1 | 8.8 | 0.96082 | 0.74469 | 0.61585
2 | 71797 | 3.09401 | 3 | 0.1 | 9.4 | 2.94352 | 0.25522 | $-0.59481$
3 | 83953 | 3.21937 | 3.3 | 0.1 | 9.9 | 2.49808 | $-0.08484$ | $-0.70967$
4 | 76738 | 3.28302 | 3.3 | 0.1 | 9.6 | 2.60513 | 0.16094 | $-0.50692$
5 | 72675 | 2.87482 | 2.8 | 0.1 | 9.3 | 2.80844 | 0.17049 | $-0.92422$

We calculated the estimates of the parameters under the method of maximum
likelihood for five distributions. Besides fitting the proposed model (2), we
fitted the Normal-Normal mixture model (NN), the Weibull-Weibull mixture model
(WW), the Normal distribution (N) and the Weibull distribution (W) to each one
of the five datasets. Table 5 presents the MLEs along with the respective
standard errors in parentheses. The global maximum of the log-likelihood
function was found using the L-BFGS-B algorithm [26] for the distribution NWW,
whereas the optimization concerning the mixture models NN and WW was performed
along the lines of the EM-algorithm presented in Nguyen et al. [27]. The
standard errors are small in all scenarios, suggesting that the estimates in
Table 5 are fairly accurate for the five distributions.

Table 5: Estimates and standard errors in parentheses. Distr. | Par. | St1 | St2 | St3 | St4 | St5
---|---|---|---|---|---|---
NWW | $k_{1}$ | 1.09455 | 0.97219 | 0.93287 | 1.16262 | 0.92095
|  | (0.0066) | (0.0047) | (0.0035) | (0.0051) | (0.0041)
| $\lambda_{1}$ | 2.44439 | 4.14273 | 5.11736 | 4.61159 | 3.26931
|  | (0.0116) | (0.0300) | (0.0151) | (0.0358) | (0.0157)
| $k_{2}$ | 1.24356 | 2.21499 | 2.74024 | 2.17717 | 2.87659
|  | (0.0090) | (0.0181) | (0.0141) | (0.0232) | (0.0135)
| $\lambda_{2}$ | 2.23302 | 4.82613 | 4.47508 | 4.74835 | 4.76046
|  | (0.0104) | (0.0129) | (0.0083) | (0.0187) | (0.0059)
NN | $\mu_{1}$ | 0.99225 | 1.25346 | 1.05246 | 1.52464 | 1.09819
|  | (0.0088) | (0.0112) | (0.0135) | (0.0152) | (0.0091)
| $\sigma_{1}$ | 0.57673 | 0.73618 | 0.64260 | 0.80711 | 0.66916
|  | (0.0051) | (0.0067) | (0.0077) | (0.0082) | (0.0056)
| $\mu_{2}$ | 2.15056 | 3.69950 | 3.71423 | 3.78985 | 3.69032
|  | (0.0121) | (0.0110) | (0.0094) | (0.0127) | (0.0118)
| $\sigma_{2}$ | 0.96181 | 1.50081 | 1.28725 | 1.42131 | 1.33344
|  | (0.0038) | (0.0051) | (0.0054) | (0.0056) | (0.0063)
| $w$ | 0.50577 | 0.24754 | 0.18591 | 0.22375 | 0.31460
|  | (0.0083) | (0.0041) | (0.0036) | (0.0054) | (0.0041)
WW | $k_{1}$ | 2.46232 | 2.61226 | 3.48661 | 2.71890 | 3.63802
|  | (0.0451) | (0.0162) | (0.0183) | (0.0146) | (0.0284)
| $\lambda_{1}$ | 2.10941 | 4.09404 | 4.26754 | 4.08741 | 4.41073
|  | (0.0089) | (0.0126) | (0.0070) | (0.0103) | (0.0110)
| $k_{2}$ | 1.27757 | 1.20035 | 1.30110 | 1.26447 | 1.36100
|  | (0.0106) | (0.0088) | (0.0069) | (0.0118) | (0.0056)
| $\lambda_{2}$ | 1.42208 | 1.54955 | 1.81537 | 1.65794 | 1.86015
|  | (0.0186) | (0.0392) | (0.0242) | (0.0526) | (0.0183)
| $w$ | 0.44618 | 0.75154 | 0.71389 | 0.83225 | 0.51531
|  | (0.0183) | (0.0075) | (0.0046) | (0.0069) | (0.0056)
N | $\mu$ | 1.56471 | 3.09401 | 3.21937 | 3.28301 | 2.87482
|  | (0.0036) | (0.0064) | (0.0054) | (0.0058) | (0.0062)
| $\sigma$ | 0.98021 | 1.71565 | 1.58052 | 1.61403 | 1.67582
|  | (0.0025) | (0.0045) | (0.0038) | (0.0041) | (0.0043)
W | $k$ | 1.58959 | 1.77193 | 2.03559 | 2.07474 | 1.65091
|  | (0.0047) | (0.0054) | (0.0059) | (0.0061) | (0.0051)
| $\lambda$ | 1.73825 | 3.44821 | 3.59712 | 3.68213 | 3.18827
|  | (0.0042) | (0.0075) | (0.0063) | (0.0066) | (0.0074)

Four information criteria were used to perform comparisons among the fitted
models. Generally, such criteria indicate that the best model is the one
presenting the lowest value, since they are related to the amount of
information lost by a given model. We used the well-known Akaike information
criterion (AIC), consistent Akaike information criterion (CAIC), Bayesian
information criterion (BIC) and Hannan-Quinn information criterion (HQIC). The
statistics of Anderson-Darling (A∗) and Cramér-von Mises (W∗) [28] were also
used to compare the fitted models. Since these statistics are measures of the
difference between the empirical distribution function and the real underlying
cdf, it is reasonable to say that the smaller their values, the better the
fit. Table 6 brings the aforementioned goodness-of-fit measures for the five
cited models fitted to each station.

Table 6: Goodness-of-fit measures. Crit. | Distr. | St1 | St2 | St3 | St4 | St5
---|---|---|---|---|---|---
AIC | NWW | 189324.5 | 272976.1 | 306366.6 | 286607.1 | 268213.1
| NN | 197221.2 | 277168.5 | 310138.0 | 289011.9 | 273287.3
| WW | 189379.2 | 273075.2 | 306370.5 | 286647.3 | 268244.1
| N | 203235.1 | 281266.5 | 315112.0 | 291251.2 | 281292.0
| W | 190666.7 | 278737.7 | 319347.3 | 291232.9 | 277572.1
CAIC | NWW | 189365.2 | 273016.8 | 306408.0 | 286648.1 | 268253.8
| NN | 197272.2 | 277219.4 | 310189.7 | 289063.1 | 273338.3
| WW | 189430.2 | 273126.1 | 306422.2 | 286698.5 | 268295.1
| N | 203255.4 | 281286.9 | 315132.7 | 291271.7 | 281312.4
| W | 190687.1 | 278758.0 | 319368.0 | 291253.4 | 277592.5
BIC | NWW | 189361.2 | 273012.8 | 306404.0 | 286644.1 | 268249.8
| NN | 197267.2 | 277214.4 | 310184.7 | 289058.1 | 273333.3
| WW | 189425.2 | 273121.1 | 306417.2 | 286693.5 | 268290.1
| N | 203253.4 | 281284.9 | 315130.7 | 291269.7 | 281310.4
| W | 190685.1 | 278756.0 | 319366.0 | 291251.4 | 277590.5
HQIC | NWW | 189335.8 | 272987.4 | 306378.0 | 286618.5 | 268224.4
| NN | 197235.3 | 277182.7 | 310152.3 | 289026.1 | 273301.4
| WW | 189393.4 | 273089.3 | 306384.8 | 286661.5 | 268258.3
| N | 203240.7 | 281272.2 | 315117.7 | 291256.9 | 281297.7
| W | 190672.4 | 278743.3 | 319353.0 | 291238.6 | 277577.7
A* | NWW | 71.78 | 22.43 | 22.60 | 19.31 | 30.98
| NN | 181.13 | 51.72 | 54.33 | 27.79 | 78.72
| WW | 81.27 | 31.05 | 32.44 | 24.62 | 50.25
| N | 531.70 | 230.18 | 303.82 | 118.73 | 494.43
| W | 257.57 | 471.01 | 1117.26 | 337.07 | 784.70
W* | NWW | 7.73 | 2.33 | 2.94 | 2.58 | 3.36
| NN | 19.49 | 5.17 | 6.60 | 3.21 | 7.61
| WW | 9.35 | 3.59 | 4.22 | 3.89 | 4.82
| N | 72.08 | 31.04 | 45.14 | 16.72 | 72.83
| W | 33.87 | 62.25 | 164.54 | 46.18 | 108.58

According to Table 6, the distributions NWW and WW present the better fits
among the competing models for all stations. The four information criteria
indicate that NWW presents a higher performance over WW concerning stations 1,
2, 4 and 5. Regarding station 3, the difference between the values of AIC of
NWW and WW is not considerable, although the AIC for NWW is slightly smaller
in such case; the same behavior states for CAIC, BIC and HQIC.

Both goodness-of-fit statistics A∗ and W∗ (see last ten rows of Table 6) agree
with the information criteria in pointing NWW and WW as the two better fits.
However, they suggest that NWW outperforms WW in fitting the datasets for all
the five stations. To get more insight into these results, plots of the fitted
densities overlapping the histograms of the wind speed data for stations 1 to
5 and the corresponding cdfs are presented in Figure 3.

Figure 3: Histograms and fitted densities.

It is worth pointing out that mixture models are commonly used to fit non-
unimodal datasets, such as those represented by the histograms of stations 2,
3 and 5. Nonetheless, the results attest that the NWW accommodates such data
better than the two mixture models presented in this study. Furthermore, NWW
has one parameter less than NN or WW do.

Finally, since the NWW distribution outperforms the competing models in
fitting the wind speed data of the Northeastern Region of Brazil, according to
different information criteria and formal goodness-of-fit statistics, we have
good reasons to recommend its use to model similar data in future works. We
also encourage practitioners of statistics to investigate the modelling
benefits of the NWW (and other submodels from the Normal-($G_{1}$,$G_{2}$)
class) with respect to data describing different phenomena usually modelled by
mixtures.

##  5 Conclusions

An alternative distribution for modelling wind speed data is proposed and some
mathematical properties of the class that generates it are discussed, like the
series representation of the pdf, the moments and the moment generating
function. The general cdf of the Normal-$(G_{1},G_{2})$ class is written as a
composition of two baselines and its submodels are identifiable as long as
both baseline cdfs are. Such is the case of the NWW distribution.

The novel model has four parameters and high flexibility. It is able to fit
right-skewed and left-skewed data and its pdf presents unimodal and bimodal
shapes.

The Monte Carlo simulation studies indicate that the MLEs of the NWW
parameters behave appropriately when the optimization is performed via the
L-BFGS-B algorithm.

The modelling gains of the NWW distribution are upheld by the satisfactory
results concerning the application to the wind speed data collected in the
Northeastern Region of Brazil. The considered information criteria and the
formal goodness-of-fit statistics of Anderson-Darling and Cramér-von Mises
suggest that the proposed model outperforms other competing distributions that
are commonly employed in wind speed modelling, especially the highly
competitive mixture model of two Weibull components.

We hope that this work may encourage the investigation of the modelling
benefits of the NWW (and other identifiable submodels from the
Normal-$(G_{1},G_{2})$ class) with respect to data describing other natural
phenomena usually modelled by mixtures.

## Appendix A

The information matrix cited in section 1.4 is given by
$J(\bm{\theta}|\bm{X})=-\nabla_{\bm{\theta}}{\nabla_{\bm{\theta}}}^{\top}\ell(\bm{\theta}|\bm{X})=-(u_{jk})_{1\leq
j\leq r+m,1\leq k\leq r+m}$ where:

$\displaystyle u_{jk}=$
$\displaystyle\sum_{i=1}^{n}\frac{1}{f_{G_{1},G_{2}}(x_{i})}\phi\left(\frac{G_{1}(x_{i})}{1-G_{1}(x_{i})}\right)\frac{1}{(1-G_{1}(x_{i}))^{2}}\left\\{\left(\frac{2}{1-G_{1}(x_{i})}-\frac{G_{1}(x_{i})}{[1-G_{1}(x_{i})]^{3}}\right)\right.$
$\displaystyle\times$
$\displaystyle\left(\frac{\partial}{\partial\theta_{k}}G_{1}(x_{i})\frac{\partial}{\partial\theta_{j}}g_{1}(x_{i})+\frac{\partial}{\partial\theta_{j}}G_{1}(x_{i})\frac{\partial}{\partial\theta_{k}}g_{1}(x_{i})\right)+\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}g_{1}(x_{i})-\frac{g_{1}(x_{i})}{(1-G_{1}(x_{i}))^{2}}$
$\displaystyle\times$
$\displaystyle\left[\frac{\partial}{\partial\theta_{j}}G_{1}(x_{i})\frac{\partial}{\partial\theta_{k}}G_{1}(x_{i})\left(\frac{3{G_{1}}^{2}(x_{i})}{(1-G_{1}(x_{i}))^{2}}+\frac{4G_{1}(x_{i})}{1-G_{1}(x_{i})}+1\right)+\frac{G_{1}(x_{i})}{1-G_{1}(x_{i})}\right.$
$\displaystyle\times$
$\displaystyle\left.\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}G_{1}(x_{i})\right]+\frac{g_{1}(x_{i})G_{1}(x_{i})}{(1-G_{1}(x_{i}))^{5}}\left(\frac{\partial}{\partial\theta_{j}}G_{1}(x_{i})+\frac{\partial}{\partial\theta_{k}}G_{1}(x_{i})\right)$
$\displaystyle+$
$\displaystyle\frac{2g_{1}(x_{i})}{(1-G_{1}(x_{i}))^{2}}\left(3-\frac{2G_{1}(x_{i})}{(1-G_{1}(x_{i}))^{2}}\right)\frac{\partial}{\partial\theta_{k}}G_{1}(x_{i})\frac{\partial}{\partial\theta_{j}}G_{1}(x_{i})+\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}G_{1}(x_{i})$
$\displaystyle\times$
$\displaystyle\left.\frac{2g_{1}(x_{i})}{1-G_{1}(x_{i})}\right\\}-\sum_{i=1}^{n}\frac{1}{{f}^{2}_{G_{1},G_{2}}(x_{i})}\phi^{2}\left(\frac{G_{1}(x_{i})}{1-G_{1}(x_{i})}\right)\frac{1}{(1-G_{1}(x_{i}))^{4}}$
$\displaystyle\times$
$\displaystyle\left[\frac{\partial}{\partial\theta_{j}}g_{1}(x_{i})+\left(\frac{2g_{1}(x_{i})}{1-G_{1}(x_{i})}-\frac{g_{1}(x_{i})G_{1}(x_{i})}{(1-G_{1}(x_{i}))^{3}}\right)\frac{\partial}{\partial\theta_{j}}G_{1}(x_{i})\right]$
$\displaystyle\times$
$\displaystyle\left[\frac{\partial}{\partial\theta_{k}}g_{1}(x_{i})+\left(\frac{2g_{1}(x_{i})}{1-G_{1}(x_{i})}-\frac{g_{1}(x_{i})G_{1}(x_{i})}{(1-G_{1}(x_{i}))^{3}}\right)\frac{\partial}{\partial\theta_{k}}G_{1}(x_{i})\right],\mathrm{for}\;1\leq
j\leq r,1\leq k\leq r;$

$\displaystyle u_{jk}=$
$\displaystyle\sum_{i=1}^{n}\frac{-1}{{f}^{2}_{G_{1},G_{2}}(x_{i})}\phi\left(\frac{G_{1}(x_{i})}{1-G_{1}(x_{i})}\right)\frac{1}{(1-G_{1}(x_{i}))^{2}}\frac{\phi(\log[1-G_{2}(x_{i})])}{1-G_{2}(x_{i})}\left(\frac{\partial}{\partial\theta_{k}}g_{1}(x_{i})\right.$
$\displaystyle+$
$\displaystyle\left.\left[\frac{2g_{1}(x_{i})}{1-G_{1}(x_{i})}-\frac{g_{1}(x_{i})G_{1}(x_{i})}{(1-G_{1}(x_{i}))^{3}}\right]\frac{\partial}{\partial\theta_{k}}G_{1}(x_{i})\right)\left(\frac{\partial}{\partial\theta_{j}}g_{2}(x_{i})+(1+\log[1-G_{2}(x_{i})])\right.$
$\displaystyle\times$
$\displaystyle\left.\frac{g_{2}(x_{i})}{1-G_{2}(x_{i})}\frac{\partial}{\partial\theta_{j}}G_{2}(x_{i})\right),\mathrm{for}\;r<j\leq
r+m,1\leq k\leq r;$

$\displaystyle u_{jk}=$
$\displaystyle\sum_{i=1}^{n}\frac{-1}{{f}^{2}_{G_{1},G_{2}}(x_{i})}\phi\left(\frac{G_{1}(x_{i})}{1-G_{1}(x_{i})}\right)\frac{1}{(1-G_{1}(x_{i}))^{2}}\frac{\phi(\log[1-G_{2}(x_{i})])}{1-G_{2}(x_{i})}\left(\frac{\partial}{\partial\theta_{j}}g_{1}(x_{i})\right.$
$\displaystyle+$
$\displaystyle\left.\left[\frac{2g_{1}(x_{i})}{1-G_{1}(x_{i})}-\frac{g_{1}(x_{i})G_{1}(x_{i})}{(1-G_{1}(x_{i}))^{3}}\right]\frac{\partial}{\partial\theta_{j}}G_{1}(x_{i})\right)\left(\frac{\partial}{\partial\theta_{k}}g_{2}(x_{i})+(1+\log[1-G_{2}(x_{i})])\right.$
$\displaystyle\times$
$\displaystyle\left.\frac{g_{2}(x_{i})}{1-G_{2}(x_{i})}\frac{\partial}{\partial\theta_{k}}G_{2}(x_{i})\right),\mathrm{for}\;1\leq
j\leq r,r<k\leq r+m;$

$\displaystyle u_{jk}=$
$\displaystyle\sum_{i=1}^{n}\frac{1}{f_{G_{1},G_{2}}(x_{i})}\frac{\phi(\log[1-G_{2}(x_{i})])}{1-G_{2}(x_{i})}\left[\frac{1+\log[1-G_{2}(x_{i})]}{1-G_{2}(x_{i})}\left(\frac{\partial}{\partial\theta_{j}}G_{2}(x_{i})\frac{\partial}{\partial\theta_{k}}g_{2}(x_{i})\right.\right.$
$\displaystyle+$
$\displaystyle\left.\frac{\partial}{\partial\theta_{k}}G_{2}(x_{i})\frac{\partial}{\partial\theta_{j}}g_{2}(x_{i})\right)+\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}g_{2}(x_{i})+\frac{\partial}{\partial\theta_{j}}G_{2}(x_{i})\frac{\partial}{\partial\theta_{k}}G_{2}(x_{i})\frac{g_{2}(x_{i})}{(1-G_{2}(x_{2}))^{2}}$
$\displaystyle\times$
$\displaystyle\left(1+\log^{2}[1-G_{2}(x_{i})]+3\log[1-G_{2}(x_{i})]\right)+(1+\log[1-G_{2}(x_{i})])\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}G_{2}(x_{i})$
$\displaystyle\times$
$\displaystyle\left.\frac{g_{2}(x_{i})}{1-G_{2}(x_{i})}\right]-\sum_{i=1}^{n}\frac{1}{{f}^{2}_{G_{1},G_{2}}(x_{i})}\frac{\phi^{2}(\log[1-G_{2}(x_{i})])}{(1-G_{2}(x_{i}))^{2}}\left[\frac{(1+\log[1-G_{2}(x_{i})])g_{2}(x_{i})}{1-G_{2}(x_{i})}\right.$
$\displaystyle\times$
$\displaystyle\left.\frac{\partial}{\partial\theta_{j}}G_{2}(x_{i})+\frac{\partial}{\partial\theta_{j}}g_{2}(x_{i})\right]\left[\frac{(1+\log[1-G_{2}(x_{i})])g_{2}(x_{i})}{1-G_{2}(x_{i})}\frac{\partial}{\partial\theta_{k}}G_{2}(x_{i})+\frac{\partial}{\partial\theta_{k}}g_{2}(x_{i})\right],$
$\displaystyle\mathrm{for}\;r<j\leq r+m,r<k\leq r+m.$

## References

  * Araújo et al. [2020] J. C. H. Araújo, W. F. d. Souza, A. J. d. A. Meireles, C. Brannstrom,  Sustainability challenges of wind power deployment in coastal Ceará state, Brazil,  Sustainability 12 (2020) 5562.
  * Perkin et al. [2015] S. Perkin, D. Garrett, P. Jensson,  Optimal wind turbine selection methodology: A case-study for Búrfell, Iceland,  Renewable Energy 75 (2015) 165–172.
  * Eltamaly [2013] A. Eltamaly,  Design and implementation of wind energy system in Saudi Arabia,  Renewable Energy 60 (2013) 42–52.
  * Ilhan and Kantar [2012] U. Ilhan, Y. M. Kantar,  Analysis of some flexible families of distribution for estimation of wind speed distributions,  Applied Energy 89 (2012) 355–367.
  * Pishgar-Komleh et al. [2015] S. Pishgar-Komleh, A. Keyhani, P. Sefeedpari,  Wind speed and power density analysis based on Weibull and Rayleigh distributions (a case study: Firouzkooh county of Iran),  Renewable and Sustainable Energy Reviews 42 (2015) 313 – 322.
  * Safari [2011] B. Safari,  Modeling wind speed and wind power distributions in Rwanda,  Renewable and Sustainable Energy Reviews 15 (2011) 925–935.
  * Weisser [2003] D. Weisser,  A wind energy analysis of Grenada: an estimation using the ‘Weibull’ density function,  Renewable Energy 28 (2003) 1803 – 1812.
  * Kollu et al. [2012] R. Kollu, S. Rayapudi, S. Narasimham, K. Pakkurthi,  Mixture probability distribution functions to model wind speed distributions,  International Journal of Energy and Environmental Engineering 3 (2012).
  * Bali and Theodossiou [2008] T. G. Bali, P. Theodossiou,  Risk measurement performance of alternative distribution functions,  Journal of Risk and Insurance 75 (2008) 411–437.
  * Hansen [1994] B. E. Hansen,  Autoregressive conditional density estimation,  International Economic Review 35 (1994) 705–730.
  * Morgan et al. [2011] E. C. Morgan, M. Lackner, R. M. Vogel, L. G. Baise,  Probability distributions for offshore wind speeds,  Energy Conversion and Management 52 (2011) 15 – 26.
  * Mohammadi et al. [2017] K. Mohammadi, O. Alavi, J. Mcgowan,  Use of Birnbaum-Saunders distribution for estimating wind speed and wind power probability distributions: A review,  Energy Conversion and Management 143 (2017) 109–122.
  * Qin et al. [2011] Z. Qin, W. Li, X. Xiong,  Estimating wind speed probability distribution using kernel density method,  Electric Power Systems Research 81 (2011) 2139–2146.
  * Hu et al. [2016] B. Hu, Y. Li, H. Yang, H. Wang,  Wind speed model based on kernel density estimation and its application in reliability assessment of generating systems,  Journal of Modern Power Systems and Clean Energy 5 (2016).
  * Han et al. [2019] Q. Han, S. Ma, T. Wang, F. Chu,  Kernel density estimation model for wind speed probability distribution with applicability to wind energy assessment in China,  Renewable and Sustainable Energy Reviews 115 (2019) 109387.
  * Jaramillo and Borja [2004] O. Jaramillo, M. Borja,  Wind speed analysis in La Ventosa, Mexico: A bimodal probability distribution case,  Renewable Energy 29 (2004) 1613–1630.
  * Chang [2011] T. Chang,  Estimation of wind energy potential using different probability density functions,  Applied Energy 88 (2011) 1848–1856.
  * Akdag et al. [2010] S. Akdag, H. Bagiorgas, G. Mihalakakou,  Use of two-component Weibull mixtures in the analysis of wind speed in the Eastern Mediterranean,  Applied Energy 87 (2010) 2566–2573.
  * Carta and Ramírez [2007] J. Carta, P. Ramírez,  Analysis of two-component mixture Weibull statistics for estimation of wind speed distributions,  Renewable Energy 32 (2007) 518–531.
  * McLachlan and Peel [2000] G. McLachlan, D. Peel, Finite Mixture Models, Wiley Interscience, 2000.
  * Brito et al. [2019] C. R. Brito, L. C. Rego, W. R. Oliveira, F. Gomes-Silva,  Method for generating distributions and classes of probability distributions: the univariate case,  Hacettepe Journal of Mathematics and Statistics 48 (2019) 897–930.
  * Alzaatreh et al. [2013] A. Alzaatreh, C. Lee, F. Famoye,  A new method for generating families of continuous distributions,  Metron 71 (2013) 63–79.
  * Mudholkar and Srivastava [1993] G. S. Mudholkar, D. K. Srivastava,  Exponentiated Weibull family for analyzing bathtub failure-rate data,  IEEE transactions on reliability 42 (1993) 299–302.
  * R Core Team [2018] R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2018\. URL: http://www.R-project.org/.
  * Von Neumann [1951] J. Von Neumann, Various techniques used in connection with random digits, Applied Mathematics Series 12, National Bureau of Standards, Washington, DC, USA, 1951.
  * Byrd et al. [1995] R. H. Byrd, P. Lu, J. Nocedal, C. Zhu,  A limited memory algorithm for bound constrained optimization,  SIAM Journal on Scientific Computing 16 (1995) 1190–1208.
  * Nguyen et al. [2018] H. D. Nguyen, D. Wang, G. J. McLachlan,  Randomized mixture models for probability density approximation and estimation,  Information Sciences 467 (2018) 135 – 148.
  * Chen and Balakrishnan [1995] G. Chen, N. Balakrishnan,  A general purpose approximate goodness-of-fit test,  Journal of Quality Technology 27 (1995) 154–161.


# The Cepheid Distance to the Narrow-Line Seyfert 1 Galaxy NGC 4051

W. Yuan Department of Physics & Astronomy, Johns Hopkins University,
Baltimore, MD 21218, USA  L. M. Macri George P. and Cynthia W. Mitchell
Institute for Fundamental Physics and Astronomy, Department of Physics and
Astronomy,
Texas A&M University, College Station, TX 77843, USA  B. M. Peterson
Department of Astronomy, The Ohio State University, 140 W 18th Ave, Columbus,
OH 43210, USA  Center for Cosmology and AstroParticle Physics, The Ohio State
University, 191 West Woodruff Ave, Columbus, OH 43210, USA  Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  A. G.
Riess Department of Physics & Astronomy, Johns Hopkins University, Baltimore,
MD 21218, USA  Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  M. M. Fausnaugh Department of Astronomy, The Ohio
State University, 140 W 18th Ave, Columbus, OH 43210, USA  Kavli Institute for
Space and Astrophysics Research, Massachusetts Institute of Technology,
77 Massachusetts Avenue, Cambridge, MA 02139, USA  S. L. Hoffmann Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
G. S. Anand Institute for Astronomy, University of Hawaii, 2680 Woodlawn
Drive, Honolulu, HI 96822, USA  M. C. Bentz Department of Physics and
Astronomy, Georgia State University, 25 Park Place, Suite 605, Atlanta, GA
30303, USA  E. Dalla Bontà Dipartimento di Fisica e Astronomia “G. Galilei,”
Università di Padova, Vicolo dell’Osservatorio 3, I-35122 Padova, Italy  INAF-
Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5 I-35122,
Padova, Italy  R. I. Davies Max Planck Institut für extraterrestrische Physik,
Postfach 1312, D-85741, Garching, Germany  G. De Rosa Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  L. Ferrarese NRC
Herzberg Astronomy and Astrophysics, National Research Council, 5071 West
Saanich Road, Victoria, BC V9E 2E7, Canada  C. J. Grier Department of
Astronomy, The Ohio State University, 140 W 18th Ave, Columbus, OH 43210, USA
Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson,
AZ 85721, USA  E. K. S. Hicks Department of Physics and Astronomy, University
of Alaska Anchorage, AK 99508, USA  C. A. Onken Research School of Astronomy
and Astrophysics, Australian National University, Canberra, ACT 2611,
Australia  Australian Research Council (ARC) Centre of Excellence in All-sky
Astrophysics (CAASTRO)  R. W. Pogge Department of Astronomy, The Ohio State
University, 140 W 18th Ave, Columbus, OH 43210, USA  Center for Cosmology and
AstroParticle Physics, The Ohio State University, 191 West Woodruff Ave,
Columbus, OH 43210, USA  T. Storchi-Bergmann Departamento de Astronomia,
Universidade Federal do Rio Grande do Sul, Av. Bento Goncalves 9500, 91501
Porto Alegre, RS, Brazil  M. Vestergaard DARK Niels Bohr Institute, University
of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark  Steward Observatory,
University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

###### Abstract

We derive a distance of $D=16.6\pm 0.3$ Mpc ($\mu=31.10\pm 0.04$ mag) to the
archetypal narrow-line Seyfert 1 galaxy NGC 4051 based on Cepheid
Period–Luminosity relations and new Hubble Space Telescope multiband imaging.
We identify 419 Cepheid candidates and estimate the distance at both optical
and near-infrared wavelengths using subsamples of precisely-measured
variables. We compare our independent photometric procedures and distance-
estimation methods to those used by the SH0ES team and find agreement to 0.01
mag. The distance we obtain suggests an Eddington ratio $\dot{m}\approx 0.2$
for NGC 4051, typical of narrow-line Seyfert 1 galaxies, unlike the seemingly-
odd value implied by previous distance estimates. We derive a peculiar
velocity of $-490\pm 34$ km s-1 for NGC 4051, consistent with the overall
motion of the Ursa Major Cluster in which it resides. We also revisit the
energetics of the NGC 4051 nucleus, including its outflow and mass accretion
rates.

##  1 Introduction

The distance to an astronomical object is one of the most fundamental
measurements that can be made, and yet it is also one of the most difficult.
Without an accurate distance, measurements of many other parameters are of
limited utility. In the case of the nearest active galactic nuclei (AGNs), the
lack of accurate and precise distance measurements to these objects precludes
a robust calibration of their physical properties, even though, compared to
the general AGN population, their detailed structure is best resolved and
their dynamics and energetics are the most easily studied. Thus, obtaining
accurate distances is the most urgent priority toward achieving an accurate
physical description of the accretion, feedback, and energetics of nearby
AGNs.

NGC 4051, one of the original six Seyfert galaxies (Seyfert, 1943) and one of
the lowest-redshift AGNs ($z=0.00234$), illustrates some of the difficulties.
NGC 4051 is the archetype of the subclass of AGNs known as “narrow-line
Seyfert 1 galaxies” (NLS1s; Osterbrock & Pogge, 1985), which are thought to be
high Eddington ratio objects; the Eddington ratio $\dot{m}$ is defined as the
ratio of mass accretion rate to the Eddington accretion rate,
$\dot{m}=\dot{M}/\dot{M}_{\mathrm{Edd}}$. Since the accretion rate is
$\dot{M}=L_{\mathrm{bol}}/\eta c^{2}$ where $L_{\mathrm{bol}}$ is the
bolometric luminosity and $\eta$ is the efficiency, as long as $\eta$ is
independent of $\dot{M}$, we can also write $\dot{m}=L_{\rm
bol}/L_{\mathrm{Edd}}$ where $L_{\mathrm{Edd}}$ is the Eddington luminosity.
To compute the bolometric luminosity, we must know the distance.
Unfortunately, redshift-independent distances for NGC 4051 are ambiguous.
Distances quoted in the literature range from 9 to 18 Mpc. If NGC 4051 is at
the lower end of this range, which is close to the redshift-based distance of
9.5 Mpc, its Eddington ratio is surprisingly low for an NLS1 galaxy. By using
the bolometric correction and weighted mean mass and luminosity from Dalla
Bontà et al. (2020), we find that at a distance of 9 Mpc the Eddington ratio
would be $\dot{m}\approx 0.07$, which is fairly typical of other local Seyfert
galaxies. However, if NGC 4051 is at the high end of this distance range, at
18 Mpc, $\dot{m}\approx 0.28$, which is more typical of NLS1s.

NGC 4051 is the first AGN for which strong evidence for inflows in molecular
gas was found using near-infrared integral field spectroscopy (Riffel et al.,
2008). The inferred geometry of the inflow depends strongly on the distance to
the galaxy, taken to be 9.3 Mpc in that work. Similarly, X-ray spectra reveal
a fast ionized outflow that loses kinetic energy by interaction with the
interstellar medium in the host galaxy (Pounds, 2014); the energetics and
timescales again strongly depend on the distance to the AGN and are
consequently systematically uncertain by a factor of several. Moreover,
spatially resolved observations of the narrow-line region reveal a strong
outflow component (Fischer et al., 2013), and again the energetics are
uncertain on account of the poorly constrained distance. The quantitative
results of these observations are important because they are probes of the
feedback mechanisms through which supermassive black hole (SMBH) accretion can
modify the interstellar medium and regulate star formation in the host galaxy.
Understanding the physical mechanisms at work in both feeding and feedback,
and determining whether or not the feedback mechanisms of theoretical models
actually work in nature, requires detailed observations of the complex inner
structure of AGNs. Nearby galaxies afford the only opportunity to test in
detail — i.e., on sub-hundred pc scales — the prescriptions used in models of
galaxy and SMBH co-evolution. It is on scales of less than 100 pc that
morphological differences between active and quiescent galaxies appear (Simões
Lopes et al., 2007; García-Burillo & Combes, 2012; Hicks et al., 2013; Davies
et al., 2014). There is observational and theoretical evidence that disk
processes are important in feeding AGNs across broad ranges of redshift and
luminosity (e.g., Hopkins & Quataert, 2010; Schawinski et al., 2012; Kocevski
et al., 2012; Riffel et al., 2013; Schnorr-Müller et al., 2014; Storchi-
Bergmann & Schnorr-Müller, 2019), so what we learn about local AGNs is
applicable to galaxies at $z>1$, where co-evolution largely occurs.

Finally, it is only in the nearest AGNs that we can compare SMBH masses
measured from stellar or gas dynamical modeling, which are distance-dependent,
with those measured by reverberation mapping (RM), which are distance-
independent. NGC 4051 is challenging as the SMBH radius of influence ($r_{\rm
BH}=GM/\sigma_{*}^{2}$, where $M$ is the BH mass and $\sigma_{*}$ is the
stellar bulge velocity dispersion) is $\sim 0.67$ pc; whether or not it is
spatially resolvable with either the Hubble or the James Webb  Space
Telescopes (HST, JWST) is unclear given the factor of 2–3 uncertainty in the
SMBH mass and the factor of two uncertainty in the distance, but it should
nevertheless be possible to get at least meaningful constraints on the central
mass with the IFU on JWST.

Our goal is to eliminate distance-dependent ambiguities by determining a
Cepheid-based distance to NGC 4051. The Cepheid Period–Luminosity Relation
(PLR) or Leavitt Law (Leavitt & Pickering, 1912) has proven to be the most
accurate and reliable empirical standard candle with a reach of $D\approx 50$
Mpc. It has been extensively applied to dozens of extragalactic systems,
including Bentz et al. (2019) and Yuan et al. (2020, hereafter Y20) who
recently measured the Cepheid-based distances to the Seyfert 1 galaxies NGC
6814 and NGC 4151, respectively.

In this work, we derive a Cepheid distance to NGC 4051 using multiband HST
observations. The rest of this paper is organized as follows. We present the
observations, data reduction, and photometry in §2 and our Cepheid selection
methods in §3. Our main result, the Cepheid distance to NGC 4051, is given in
§4 and we discuss its implications in §5.

##  2 Observations, Data Reduction, and Photometry

We observed NGC 4051 and carried out the subsequent analysis following the
same strategy used by Y20 for NGC 4151. We obtained optical and near-infrared
(NIR) HST images of a field that covers a large area of the disk (see Figure
1). We briefly summarize the observations, data reduction, and photometry
procedure below and refer interested readers to Y20 for further details.

Figure 1: Pseudo-color image of NGC 4051 using our master frames in F350LP,
F555W, F814W, and F160W as the luminance layer and the blue, green, and red
channels, respectively. A square-root scaling was applied to all layers to
increase the visibility. Color is only shown where all layers overlap, while
an inverse grayscale is used in areas not imaged with F160W. The compass lines
are $20\arcsec$ in length. The black diamond-shaped area near the center is
due to the gap between the WFC3/UVIS CCDs.

The NGC 4051 observations include 12, 3, 3, and 6 epochs in F350LP, F555W,
F814W, and F160W, respectively. The observation baseline spans 73 days, during
which the orientation angle of HST changed once. The observation log and
exposure times are given in Table 1. We note that one of the two F555W images
obtained in epoch 9 had a significantly worse FWHM, possibly due to increased
jitter, and thus was excluded from further analysis.

Table 1: Observation Log

Epoch | MJD | Dither$\times$Exposure Time (seconds)
---|---|---
|  | F350LP | F555W | F814W | F160W
1 | 58089.2 | 3$\times$350 |  $\cdots$ |  $\cdots$ | 2$\times$553
2 | 58107.8 | 3$\times$350 | 2$\times$550 |  $\cdots$ |  $\cdots$
3 | 58111.9 | 3$\times$350 |  $\cdots$ | 2$\times$550 |  $\cdots$
4 | 58117.4 | 3$\times$350 |  $\cdots$ |  $\cdots$ | 2$\times$553
5 | 58122.6 | 3$\times$350 | 2$\times$550 |  $\cdots$ |  $\cdots$
6 | 58129.0 | 3$\times$350 |  $\cdots$ |  $\cdots$ | 2$\times$553
7 | 58135.6 | 3$\times$350 |  $\cdots$ | 2$\times$550 |  $\cdots$
8 | 58141.0 | 3$\times$350 |  $\cdots$ |  $\cdots$ | 2$\times$553
9 | 58145.3 | 3$\times$350 | 2$\times$550a |  $\cdots$ |  $\cdots$
10 | 58149.6 | 3$\times$350 |  $\cdots$ |  $\cdots$ | 2$\times$553
11 | 58155.2 | 3$\times$350 |  $\cdots$ | 2$\times$550 |  $\cdots$
12 | 58162.1 | 3$\times$350 |  $\cdots$ |  $\cdots$ | 2$\times$553

Note. — $a$: The first image of this set was excluded due to its poor FWHM,
possibly due to increased jitter.

Table 2: Secondary Standards

R.A. | Dec. | F350LP | F555W | F814W
---|---|---|---|---
[J2000] | [mag(mmag)]
180.76803 | 44.53942 | 23.520(13) | 23.604(18) | 22.862(16)
180.77227 | 44.53125 | 23.943(15) | 24.202(21) | 23.034(14)
180.77507 | 44.52808 | 23.601(23) | 23.707(19) | 23.372(30)
180.77511 | 44.53187 | 23.566(17) | 23.772(20) | 22.798(20)
180.77630 | 44.53949 | 23.089(21) | 23.228(20) | 23.009(30)
180.77693 | 44.54960 | 23.822(19) | 24.041(23) | 22.981(21)
180.78067 | 44.54180 | 23.060(18) | 23.309(23) | 22.898(24)
180.78152 | 44.54246 | 22.434(11) | 22.603(13) | 21.713(11)
180.78185 | 44.51941 | 22.824(20) | 22.894(15) | 22.462(28)
180.78258 | 44.54943 | 23.771(18) | 23.959(24) | 23.325(31)

Note. — Coordinates are based on the WCS solution of the first F350LP image.
Uncertainties are given in parentheses and are expressed in units of mmag.
This table is available in its entirety in machine-readable form.

We registered and drizzle-combined all images using the AstroDrizzle v2.2.6
package. We created master frames in each band by combining the respective
images obtained in all epochs. We also created drizzled images of each F350LP
epoch for variable detection. We modeled the surface brightness gradients
across all images and subtracted them using the same methods as described by
Y20.

We carried out point-spread function (PSF) photometry on the NGC 4051 images
following identical procedures to those used in the case of NGC 4151 (Y20). We
firstly derived a source list through a two-pass source detection on the
F350LP master frame using DAOPHOT/ALLSTAR (Stetson, 1987), then performed
ALLFRAME (Stetson, 1994) time-series photometry on the 12 epochs of images
obtained through that filter. The frame-to-frame magnitude offsets were
calculated based on 69 carefully-selected secondary standards which are listed
in Table 2. We transformed the coordinates of the objects in the F350LP source
list to carry out fixed-position PSF photometry in the F555W, F814W & F160W
master frames (no time-series measurements were obtained, only mean
magnitudes). We performed an additional detection step on the F160W master
frame after removing the stars that were present in the F350LP source list.

##  3 Cepheid Identification and Characterization

We applied the same Cepheid selection procedures as described by Y20 and
identified 419 Cepheids in NGC 4051. We measured their periods, amplitudes,
and crowding bias corrections.

###  3.1 Cepheid identification

Our Cepheid identification included four steps: variability index cut
(Stetson, 1996, adjusted for magnitude dependence), template fitting, period
and amplitude cuts, and visual inspection. We firstly selected 4519 variable
candidates with Stetson index $L>0.65$. This threshold was slightly higher
than the initial choice for NGC 4151 ($L>0.5$) for two reasons: (1) the
minimum $L$ value for NGC 4151 Cepheids that were eventually used for distance
determination was $\sim$0.9; (2) the NGC 4051 field contains more Cepheid-
yielding star-forming regions, and thus the number of variables is not a
limiting factor for the NGC 4051 distance determination. The $L$–magnitude
distribution for NGC 4051 and our updated $L$ threshold are shown in Figure 2.
Secondly, we fit the Yoachim et al. (2009) Cepheid templates to the F350LP
light curves of these variable candidates using the method described by Y20.
Thirdly, we restricted the sample based on the best-fit periods and amplitudes
as done in Y20; this yielded 2804 objects with $P<90$d and peak-to-trough
amplitudes between 0.4 and 1.7 mag. Finally, we visually inspected the F350LP
light curves of these 2804 objects and identified 419 Cepheids. Figure 3 shows
6 examples to illustrate the light curve quality of the selected variables.
The light curves of all 419 objects are given in Table 3.

Figure 2: Adjusted $L$ index against F350LP magnitudes for all the sources.
The red line indicates our threshold for variable selection. Magenta points
indicate the 419 Cepheids described in §3.1. Figure 3: Example F350LP light
curves of selected Cepheids spanning the period range. The blue curves
indicate the best-fit Cepheid templates. We plot two cycles of pulsation for
visualization purposes. Table 3: F350LP light curves

Cepheid | ta | mag | $\sigma$
---|---|---|---
ID | [day] | [mag]
9734 | 0.729 | 26.518 | 0.062
9734 | 19.282 | 26.473 | 0.070
9734 | 23.394 | 26.610 | 0.055
9734 | 28.883 | 26.705 | 0.079
9734 | 34.113 | 26.199 | 0.045
9734 | 40.468 | 26.267 | 0.067
9734 | 47.088 | 26.703 | 0.111
9734 | 52.458 | 26.353 | 0.053
9734 | 56.823 | 26.135 | 0.051
9734 | 61.067 | 26.201 | 0.051
9734 | 66.699 | 26.537 | 0.051
9734 | 73.602 | 26.595 | 0.078

Note. — $a$: JD$-$2458089\. This table is available in its entirety in
machine-readable form.

###  3.2 Crowding correction and magnitude calibration

We applied crowding, aperture, and phase corrections to the Cepheid
magnitudes, then calibrated our photometric measurements on the Vega system.

While we adopted the same crowding correction method as for NGC 4151, we
slightly updated our procedure in F160W relative to Y20. In our previous work,
we only measured the F160W magnitude and crowding bias for Cepheids with
$P>25$ d. We removed such a restriction for this analysis; 379 of the
optically discovered Cepheids were located in the smaller NIR field of view.
Of these, 278 were detected and photometered. We determined their F160W
crowding correction $\Delta\mathrm{{\it F160W}}$ by carrying out three passes
of the artificial star test procedure, compared to only one pass for NGC 4151.
In the first pass, we determined the artificial star magnitudes using a PLR
based on Cepheids with $P>40$ d. In the next two passes, we based our
calculation on Cepheids with $\sigma(\Delta\mathrm{{\it F160W}})<0.25$ mag as
determined from the preceding pass. The restriction of $P>40$ d and
$\sigma(\Delta\mathrm{{\it F160W}})<0.25$ mag yielded a more accurate
estimation of the input artificial star magnitudes, and our tests converged
faster.

Figure 4: Color–magnitude diagram of optically-selected Cepheids (filled
symbols). Representative errorbars are shown on the right. Gray dots are a
subsample of field stars. The area between the red lines shows the Cepheid
instability strip for $P\\!>\\!10$ d, derived from LMC observations (Riess et
al., 2019). The arrow indicates the effect of $A_{V}=1$ mag. Figure 5: F555W
(top), F814W (middle), and $W^{I}$ (bottom) PLRs for optically selected
Cepheids (filled circles). Rejected variables are indicated by open circles.
The solid and dashed blue lines indicate the best-fit PLRs and $\pm 1\sigma$
scatter, respectively. Table 4: Cepheid Properties

ID | R.A. | Dec. | $P$ | Amp. | F350LP | F555W | F814W | F160W | $SB^{c}$ | Fd |  $\Delta$F555W |  $\Delta$F814W |  $\Delta$F160W
---|---|---|---|---|---|---|---|---|---|---|---|---|---
| [J2000]a | [day] | [mag] | [mag (mmag)]b |  |  | [mmag]e
82977 | 180.769100 | 44.538239 | 24.407 | 0.938 | 26.236(103) | 26.742(132) | 25.362(150) | 24.062(601) | 308.0 |  | -7(106) | 3(120) | 518(587)
88413 | 180.769205 | 44.537613 | 13.390 | 1.157 | 26.617(133) | 26.843(166) | 26.142(239) | 26.284(573) | 295.3 |  | 57(145) | 59(205) | 1308(508)
142230 | 180.769992 | 44.531406 | 14.785 | 0.885 | 26.257(158) | 26.506(181) | 25.572(184) |  $\cdots$ | 253.9 |  | 33(148) | -3(165) |  $\cdots$
129685 | 180.771219 | 44.533581 | 14.021 | 0.479 | 26.440(113) | 26.759(145) | 25.610(178) | 25.161(638) | 251.5 |  | 13(127) | 28(152) | 926(614)
125827 | 180.771822 | 44.534374 | 76.796 | 0.436 | 25.147( 36) | 25.584( 69) | 24.152( 46) | 22.963(222) | 275.9 | H | 7( 37) | 9( 36) | 5(212)

Note. — $a$: Coordinates from first F350LP image. $b$: Fully-calibrated Vega-
system magnitudes, including crowding corrections and their uncertainties.
$c$: F160W-band local surface brightness in units of counts/s/sq arcsec. $d$:
Cepheids used for $W^{I}$ PLR (I), for both $W^{I}$ and $W^{H}$ PLRs (H), or
rejected from optical PLRs (o). $e$: Crowding corrections and uncertainties.
This table is available in its entirety in machine-readable form.

Unlike Y20, we derived aperture corrections using field stars (listed in Table
2) for F350LP, F555W, and F814W. In the case of F160W, we adopted the aperture
correction value derived by Y20 due to the insufficient number of suitable
stars in the respective master frame. All magnitudes presented in this paper
are in the Vega system and were calibrated using the current STScI zeropoints
(see §4.2 of Y20 for the detailed method). The adopted Vega zeropoints for
“infinite aperture” magnitudes are 26.817, 25.843, 24.712, and 24.689 for
F350LP, F555W, F814W, and F160W, respectively.

We present the Cepheid periods, amplitudes, fully calibrated magnitudes, and
crowding corrections in Table 4. We adopted the Riess et al. (2019) PLR slopes
of -2.76, -2.96, and -3.31 in F555W, F814W, and $W^{I}$ (defined below),
respectively. We iteratively fit these PLRs to our sample, rejecting the
single largest $>$3$\sigma$ outlier in each step until convergence. We arrived
at a final optically-selected sample of 397 variables. Their color–magnitude
diagram and optical PLRs are shown in Figure 4 and Figure 5, respectively.

##  4 Results

We determined a distance to NGC 4051 of $D=16.62\pm 0.32$ Mpc, equivalent to a
distance modulus of $\mu=31.103\pm 0.042$ mag, by applying the reddening-free
Wesenheit indices (Madore, 1982; Riess et al., 2019)

$\displaystyle W^{I}$ $\displaystyle=$ $\displaystyle\mathrm{{\it
F814W}}-1.3(\mathrm{{\it F555W}}-\mathrm{{\it F814W}}),$ (1) $\displaystyle
W^{H}$ $\displaystyle=$ $\displaystyle\mathrm{{\it F160W}}-0.386(\mathrm{{\it
F555W}}-\mathrm{{\it F814W}})$ (2)

to Cepheid subsamples with precise measurements. In this section, we present
the details of this selection, the corresponding distance determination, and a
comparison to results obtained with an alternative procedure (Riess et al.,
2016).

###  4.1 Selection of subsamples for distance determination

We selected high-quality subsamples of Cepheids for distance determination
purposes, in order to obtain tight Cepheid PLRs and reduce the noise from
poorly-measured outliers that were not detected in §3.2. The leading source of
scatter in the $W^{I}$ index is the uncertainty in the F555W crowding
correction $\sigma(\Delta\mathrm{{\it F555W}})$. Thus, we empirically
determined a cut of $\sigma(\Delta\mathrm{{\it F555W}})<0.1$ mag based on the
weighted scatter of the $W^{I}$ PLR residuals, as shown in the upper panel of
Figure 6. We applied an iterative 3$\sigma$ outlier rejection to obtain a
subsample of 123 Cepheids for the optical distance determination. Similarly,
we then applied a further empirical cut of $\sigma(\Delta\mathrm{{\it
F160W}})<0.3$ mag to arrive at a “NIR distance subsample” of 47 Cepheids, as
shown in the lower panel of Figure 6. The weighted scatter of the full $W^{I}$
and $W^{H}$ samples are 0.35 and 0.41 mag, respectively, while the
corresponding values for the distance-determination subsamples are 0.23 and
0.24 mag, respectively.

Figure 6: Optical (upper) and NIR (lower) subsamples used for distance
determination, selected through limits on crowding correction uncertainty and
iterative 3$\sigma$ outlier rejection. The selected Cepheids, rejected
outliers, and excluded Cepheids are indicated by black filled circles, red
filled circles, and gray open circles, respectively. The blue curves show the
weighted scatter of PLR residuals as a function of maximum crowding correction
uncertainty, while the green vertical lines indicate our adopted cuts. Figure
7: 3-D plot showing the correlation between the pairs among
$\sigma(\Delta\mathrm{{\it F160W}})$, period, and local surface brightness
(SB). The NIR subsample used for distance determination is indicated by red
points. An interactive version of this plot is available in the online
article. Figure 8: Optical (left) and NIR (right) grayscale images of NGC 4051
showing the locations of Cepheids from various subsamples. Blue, white and red
circles denote Cepheids used for distance determination, excluded variables,
and rejected objects, respectively (matching the black, red and open symbols
used in Fig. 6). Magenta crosses in the right panel indicate optically-
selected Cepheids that were not detected in the NIR frame. Figure 9: Optical
(left) and NIR (right) Wesenheit PLRs for the Cepheid subsamples used for
distance determination. The solid lines indicate the best-fit PLRs, while the
dashed lines indicate the 1$\sigma$ standard deviations.

These selection criteria are different from those adopted by Y20, where we
applied cuts on period and F160W surface brightness. Indeed, the crowding
correction uncertainty is strongly correlated with period (or luminosity) and
local surface brightness, as shown in Figure 7. We had noted in §4.3 of Y20
that a cut on $\sigma(\Delta\mathrm{{\it F160W}})$ yielded results consistent
with those obtained via a cut on local surface brightness. Fig. 8 shows the
positions of all Cepheids, color-coded by subsample. We expect that a long-
period or low-crowding Cepheid subsample will be minimally impacted by
systematics.

###  4.2 The distance to NGC 4051

We present the optical and NIR reddening-free Wesenheit PLRs for the distance-
determination Cepheid subsamples in Figure 9. The PLR fit results are
summarized in Table 5. The LMC-relative distance moduli are based on the
Cepheid PLRs reported by Riess et al. (2019) and include a correction of
$+0.021$ mag to our F160W photometry to account for the difference in the
adopted Vega-based zeropoint. The $W^{I}$ and $W^{H}$ PLRs yielded
$\Delta\mu_{I}=12.602\pm 0.040$ and $\Delta\mu_{H}=12.672\pm 0.056$ mag,
respectively. These two measurements are consistent within 1$\sigma$, and
their combination yields $\Delta\mu=12.626\pm 0.033$ mag. We adopt the LMC
distance modulus of $\mu_{\mathrm{LMC}}=18.477\pm 0.026$ mag (Pietrzyński et
al., 2019) to arrive at a distance modulus for NGC 4051 of $\mu=31.103\pm
0.042$ mag, equivalent to a distance of $D=16.62\pm 0.32$ Mpc.

We present our error budget in Table 6 using the same format as for NGC 4151
in Y20. Compared to our earlier work, we adopted a smaller and conservative
magnitude calibration error of 0.01 mag for the optical PLR because we derived
the F814W and F555W aperture corrections using field stars. The phase
correction errors were also reduced as a result of the increased number of
Cepheids in NGC 4051. The other entries are identical to NGC 4151. We
estimated the uncertainty due to possible metallicity effects by using the
[O/H] dispersion in 21 Cepheid hosts and the terms reported by Riess et al.
(2016). We refer interested readers to §4.4 of Y20 for a detailed discussion
of those systematics. We did not find a distance dependence on the minimum
period cut for either the optical or NIR Cepheid subsamples, as shown in
Figure 10.

Figure 10: Distance dependence on minimum period cut for the optical (upper)
and NIR (lower) Cepheid subsamples.

We compared our distance determination to literature values, as compiled by
the NASA/IPAC Extragalactic Database (NED). While the uncertainties in most of
those previous determinations are too large to provide a statistically
meaningful comparison against our result, Yoshii et al. (2014) derived a
distance with comparable precision as ours. They used the AGN time-lag method
and dust reverberation modeling to obtain a distance of $D=17.6\pm 0.6$ Mpc.
Their result is consistent with ours within 1.4$\sigma$, though we note that
their distance determination for NGC 4151 is not consistent with ours, as
reported in Y20.

Table 5: Wesenheit PLRs used for distance determination

Index | Expression | Slope | Intercept | Weighted | $N$
---|---|---|---|---|---
|  |  |  | LMC | N 4051 | scatter |
$W^{I}$ | ${\it F814W}-1.3$ | $({\it F555W}-{\it F814W})$ | -3.31 | 15.933 | 28.535 | 0.23 | 123
$W^{H}$ | ${\it F160W}-0.386$ | $({\it F555W}-{\it F814W})$ | -3.26 | 15.915a | 28.587b | 0.24 | 47

Note. — $a$: Corrected for the F160W-band count-rate nonlinearity; see Y20 for
details. $b$: Adjusted by 0.021 mag for the difference in adopted Vega
zeropoints between Riess et al. (2019) and this study.

Table 6: Error Budget Source | $W^{I}$ | $W^{H}$
---|---|---
| (mag) | (mag)
PLR scatter | 0.021 | 0.035
PLR slope | $\cdots$ | 0.033
Magnitude calibration | 0.01 | 0.02
Metallicity | 0.03 | 0.02
Phase correction | 0.009 | 0.004
Reddening law | 0.01 | $\cdots$
Subtotal | 0.040 | 0.056
Subtotal combined | 0.033
LMC distance | 0.026
Total | 0.048 | 0.062
Total combined | 0.042
Figure 11: Top: NIR Cepheid PLRs using SH0ES F160W photometry and SH0ES
selection criteria. The resultant distance modulus and weighted PLR scatter
are shown in the top-left of the panel. Middle: Same as top but for our NIR
calibration sample. Bottom: Error distribution of the SH0ES artificial star
tests. For each artificial star, we subtracted the mean bias correction and
divided by the standard deviation of all artificial stars associated with the
same Cepheid.

###  4.3 Comparison to SH0ES procedure

The SH0ES project (Riess et al., 2016) measures Cepheid F160W magnitudes
following a different procedure, which gives us an opportunity to compare the
two approaches. We carried out SH0ES photometry on the NGC 4051 F160W master
image and derived NIR PLRs for a sample of Cepheids that were selected using
the SH0ES criteria. Both methods use fixed-position PSF photometry in F160W,
with positions derived from the higher signal-to-noise ratio and better-
sampled optical master frame. Likewise, they both use artificial star tests to
correct for crowding. Unlike this study, where we fixed positions for all
sources appearing in both F350LP and F160W images, the SH0ES approach only
fixes the positions for the Cepheids.

In addition, the SH0ES procedure includes two features to derive accurate
photometry and Gaussian errors in magnitude space, as demonstrated in
artificial star tests (see the bottom panel of Figure 11): (1) a large
“protection” radius of $0.75\times$ FWHM to avoid spurious de-blends; (2) the
selection of artificial stars that match the apparent source displacement of
the Cepheid relative to the predicted optical position (a secondary
consequence of blending, Riess et al., 2009). As a result, the SH0ES approach
yielded unbiased photometry for a broader range of local crowding levels and
thus resulted in a larger sample of Cepheids for distance determination. As
shown in the top panel of Figure 11, the SH0ES criteria yielded 182 Cepheids
and a distance modulus of 31.16 mag. We also determined a distance modulus
using the SH0ES-based photometry for 45 Cepheids that are in common with our
NIR subsample, which yielded 31.17 mag (middle panel of Figure 11). Both are
in excellent agreement (to within $0.01$ mag) with our results. We expect a
statistical uncertainty between methods mostly due to the finite number of
unique samples of artificial stars, $\sim$0.006 mag, and a systematic
uncertainty arising from the choice of different PSFs, sky background regions,
fitting techniques, etc., $\sim$0.02 mag.

##  5 Discussion

###  5.1 Eddington ratio in NGC 4051

Peterson et al. (2004) and Collin et al. (2006) recognized that there might be
a problem with the redshift-based distance to NGC 4051, and it was these
concerns that led us to this investigation. Peterson et al. (2004, their
Figure 16) and Collin et al. (2006, their Figure 6) showed that the Eddington
ratio in NGC 4051 was surprisingly low for an NLS1 galaxy, one of only two
such objects with $\dot{m}<0.1$ (the other is MCG-06-30-15; Bentz et al.,
2016). Now, at the corrected distance of 16.6 Mpc, the Eddington ratio is
$\dot{m}\approx 0.2$, a value more typical of NLS1 galaxies. In Figure 12, we
update the mass–luminosity diagram of Peterson et al. (2004), where we see
that the improved and increased distance places NGC 4051 in the same locus as
other NLS1s, with $\dot{m}\gtrsim 0.1$, instead of the seemingly-odd value
implied by previous distance estimates.

Figure 12: Black hole mass versus bolometric luminosity for AGNs with
reverberation-based masses using the H$\beta$ emission line. Shown in blue are
sources gathered from the literature by Dalla Bontà et al. (2020, their Table
A1); shown in green are SDSS-RM AGNs with H$\beta$-based reverberation masses
(Grier et al., 2017) as compiled by Dalla Bontà et al. (2020, their Table A2).
NLS1s are shown as open circles. The updated measurement for NGC 4051 is
represented by the open star in the lower left of the diagram.

###  5.2 Distance and motion

NGC 4051 is a member of what is often referred to as the Ursa Major Cluster
(Tully et al., 1996). CosmicFlows-3 (Tully et al., 2016) reports an overall
distance to Ursa Major of $17.2\pm 1.0$ Mpc, which agrees with our Cepheid
distance to NGC 4051 within the quoted uncertainty. Over the years, it has
become clear that the Ursa Major complex is likely composed of several
distinct bound groups (Karachentsev et al., 2013; Kourkchi & Tully, 2017).
Both Karachentsev et al. (2013) and Kourkchi & Tully (2017) place NGC 4051
into a group with the lenticular galaxy NGC 4111 as the brightest member.
Kourkchi & Tully (2017) list 14 members in the group, and have redshift-
independent distances for five of these galaxies. With this information, they
derive a weighted average group distance of 14.36 $\pm$ 1.01 Mpc, which is
broadly consistent with our derived Cepheid distance to NGC 4051 when taking
into account the size of the group itself.

The observed velocity of NGC 4051 in the reference frame of the Local Group is
$v_{\mathrm{LG}}=740\pm 3$ km s-1 (Karachentsev & Makarov, 1996). Adopting a
value of the Hubble constant of $74.0\pm 1.4$ km s-1 Mpc-1 (Riess et al.,
2019) allows us to calculate a peculiar velocity of $v_{\mathrm{pec}}=-490\pm
34$ km s-1. While seemingly high, we note that the vast majority of galaxies
in the Ursa Majoris complex exhibit a negative peculiar velocity, with an
average peculiar velocity with respect to the Local Group of
$v_{\mathrm{pec}}=-337$ km s-1 (Karachentsev et al., 2013). Similarly, Pak et
al. (2014) identified 166 Ursa Major members and found that they lie in a
redshift range of 620–1208 km s-1. Considering that all these galaxies are
roughly equidistant, the peculiar velocity of NGC 4051 lies well within the
inferred peculiar velocity range of other Ursa Major galaxies. When placed in
this context, our derived distance to NGC 4051 is consistent with
expectations.

###  5.3 Revisiting the energetics of the nucleus

It seems that the most commonly assumed distance for NGC 4051 is the redshift-
based distance of 9.3 Mpc, which is considerably smaller than the value
reported here. This has a significant effect on linear distances within the
galaxy, which now become a factor of 1.8 larger, and luminosity measures,
which increase by a factor of 3.2.

Barbosa et al. (2009), for example, estimated the mass outflow rate and power
along with three other Seyfert galaxies. Repeating their outflow calculation
and increasing the diameter of the eastern outflow from 45 pc to 80 pc results
in an outflow rate $\dot{M}=2.2\,n_{100}f\,\mbox{$M_{\odot}$}\,{\rm yr}^{-1}$,
where $n_{100}$ is the electron density in units of $100\,{\rm cm}^{-3}$ and
$f$ is the volume filling factor, for which Barbosa et al. (2009) suggest a
value of $10^{-3}$. We similarly correct the diameter of the western outflow
from 44.5 pc to 79.2 pc, thus increasing the outflow rate to
$1.4\,n_{100}f\,\mbox{$M_{\odot}$}\,{\rm yr}^{-1}$. The total outflow rate is
the sum of these, $\dot{M}=3.6\,n_{100}f\,\mbox{$M_{\odot}$}\,{\rm yr}^{-1}$.
Similarly, recalculation of the kinetic luminosity gives

$P=\frac{1}{2}\dot{M}V^{2}=2.6\times 10^{40}n_{\rm
100}f\,\mbox{erg\,s${}^{-1}$},$ (3)

where $V$ is the outflow velocity.

Barbosa et al. (2009) compare the mass outflow rate to the mass accretion rate
required to sustain the luminosity, but unfortunately using a bolometric
luminosity that is much too low. The luminosity they use traces back to a
Chandra X-ray observation obtained in 2000 April, when NGC 4051 was in an
extremely faint X-ray state (Collinge et al., 2001); NGC 4051 has some history
of the X-rays nearly turning off for long periods of time (Uttley et al.,
1999; Peterson et al., 2000). We instead use the bolometric luminosity from
Dalla Bontà et al. (2020), which adjusted to our distance of 16.6 Mpc becomes
$\log L_{\rm bol}=43.38$. It is worth noting at this point that the kinetic
luminosity of the outflow is some six orders of magnitude less than the
radiative luminosity for $f\sim 10^{-3}$.

The corrected mass accretion rate is thus

$\dot{M}_{\rm acc}=\frac{L_{\rm
bol}}{\eta\,c^{2}}=0.0042\,\mbox{$M_{\odot}$}\,{\rm yr}^{-1},$ (4)

where we have taken the efficiency to have the nominal value $\eta=0.1.$
Barbosa et al. (2009) argue that if the total outflow rate exceeds the
accretion rate, then the outflow must be comprised of accelerated
circumnuclear material. In our revised calculation, this condition is met if
$f\gtrsim 10^{-3}n_{100}^{-1}$.

On the other hand, the analysis of high-velocity outflows in NGC 4051 by
Pounds (2014) does not need reassessment as a distance of 15.2 Mpc, based on a
Tully–Fisher measurement, was assumed (Pounds & Vaughan, 2011a, b; Pounds &
King, 2013). Adjustment to a distance of 16.6 Mpc results in little change.

##  6 Summary

We analyzed HST time-series observations for the Seyfert 1 galaxy NGC 4051 and
identified 419 Cepheids in this system. Using subsamples of Cepheids that were
precisely measured at optical and NIR wavelengths, we derived a distance of
$D=16.6\pm 0.3$ Mpc, equivalent to a distance modulus of $\mu=31.10\pm 0.04$
mag. We calculated the Eddington ratio of the AGN residing in this galaxy to
be $\dot{m}\approx 0.2$, thus confirming NGC 4051 as a typical NLS1 in terms
of the AGN mass–luminosity ratio. We obtained a peculiar velocity of $-490\pm
34$ km s-1 for NGC 4051 and compared it to those of the other Ursa Major
member galaxies. We revisited the energetics of the NGC 4051 nucleus based on
Barbosa et al. (2009) and derived an outflow rate of $\sim 3.6\times
10^{-3}\,\mbox{$M_{\odot}$}\,{\rm yr}^{-1}$ and an outflow kinetic luminosity
of $\sim 2.6\times 10^{37}\,{\rm erg}\,{\rm s}^{-1}$ if $n_{100}f\sim 10^{-3}$
was assumed, as well as a mass accretion rate of $\sim 4.2\times
10^{-3}\,\mbox{$M_{\odot}$}\,{\rm yr}^{-1}$.

We thank K.Z. Stanek for helpful advice early in this project. MMF and BMP
thank Tod Lauer for an enlightening conversation. Support for HST program
GO-13765 was provided by NASA through a grant from the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS5-26555. MCB gratefully
acknowledges support from the NSF through grant AST-2009230. MV gratefully
acknowledges support from the Independent Research Fund Denmark via grant
number DFF 8021-00130. This research has made use of the NASA/IPAC
Extragalactic Database (NED), which is funded by the National Aeronautics and
Space Administration and operated by the California Institute of Technology.

## References

  * Barbosa et al. (2009) Barbosa, F. K. B., Storchi-Bergmann, T., Cid Fernandes, R., Winge, C., & Schmitt, H. 2009, MNRAS, 396, 2, doi: 10.1111/j.1365-2966.2009.14485.x
  * Bentz et al. (2016) Bentz, M. C., Cackett, E. M., Crenshaw, D. M., et al. 2016, ApJ, 830, 136, doi: 10.3847/0004-637X/830/2/136
  * Bentz et al. (2019) Bentz, M. C., Ferrarese, L., Onken, C. A., Peterson, B. M., & Valluri, M. 2019, ApJ, 885, 161, doi: 10.3847/1538-4357/ab48fb
  * Collin et al. (2006) Collin, S., Kawaguchi, T., Peterson, B. M., & Vestergaard, M. 2006, A&A, 456, 75, doi: 10.1051/0004-6361:20064878
  * Collinge et al. (2001) Collinge, M. J., Brandt, W. N., Kaspi, S., et al. 2001, ApJ, 557, 2, doi: 10.1086/321635
  * Dalla Bontà et al. (2020) Dalla Bontà, E., Peterson, B. M., Bentz, M. C., et al. 2020, arXiv e-prints, arXiv:2007.02963.  https://arxiv.org/abs/2007.02963
  * Dalla Bontà et al. (2020) Dalla Bontà, E., Peterson, B. M., Bentz, M. C., et al. 2020, ApJ, 903, 112, doi: 10.3847/1538-4357/abbc1c
  * Davies et al. (2014) Davies, R. I., Maciejewski, W., Hicks, E. K. S., et al. 2014, ApJ, 792, 101, doi: 10.1088/0004-637X/792/2/101
  * Fischer et al. (2013) Fischer, T. C., Crenshaw, D. M., Kraemer, S. B., & Schmitt, H. R. 2013, The Astrophysical Journal Supplement Series, 209, 1, doi: 10.1088/0067-0049/209/1/1
  * García-Burillo & Combes (2012) García-Burillo, S., & Combes, F. 2012, in Journal of Physics Conference Series, Vol. 372, Journal of Physics Conference Series, 012050, doi: 10.1088/1742-6596/372/1/012050
  * Grier et al. (2017) Grier, C. J., Trump, J. R., Shen, Y., et al. 2017, ApJ, 851, 21, doi: 10.3847/1538-4357/aa98dc
  * Hicks et al. (2013) Hicks, E. K. S., Davies, R. I., Maciejewski, W., et al. 2013, ApJ, 768, 107, doi: 10.1088/0004-637X/768/2/107
  * Hopkins & Quataert (2010) Hopkins, P. F., & Quataert, E. 2010, MNRAS, 405, L41, doi: 10.1111/j.1745-3933.2010.00855.x
  * Karachentsev & Makarov (1996) Karachentsev, I. D., & Makarov, D. A. 1996, AJ, 111, 794, doi: 10.1086/117825
  * Karachentsev et al. (2013) Karachentsev, I. D., Nasonova, O. G., & Courtois, H. M. 2013, MNRAS, 429, 2264, doi: 10.1093/mnras/sts494
  * Kocevski et al. (2012) Kocevski, D. D., Faber, S. M., Mozena, M., et al. 2012, ApJ, 744, 148, doi: 10.1088/0004-637X/744/2/148
  * Kourkchi & Tully (2017) Kourkchi, E., & Tully, R. B. 2017, ApJ, 843, 16, doi: 10.3847/1538-4357/aa76db
  * Leavitt & Pickering (1912) Leavitt, H. S., & Pickering, E. C. 1912, Harvard College Observatory Circular, 173, 1
  * Madore (1982) Madore, B. F. 1982, ApJ, 253, 575, doi: 10.1086/159659
  * Osterbrock & Pogge (1985) Osterbrock, D. E., & Pogge, R. W. 1985, ApJ, 297, 166, doi: 10.1086/163513
  * Pak et al. (2014) Pak, M., Rey, S.-C., Lisker, T., et al. 2014, MNRAS, 445, 630, doi: 10.1093/mnras/stu1722
  * Peterson et al. (2000) Peterson, B. M., McHardy, I. M., Wilkes, B. J., et al. 2000, ApJ, 542, 161, doi: 10.1086/309518
  * Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682, doi: 10.1086/423269
  * Pietrzyński et al. (2019) Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200, doi: 10.1038/s41586-019-0999-4
  * Pounds (2014) Pounds, K. 2014, Space Sci. Rev., 183, 339, doi: 10.1007/s11214-013-0008-4
  * Pounds & King (2013) Pounds, K. A., & King, A. R. 2013, MNRAS, 433, 1369, doi: 10.1093/mnras/stt807
  * Pounds & Vaughan (2011a) Pounds, K. A., & Vaughan, S. 2011a, MNRAS, 413, 1251, doi: 10.1111/j.1365-2966.2011.18211.x
  * Pounds & Vaughan (2011b) —. 2011b, MNRAS, 415, 2379, doi: 10.1111/j.1365-2966.2011.18866.x
  * Riess et al. (2019) Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85, doi: 10.3847/1538-4357/ab1422
  * Riess et al. (2009) Riess, A. G., Macri, L., Casertano, S., et al. 2009, ApJ, 699, 539, doi: 10.1088/0004-637X/699/1/539
  * Riess et al. (2016) Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56, doi: 10.3847/0004-637X/826/1/56
  * Riffel et al. (2008) Riffel, R., Pastoriza, M. G., Rodríguez-Ardila, A., & Maraston, C. 2008, MNRAS, 388, 803, doi: 10.1111/j.1365-2966.2008.13440.x
  * Riffel et al. (2013) Riffel, R. A., Storchi-Bergmann, T., & Winge, C. 2013, MNRAS, 430, 2249, doi: 10.1093/mnras/stt045
  * Schawinski et al. (2012) Schawinski, K., Simmons, B. D., Urry, C. M., Treister, E., & Glikman, E. 2012, MNRAS, 425, L61, doi: 10.1111/j.1745-3933.2012.01302.x
  * Schnorr-Müller et al. (2014) Schnorr-Müller, A., Storchi-Bergmann, T., Nagar, N. M., & Ferrari, F. 2014, MNRAS, 438, 3322, doi: 10.1093/mnras/stt2440
  * Seyfert (1943) Seyfert, C. K. 1943, ApJ, 97, 28, doi: 10.1086/144488
  * Simões Lopes et al. (2007) Simões Lopes, R. D., Storchi-Bergmann, T., de Fátima Saraiva, M., & Martini, P. 2007, ApJ, 655, 718, doi: 10.1086/510064
  * Stetson (1987) Stetson, P. B. 1987, Publications of the Astronomical Society of the Pacific, 99, 191, doi: 10.1086/131977
  * Stetson (1994) —. 1994, Publications of the Astronomical Society of the Pacific, 106, 250, doi: 10.1086/133378
  * Stetson (1996) —. 1996, Publications of the Astronomical Society of the Pacific, 108, 851, doi: 10.1086/133808
  * Storchi-Bergmann & Schnorr-Müller (2019) Storchi-Bergmann, T., & Schnorr-Müller, A. 2019, Nature Astronomy, 3, 48, doi: 10.1038/s41550-018-0611-0
  * Tully et al. (2016) Tully, R. B., Courtois, H. M., & Sorce, J. G. 2016, AJ, 152, 50, doi: 10.3847/0004-6256/152/2/50
  * Tully et al. (1996) Tully, R. B., Verheijen, M. A. W., Pierce, M. J., Huang, J.-S., & Wainscoat, R. J. 1996, AJ, 112, 2471, doi: 10.1086/118196
  * Uttley et al. (1999) Uttley, P., McHardy, I. M., Papadakis, I. E., Guainazzi, M., & Fruscione, A. 1999, MNRAS, 307, L6, doi: 10.1046/j.1365-8711.1999.02801.x
  * Yoachim et al. (2009) Yoachim, P., McCommas, L. P., Dalcanton, J. J., & Williams, B. F. 2009, AJ, 137, 4697, doi: 10.1088/0004-6256/137/6/4697
  * Yoshii et al. (2014) Yoshii, Y., Kobayashi, Y., Minezaki, T., Koshida, S., & Peterson, B. A. 2014, ApJ, 784, L11, doi: 10.1088/2041-8205/784/1/L11
  * Yuan et al. (2020) Yuan, W., Fausnaugh, M. M., Hoffmann, S. L., et al. 2020, ApJ, 902, 26, doi: 10.3847/1538-4357/abb377