The upper-crossing/solution (US) algorithm for root-finding with strongly
stable convergence

Xunjian LI† and Guo-Liang TIAN†,∗

Department of Statistics and Data Science, Southern University of Science and
Technology,

Shenzhen 518055, Guangdong Province, P. R. China Email:
<EMAIL_ADDRESS>

†Equally contributed

∗Corresponding author’s Email<EMAIL_ADDRESS>
Abstract. In this paper, we propose a new and broadly applicable root–finding
method, called as the upper-crossing/solution (US) algorithm, which belongs to
the category of non-bracketing (or open domain) methods. The US algorithm is a
general principle for iteratively seeking the unique root $\theta^{*}$ of a
non-linear equation $g(\theta)=0$ and its each iteration consists of two
steps: an upper-crossing step (U-step) and a solution step (S-step), where the
U-step finds an upper-crossing function or a $U$-function
$U(\theta|\theta^{(t)})$ [whose form depends on $\theta^{(t)}$ being the
$t$-th iteration of $\theta^{*}$] based on a new notion of so-called changing
direction inequality, and the S-step solves the simple $U$-equation
$U(\theta|\theta^{(t)})=0$ to obtain its explicit solution $\theta^{(t+1)}$.
The US algorithm holds two major advantages: (i) It strongly stably converges
to the root $\theta^{*}$; and (ii) it does not depend on any initial values,
in contrast to Newton’s method. The key step for applying the US algorithm is
to construct one simple $U$-function $U(\theta|\theta^{(t)})$ such that an
explicit solution to the $U$-equation $U(\theta|\theta^{(t)})=0$ is available.
Based on the first–, second– and third–derivative of $g(\theta)$, three
methods are given for constructing such $U$-functions. We show various
applications of the US algorithm in such as calculating quantile in continuous
distributions, calculating exact $p$-values for skew null distributions, and
finding maximum likelihood estimates of parameters in a class of
continuous/discrete distributions. The analysis of the convergence rate of the
US algorithm and some numerical experiments are also provided. Especially,
because of the property of strongly stable convergence, the US algorithm could
be one of the powerful tools for solving an equation with multiple roots.

Keywords: Changing direction inequality; Multiple roots; Non-bracketing
methods; Strongly stable convergence; US algorithm.

##  1 . Introduction

Solving roots of an equation or zeros of a function is the most fundamental
problem in computational mathematics. Many problems in science and engineering
can be formulated as: Given a continuous function $g(x)$ of a single variable
defined in the real line ${\mathbb{R}}$, find the value $x^{*}\in{\mathbb{R}}$
such that $g(x^{*})=0$. These problems are usually called root finding
problems. Calculating the maximum likelihood estimate (MLE) of a parameter in
classical statistics and computing the posterior mode in Bayesian statistics
can be viewed as finding a solution to a non-linear equation. Let $g(\theta)$
be a non-linear function of a single variable and there exist a unique root
$\theta^{*}$ to the following equation:

$\displaystyle g(\theta)=0,\quad\theta\in\Theta\subseteq{\mathbb{R}},$ (1.1)

where $\Theta$ is the parameter space or an interval in the real line
${\mathbb{R}}$.

There exist many methods for seeking the root of (1.1) in numerical analysis
and they are roughly categorized into two classes:

  * (i)

Closed domain (bracketing) methods, for example, the bisection method and the
false position method (Wood 1992, Chabert 1999 p.83–112).

  * (ii)

Open domain (non-bracketing) methods such as Newton’s method, the secant
method, and Muller’s method (Wolfe 1959, Boyd & Vandenberghe 2004 p.484–496,
Costabile et al. 2006).

The bisection method (Wood 1992) can guarantee to obtain a solution since the
root is trapped in the closed interval, while it converges slowly since the
information about actual functions behavior is not utilized. By using
information about the function $g(\theta)$, the false position method (Chabert
1999) generally converges more rapidly than the bisection method, however it
does not give a bound on the error of the root. Non-bracketing methods are not
as robust as bracketing methods and may actually diverge. However, they use
information about the function itself to refine the approximations of the
root. Thus, they are considerably more efficient than bracketing ones.

It is well known that Newton’s method or the Newton–Raphson (NR) algorithm is
a golden standard for seeking the root of (1.1) with the following NR
iteration:

$\displaystyle\theta^{(t+1)}=\theta^{(t)}-\frac{g(\theta^{(t)})}{g^{\prime}(\theta^{(t)})},$
(1.2)

where $\theta^{(t)}$ denotes the $t$-th iteration/approximation of the root
$\theta^{*}$ and $g^{\prime}(\theta)$ is the first-order derivative of
$g(\theta)$. Although the NR algorithm converges quadratically, it is
sensitive to initial values (Lindstrom & Bates 1988). This property hinders
its applications in repeated use of the NR algorithm such as in bootstrapping
calculation of the standard deviation of an MLE. One advantage of Newton’s
method is that it is an optimization technique as well as an root–finding
technique. In computational statistics, we know that the
expectation–maximization (EM) algorithm (Dempster et al. 1977) and
minorization–maximization (MM) algorithm (Lange et al. 2000, Hunter & Lange
2004) are two popular optimization tools for calculating MLEs of parameters
with monotonic convergence, rather than two root–finding tools.

To overcome the drawback of sensitiveness to initial values associated with
the NR algorithm and to retain a similar property to the monotonic convergence
associated with the EM/MM algorithms, we in this paper propose a new
root–finding method, called as the upper-crossing/solution (US) algorithm,
which belongs to the category of open domain methods. The US algorithm is a
general principle for iteratively seeking the root $\theta^{*}$ of equation
(1.1) with strongly stable convergence (its definition will be given in the
beginning of Section 2), and its each iteration consists of an upper-crossing
step (U-step) and a solution step (S-step), where

  * U-step:

Find a $U$-function $U(\theta|\theta^{(t)})$ satisfying (2.5);

  * S-step:

Solve the $U$-equation to obtain its root $\theta^{(t+1)}$ as showed by (2.6).

The two-step process is repeated until convergence occurs.

The rest of the paper is organized as follows. In Section 2, first we define
two new relation symbols on the changing direction inequality, and the concept
of weakly/strongly stable convergence. Next we introduce new notions such as
upper-crossing function or $U$-function and $U$-equation. Finally each US
iteration with two steps (i.e., U-step and S-step) is formulated and the
strongly stable convergence of the US algorithm is proved. Three methods for
constructing $U$-functions are given in Section 3. Section 4 presents various
applications such as calculating quantile in continuous distributions,
calculating exact $p$-values for skew null distributions, and finding MLEs of
parameters in a class of continuous/discrete distributions. Section 5 provides
the analysis of the convergence rate of the US algorithm. Numerical
experiments and some discussions are given in Sections 6–7. More technical
details are put in three Appendices. Here we first present one simple example
to give the flavour of the US algorithm.

Example 1 (Root of an equation involving a trigonometric function). Suppose
that we want to find the root of the equation $g(x)\triangleq\cos(0.5\pi
x)-x=0$ for all $x\in{\mathbb{R}}$. From (3.1), note that

$g^{\prime}(x)=-0.5\pi\sin(\pi x/2)-1\geqslant-0.5\pi-1\triangleq b_{1},$

then the U-step of the US iteration is to find an upper-crossing function or a
$U$-function $U(x|x^{(t)})$, in this example, it is given by (3.6); i.e.,

$U(x|x^{(t)})=g(x^{(t)})+b_{1}(x-x^{(t)})=g(x^{(t)})-(0.5\pi+1)x+(0.5\pi+1)x^{(t)}.$

And the S-step is to solve the $U$-equation: $U(x|x^{(t)})=0$, and to obtain
its root $x^{(t+1)}$ as

$\displaystyle
x^{(t+1)}=x^{(t)}-\frac{g(x^{(t)})}{b_{1}}=\frac{(0.5\pi+1)x^{(t)}+g(x^{(t)})}{0.5\pi+1}.$
(1.3)

Figure 1 plots both $g(x)$ and $U(x|x^{(t)})$ with the root $x^{*}=0.594612$.
Using the initial value $x^{(0)}=-1$, Figure 1(a) shows that the sequence
$\\{x^{(t)}\\}_{t=0}^{\infty}$ generated by (1.3) strongly stably converges to
$x^{*}$ satisfying $x^{(0)}<x^{(1)}<\cdots<x^{(t)}<\cdots\leqslant x^{*}$;
while Figure 1(b) shows that with $x^{(0)}=2$, the sequence
$\\{x^{(t)}\\}_{t=0}^{\infty}$ generated by (1.3) strongly stably converges to
$x^{*}$ satisfying $x^{*}\leqslant\cdots<x^{(t)}<\cdots<x^{(1)}<x^{(0)}$.

Table 1 displayed these numerical results of the US iterations for two
different initial values $x^{(0)}=-1,2$. With the rate of convergence being
0.1069, the US algorithm moved for 10 and 8 steps, respectively. $\|$

(a) (b)

Figure 1. Plots of $g(x)=\cos(0.5\pi x)-x$ and
$U(x|x^{(t)})=-(0.5\pi+1)x+(0.5\pi+1)x^{(t)}+g(x^{(t)})$. The true root of the
equation $g(x)=0$ is $x^{*}=0.594612$. (a) With the initial value
$x^{(0)}=-1$, the sequence $\\{x^{(t)}\\}_{t=0}^{\infty}$ generated by (1.3)
strongly stably converges to $x^{*}$ satisfying
$x^{(0)}<x^{(1)}<\cdots<x^{(t)}<\cdots\leqslant x^{*}$; (b) With the initial
value $x^{(0)}=2$, the sequence $\\{x^{(t)}\\}_{t=0}^{\infty}$ generated by
(1.3) strongly stably converges to $x^{*}$ satisfying
$x^{*}\leqslant\cdots<x^{(t)}<\cdots<x^{(1)}<x^{(0)}$.

Table 1. The US algorithm for solving the equation: $\cos(0.5\pi x)-x=0$ with
two different initial values $x^{(0)}=-1,2$

$t$ | $x^{(t)}$ | $g(x^{(t)})$ | $\varepsilon^{(t)}$ | Rate | $x^{(t)}$ | $g(x^{(t)})$ | $\varepsilon^{(t)}$ | Rate
---|---|---|---|---|---|---|---|---
0 |  $-$1 | 1 |  $-$1.59461 | .7561 | 2 |  $-$3 | 1.40539 | .1697
1 |  $-$.611015 | 1.18471 |  $-$1.20563 | .6178 | .833046 |  $-$.573792 | .238435 | .0639
2 |  $-$.150180 | 1.12248 |  $-$.744791 | .4138 | .609850 |  $-$.034652 | .015239 | .1154
3 |  .286449 | .614018 |  $-$.308163 | .2249 | .596371 |  $-$.003983 | .001759 | .1192
4 |  .525293 | .153170 |  $-$.069319 | .1405 | .594821 |  $-$.000475 | .000210 | .1197
5 |  .584874 | .021967 |  $-$.009738 | .1225 | .594637 |  $-$.000057 | .000025 | .1197
6 |  .593418 | .002700 |  $-$.001193 | .1201 | .594615 |  $-$.000007 | .000003 | .1196
7 |  .594468 | .000324 |  $-$.000143 | .1198 | .594612 |  $-$.000001 | .000000 | .1182
8 |  .594594 | .000039 |  $-$.000017 | .1196 | .594612 |  .000000 | .000000 | .1069
9 |  .594610 | .000005 |  $-$.000002 | .1182 | – | – | – | –
10 |  .594611 | .000001 |  .000000 | .1069 | – | – | – | –

Note: $\varepsilon^{(t)}\triangleq x^{(t)}-x^{*}$ is the difference between
the current iteration and the root; Rate
$\triangleq|\varepsilon^{(t+1)}/\varepsilon^{(t)}|$ is the rate of
convergence.

##  2 . The US algorithm with strongly stable convergence

First, we define two new relation symbols on the changing direction (CD)
inequalities, “$\stackrel{{\scriptstyle{\rm sgn}(a)}}{{\geqslant}}$” and
“$\stackrel{{\scriptstyle{\rm sgn}(a)}}{{\leqslant}}$”, as follows: For two
functions $h_{1}(x)$ and $h_{2}(x)$ with the same domain ${\mathbb{X}}$,

$h_{1}(x)\stackrel{{\scriptstyle{\rm sgn}(a)}}{{\geqslant}}h_{2}(x)\mbox{
means}\left\\{\\!\\!\begin{array}[]{ll}h_{1}(x)\geqslant h_{2}(x),&\mbox{if }\
a>0,\\\\[5.69054pt] h_{1}(x)=h_{2}(x),&\mbox{if }\ a=0,\\\\[5.69054pt]
h_{1}(x)\leqslant h_{2}(x),&\mbox{if }\ a<0,\end{array}\right.$

and

$h_{1}(x)\stackrel{{\scriptstyle{\rm sgn}(a)}}{{\leqslant}}h_{2}(x)\mbox{
means}\left\\{\\!\\!\begin{array}[]{ll}h_{1}(x)\leqslant h_{2}(x),&\mbox{if }\
a>0,\\\\[5.69054pt] h_{1}(x)=h_{2}(x),&\mbox{if }\ a=0,\\\\[5.69054pt]
h_{1}(x)\geqslant h_{2}(x),&\mbox{if }\ a<0.\end{array}\right.$

For example, $x^{2n-1}\stackrel{{\scriptstyle{\rm sgn}(x)}}{{\geqslant}}0$ for
any positive integer $n$; $\sin(x)\stackrel{{\scriptstyle{\rm
sgn}(x-\pi)}}{{\geqslant}}0$ for any $x\in(0,2\pi)$; and
$\cos(x)\stackrel{{\scriptstyle{\rm sgn}(x-\pi/2)}}{{\leqslant}}0$ for any
$x\in(0,\pi)$ are three illustrative examples.

Next, a sequence $\\{\theta^{(t)}\\}_{t=0}^{\infty}$ is called to weakly
stably converge to a fixed point $\theta^{*}$, if

$|\theta^{(t+1)}-\theta^{*}|<|\theta^{(t)}-\theta^{*}|,\qquad\forall\;t=0,1,\ldots,\infty.$

A sequence $\\{\theta^{(t)}\\}_{t=0}^{\infty}$ is called to strongly stably
converge to a fixed point $\theta^{*}$, if

$\theta^{(0)}<\theta^{(1)}<\cdots<\theta^{(t)}<\cdots\leqslant\theta^{*}\quad\quad\mbox{or}\quad\quad\theta^{*}\leqslant\cdots<\theta^{(t)}<\cdots<\theta^{(1)}<\theta^{(0)}.$

An algorithm is said to possess weakly (or strongly) stable convergence if the
sequence $\\{\theta^{(t)}\\}_{t=0}^{\infty}$ generated by this algorithm
weakly (or strongly) stably converges to a fixed point $\theta^{*}$. It is
clear that a strongly stable convergence implies a weakly stable convergence.

###  2.1 Definition of a $U$-function

Suppose that directly solving the unique root $\theta^{*}$ of the non-linear
equation (1.1) is very difficult. Without loss of generality, we assume that

$\displaystyle g(\theta)>0,\quad\forall\;\theta<\theta^{*}\quad\mbox{and}\quad
g(\theta)<0,\quad\forall\;\theta>\theta^{*}.$ (2.1)

If $g(\theta)<0$ when $\theta<\theta^{*}$, then we can multiplies $-1$ on both
sides of (1.1) and obtain a new equation $-g(\theta)=0$, which satisfies the
assumption (2.1) and has the same root $\theta^{*}$.

The idea of the US algorithm is as follows. Let $\theta^{(t)}$ be the $t$-th
iteration of $\theta^{*}$ and $U(\theta|\theta^{(t)})$ be a real-valued
function of $\theta$ whose form depends on $\theta^{(t)}$. The function
$U(\theta|\theta^{(t)})$ is called a $U$-function for $g(\theta)$ at
$\theta=\theta^{(t)}$ if the following three conditions are satisfied:

$\displaystyle g(\theta)$ $\displaystyle\leqslant$ $\displaystyle
U(\theta|\theta^{(t)}),\qquad\mbox{if }\ \;\theta<\theta^{(t)},$ (2.2)
$\displaystyle g(\theta^{(t)})$ $\displaystyle=$ $\displaystyle
U(\theta^{(t)}|\theta^{(t)})\hskip 18.49411pt\mbox{if }\
\;\theta=\theta^{(t)},$ (2.3) $\displaystyle g(\theta)$
$\displaystyle\geqslant$ $\displaystyle U(\theta|\theta^{(t)}),\qquad\mbox{if
}\ \;\theta>\theta^{(t)}.$ (2.4)

The three conditions can be equivalently rewritten as

$\displaystyle g(\theta)\stackrel{{\scriptstyle{\rm
sgn}(\theta-\theta^{(t)})}}{{\geqslant}}U(\theta|\theta^{(t)}),\quad\forall\;\theta,\theta^{(t)}\in\Theta.$
(2.5)

The conditions (2.2)–(2.3) indicate that the $U(\theta|\theta^{(t)})$ function
majorizes the objective function $g(\theta)$ at $\theta=\theta^{(t)}$ for all
$\theta\leqslant\theta^{(t)}$, and the conditions (2.4)–(2.3) show that the
$U(\theta|\theta^{(t)})$ function minorizes $g(\theta)$ at
$\theta=\theta^{(t)}$ for all $\theta\geqslant\theta^{(t)}$. The geometric
interpretations of the objective function $g(\theta)$ (solid curve) and the
$U$-function $U(\theta|\theta^{(t)})$ (dotted curve) are showed in Figure 1.

The assumption (2.1) implies $0\stackrel{{\scriptstyle{\rm
sgn}(\theta-\theta^{*})}}{{\geqslant}}g(\theta)$ for all
$\theta,\theta^{*}\in\Theta$. By comparing this CD inequality with (2.5), we
conclude that $g(\theta)$ is a $U$-function for 0 at $\theta=\theta^{*}$.
Other properties on $U$-functions are listed in Property 1 below.

Property 1 (Basic properties on $U$-functions).

  * (a)

Let $g(\theta)$ be defined in $\Theta$, then $g(\theta)$ is a $U$-function for
itself at any inner point in $\Theta$;

  * (b)

If $U_{1}(\theta|\theta^{(t)})$ and $U_{2}(\theta|\theta^{(t)})$ are two
$U$-functions for $g(\theta)$ at $\theta=\theta^{(t)}$, then
$[U_{1}(\theta|\theta^{(t)})+U_{2}(\theta|\theta^{(t)})]/2$ is a $U$-function
for $g(\theta)$ at $\theta=\theta^{(t)}$;

  * (c)

For two functions $g_{1}(\theta)$ and $g_{2}(\theta)$ defined in $\Theta$, if
$U_{1}(\theta|\theta^{(t)})$ is a $U$-function for $g_{1}(\theta)$ at
$\theta=\theta^{(t)}$, then $U_{1}(\theta|\theta^{(t)})+g_{2}(\theta)$ is a
$U$-function for $g_{1}(\theta)+g_{2}(\theta)$ at $\theta=\theta^{(t)}$. $\|$

###  2.2 The $U$-equation and the US algorithm

If a $U(\theta|\theta^{(t)})$ function satisfying (2.5) could be found, then
one need only solve the simple U-equation: $U(\theta|\theta^{(t)})=0$ to
obtain its solution $\theta^{(t+1)}$, denoted by

$\displaystyle\theta^{(t+1)}=\mbox{sol}\,\\{U(\theta|\theta^{(t)})=0,\;\forall\;\theta,\theta^{(t)}\in\Theta\\},$
(2.6)

which defines a US algorithm. Therefore, the US algorithm is a general
principle for iteratively seeking the root $\theta^{*}$ with strongly stable
convergence, and its each iteration consists of an upper-crossing step
(U-step) and a solution step (S-step), where

  * U-step:

Find a $U$-function $U(\theta|\theta^{(t)})$ satisfying (2.5);

  * S-step:

Solve the $U$-equation to obtain its root $\theta^{(t+1)}$ as showed by (2.6).

The two-step process is repeated until convergence occurs.

In general, the solution $\theta^{(t+1)}$ in (2.6) has an explicit expression.
Usually, the $U$-equation $U(\theta|\theta^{(t)})=0$ could be any equation
with analytic solution, e.g., a quadratic equation, even a linear equation.
Therefore, the US principle is iteratively solving a sequence of tractable
surrogate equations $U(\theta|\theta^{(t)})=0$ instead of solving the
intractable original non-linear equation $g(\theta)=0$.

Theorem 1 below states that the US algorithm holds two characteristics: (i) It
strongly stably converges to the root $\theta^{*}$; (ii) It does not depend on
any initial values in $\Theta$, in contrast to Newton’s method. Its proof is
put in Appendix A.

Theorem 1 (Strongly stable convergence). Let $\theta^{*}$ be the unique root
of the equation (1.1), and $U(\theta|\theta^{(t)})$ be one $U$-function for
$g(\theta)$ at $\theta=\theta^{(t)}$ in $\Theta$. Then, for any initial value
$\theta^{(0)}$ in $\Theta$, the sequence $\\{\theta^{(t+1)}\\}_{t=0}^{\infty}$
specified by (2.6) strongly stably converges to $\theta^{*}$. $\|$

##  3 . Three methods for constructing $U$-functions

For the same objective function $g(\theta)$, it is possible that there exist
several $U$-functions. In other words, for the equation (1.1), it may exist
more than one US algorithm. We know that the key for a US algorithm is how to
find a $U$-function. In this section, we propose three general methods to
construct $U$-functions based on the objective function $g(\theta)$.

###  3.1 The FLB function method

Suppose the first–order derivative $g^{\prime}(\cdot)$ exists and it is
bounded by some function $b(\cdot)$; i.e.,

$\displaystyle g^{\prime}(\theta)\geqslant
b(\theta),\quad\forall\;\theta\in\Theta,$ (3.1)

where $b(\theta)$ is called the first-derivative lower bound (FLB) function.
It is easy to verify that

$\displaystyle U(\theta|\theta^{(t)})\triangleq
g(\theta^{(t)})+\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z,\quad\forall\;\theta,\theta^{(t)}\in\Theta,$
(3.2)

is a $U$-function for $g(\theta)$ at $\theta=\theta^{(t)}$. In fact, we have

$\displaystyle g(\theta)-U(\theta|\theta^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle[g(\theta)-g(\theta^{(t)})]-\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z\;=\;\int_{\theta^{(t)}}^{\theta}g^{\prime}(z)\,\mbox{d}z-\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z$
$\displaystyle=$
$\displaystyle\int_{\theta^{(t)}}^{\theta}[g^{\prime}(z)-b(z)]\,\mbox{d}z\;\left\\{\begin{array}[]{l}\leqslant
0,\quad\mbox{if }\ \theta<\theta^{(t)},\mbox{ indicating
}(\ref{eqnfive2.2}),\\\\[11.38109pt] =0,\quad\mbox{if }\
\theta=\theta^{(t)},\mbox{ indicating }(\ref{eqnfive2.3}),\\\\[11.38109pt]
\geqslant 0,\quad\mbox{if }\ \theta>\theta^{(t)},\mbox{ indicating
}(\ref{eqnfive2.4})\end{array}\right.$
$\displaystyle\stackrel{{\scriptstyle{\rm
sgn}(\theta-\theta^{(t)})}}{{\geqslant}}$ $\displaystyle
0,\quad\forall\;\theta,\theta^{(t)}\in\Theta,$

where we have used the result
$g^{\prime}(z)-b(z)\stackrel{{\scriptstyle(\ref{eqnfive3.1})}}{{\geqslant}}0$
for all $z\in\Theta$ in the above CD inequality.

####  3.1.1 Two special cases of the FLB function $b(\theta)$

Let $b(\theta)=b_{1}+b_{2}\theta$ be a linear function. Then the $U$-function
defined by (3.2) becomes a quadratic function

$\displaystyle
U(\theta|\theta^{(t)})=\frac{1}{2}b_{2}\theta^{2}+b_{1}\theta+c^{(t)}_{1}\quad\mbox{with}\quad
c^{(t)}_{1}\triangleq
g(\theta^{(t)})-\frac{1}{2}b_{2}\theta^{(t)2}-b_{1}\theta^{(t)},$ (3.4)

if $b_{2}\neq 0$. From (2.6), the $U$-equation $U(\theta|\theta^{(t)})=0$ has
two solutions

$\displaystyle\theta^{(t+1)}=\frac{-b_{1}\pm\sqrt{b_{1}^{2}-2b_{2}c^{(t)}_{1}}}{b_{2}},$
(3.5)

but only one is the desired US iteration.

In particular, if $b_{2}=0$ (i.e., $g^{\prime}(\theta)\geqslant b_{1}$ for all
$\theta\in\Theta$, where $b_{1}$ is called the FLB constant), then (3.4)
becomes to a linear function

$\displaystyle
U(\theta|\theta^{(t)})=g(\theta^{(t)})+b_{1}(\theta-\theta^{(t)})$ (3.6)

and the corresponding US iteration is given by

$\displaystyle\theta^{(t+1)}=\theta^{(t)}-\frac{g(\theta^{(t)})}{b_{1}},\quad
b_{1}\neq 0,$ (3.7)

which can be viewed as the NR iteration (1.2) by replacing
$g^{\prime}(\theta^{(t)})$ with $b_{1}$.

####  3.1.2 Motivation for constructing the $U$-function (3.2)

Under the assumption (3.1), the key for constructing the US algorithm is how
to obtain $U(\theta|\theta^{(t)})$ function specified by (3.2). Our motivation
is from the mean value theorem (i.e., the first–order Taylor expansion). It
states that there exists a point $z^{*}$ between $\theta$ and $\theta^{(t)}$
such that

$\displaystyle
g(\theta)-g(\theta^{(t)})=(\theta-\theta^{(t)})g^{\prime}(z^{*})=\int_{\theta^{(t)}}^{\theta}g^{\prime}(z)\,\mbox{d}z.$
(3.8)

Rewriting (3.8), we obtain

$\displaystyle g(\theta)$ $\displaystyle=$ $\displaystyle
g(\theta^{(t)})+\int_{\theta^{(t)}}^{\theta}g^{\prime}(z)\,\mbox{d}z$
$\displaystyle\stackrel{{\scriptstyle{\rm
sgn}(\theta-\theta^{(t)})}}{{\geqslant}}$ $\displaystyle
g(\theta^{(t)})+\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z\;\triangleq\;U(\theta|\theta^{(t)}),$

where the result
$g^{\prime}(z)-b(z)\stackrel{{\scriptstyle(\ref{eqnfive3.1})}}{{\geqslant}}0$
is used in the above CD inequality.

###  3.2 The SLUB constants method

Suppose that $g^{\prime\prime}(\cdot)$ exists and there are two constants
$\\{b_{21},b_{22}\\}$ (which are called the second-derivative lower–upper
bound (SLUB) constants) such that

$\displaystyle b_{21}\leqslant g^{\prime\prime}(\theta)\leqslant
b_{22},\quad\forall\;\theta\in\Theta.$ (3.9)

Let $I(\cdot)$ denote the indicator function and define a step constant

$\displaystyle b_{2}(\theta|\theta^{(t)})\triangleq b_{22}\cdot
I(\theta\leqslant\theta^{(t)})+b_{21}\cdot I(\theta>\theta^{(t)}).$ (3.10)

It is easy to verify that the quadratic function

$\displaystyle U(\theta|\theta^{(t)})\triangleq
g(\theta^{(t)})+g^{\prime}(\theta^{(t)})(\theta-\theta^{(t)})+\frac{1}{2}b_{2}(\theta|\theta^{(t)})(\theta-\theta^{(t)})^{2},\quad\forall\;\theta,\theta^{(t)}\in\Theta,$
(3.11)

is a $U$-function for $g(\theta)$ at $\theta=\theta^{(t)}$.

Note that the second–order Taylor’s expansion

$g(\theta)=g(\theta^{(t)})+g^{\prime}(\theta^{(t)})(\theta-\theta^{(t)})+\frac{1}{2}g^{\prime\prime}(\tilde{\theta})(\theta-\theta^{(t)})^{2},\quad\forall\;\theta,\theta^{(t)}\in\Theta,$

where $\tilde{\theta}$ lies between $\theta$ and $\theta^{(t)}$, we only need
to prove that $b_{2}(\theta|\theta^{(t)})(\theta-\theta^{(t)})^{2}$ is a
$U$-function for $g^{\prime\prime}(\tilde{\theta})(\theta-\theta^{(t)})^{2}$
at $\theta=\theta^{(t)}$. In fact, by using (3.9) and (3.10), we have

$\displaystyle
g^{\prime\prime}(\tilde{\theta})(\theta-\theta^{(t)})^{2}\left\\{\begin{array}[]{ll}\leqslant
b_{22}(\theta-\theta^{(t)})^{2},&\mbox{if }\ \theta<\theta^{(t)},\mbox{
indicating }(\ref{eqnfive2.2}),\\\\[5.69054pt] =0,&\mbox{if }\
\theta=\theta^{(t)},\mbox{ indicating }(\ref{eqnfive2.3}),\\\\[5.69054pt]
\geqslant b_{21}(\theta-\theta^{(t)})^{2},&\mbox{if }\
\theta>\theta^{(t)},\mbox{ indicating }(\ref{eqnfive2.4})\end{array}\right.$
$\displaystyle\stackrel{{\scriptstyle{\rm
sgn}(\theta-\theta^{(t)})}}{{\geqslant}}$ $\displaystyle
b_{2}(\theta|\theta^{(t)})(\theta-\theta^{(t)})^{2},\quad\forall\;\theta,\theta^{(t)}\in\Theta.$

###  3.3 The TLB constant method

Suppose that $g^{\prime\prime\prime}(\theta)$ exists and is bounded by some
constant $b_{3}$; i.e.,

$\displaystyle g^{\prime\prime\prime}(\theta)\geqslant
b_{3},\quad\forall\;\theta\in\Theta,$ (3.13)

where $b_{3}$ is called the third-derivative lower bound (TLB) constant. It is
easy to verify that the cubic function (if $b_{3}\neq 0$)

$\displaystyle U(\theta|\theta^{(t)})$ $\displaystyle\triangleq$
$\displaystyle
g(\theta^{(t)})+g^{\prime}(\theta^{(t)})(\theta-\theta^{(t)})+\frac{1}{2}g^{\prime\prime}(\theta^{(t)})(\theta-\theta^{(t)})^{2}$
(3.14)
$\displaystyle+\;\frac{1}{6}b_{3}(\theta-\theta^{(t)})^{3},\quad\forall\;\theta,\theta^{(t)}\in\Theta,$

is a $U$-function for $g(\theta)$ at $\theta=\theta^{(t)}$. If $b_{3}=0$, then
(3.14) becomes a quadratic function.

In fact, by applying the third–order Taylor expansion

$\displaystyle
g(\theta)=g(\theta^{(t)})+g^{\prime}(\theta^{(t)})(\theta-\theta^{(t)})+\frac{1}{2}g^{\prime\prime}(\theta^{(t)})(\theta-\theta^{(t)})^{2}+\frac{1}{6}g^{\prime\prime\prime}(\tilde{\theta})(\theta-\theta^{(t)})^{3},$
(3.15)

where $\tilde{\theta}$ lies between $\theta$ and $\theta^{(t)}$, we obtain

$\displaystyle g(\theta)-U(\theta|\theta^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.12})\ \&\
(\ref{eqnfive3.13})}}{{=}}$
$\displaystyle\frac{(\theta-\theta^{(t)})^{3}}{6}[g^{\prime\prime\prime}(\tilde{\theta})-b_{3}]\left\\{\begin{array}[]{ll}\leqslant
0,&\mbox{if }\ \theta<\theta^{(t)},\mbox{ indicating
}(\ref{eqnfive2.2}),\\\\[5.69054pt] =0,&\mbox{if }\ \theta=\theta^{(t)},\mbox{
indicating }(\ref{eqnfive2.3}),\\\\[5.69054pt] \geqslant 0,&\mbox{if }\
\theta>\theta^{(t)},\mbox{ indicating }(\ref{eqnfive2.4})\end{array}\right.$
$\displaystyle\stackrel{{\scriptstyle{\rm
sgn}(\theta-\theta^{(t)})}}{{\geqslant}}$ $\displaystyle
0,\quad\forall\;\theta,\theta^{(t)}\in\Theta,$

where we have used the result
$g^{\prime\prime\prime}(\theta)-b_{3}\stackrel{{\scriptstyle(\ref{eqnfive3.11})}}{{\geqslant}}0$
for all $\theta\in\Theta$ in the above CD inequality.

##  4 . Applications

###  4.1 Calculation of quantile in continuous distributions

Let $f(x|{\boldsymbol{\theta}})$ and $F(x|{\boldsymbol{\theta}})$ respectively
denote the pdf and cdf of a population distribution. Given
${\boldsymbol{\theta}}$ and a real number $p\in(0,1)$, the $p$-th quantile
$\xi_{p}$ of the distribution can be calculated as

$\xi_{p}=\mbox{sol}\,\Big{\\{}g(x|{\boldsymbol{\theta}})\triangleq
p-F(x|{\boldsymbol{\theta}})=0,\;\forall\;x\in{\mathbb{X}}\Big{\\}},$

where we need to solve an integral equation. In this subsection, we will
demonstrate that the US algorithm can be employed to iteratively find the
$p$-th quantile $\xi_{p}$.

Assume that
$g^{\prime}(x|{\boldsymbol{\theta}})\triangleq\,\mbox{d}g(x|{\boldsymbol{\theta}})/\,\mbox{d}x$
is bounded by some function $b(x|{\boldsymbol{\theta}})$; i.e.,

$\displaystyle
g^{\prime}(x|{\boldsymbol{\theta}})=-f(x|{\boldsymbol{\theta}})\geqslant
b(x|{\boldsymbol{\theta}}),\quad\forall\;x\in{\mathbb{X}}.$ (4.1)

Then, from (2.6), the US iteration is to solve

$\displaystyle
x^{(t+1)}=\mbox{sol}\left\\{U(x|x^{(t)},{\boldsymbol{\theta}})\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}g(x^{(t)}|{\boldsymbol{\theta}})+\int_{x^{(t)}}^{x}b(z|{\boldsymbol{\theta}})\,\mbox{d}z=0,\;\forall\;x,x^{(t)}\in{\mathbb{X}}\right\\}.$
(4.2)

In this subsection, we will consider two cases: (i) The mode $x_{\rm mod}$ of
the density $f(x|{\boldsymbol{\theta}})$ exists; (ii) The mode $x_{\rm mod}$
of $f(x|{\boldsymbol{\theta}})$ does not exist.

####  4.1.1 The mode of the density exists

Let $x_{\rm mod}$ denote the mode of the density $f(x|{\boldsymbol{\theta}})$,
then we have $-f(x|{\boldsymbol{\theta}})\geqslant-f(x_{\rm
mod}|{\boldsymbol{\theta}})$ for all $x\in{\mathbb{X}}$. From (4.2), the US
iteration is to solve

$\displaystyle x^{(t+1)}$ $\displaystyle=$
$\displaystyle\mbox{sol}\Big{\\{}g(x^{(t)}|{\boldsymbol{\theta}})-f(x_{\rm
mod}|{\boldsymbol{\theta}})(x-x^{(t)})=0,\;\forall\;x,x^{(t)}\in{\mathbb{X}}\Big{\\}}$
(4.3) $\displaystyle=$ $\displaystyle
x^{(t)}-\frac{F(x^{(t)}|{\boldsymbol{\theta}})-p}{f(x_{\rm
mod}|{\boldsymbol{\theta}})}.$

Example 2 (Quantile of the normal distribution). To find the $p$-th quantile
of the normal distribution with mean $\mu$ and variance $\sigma^{2}$, we
denote its cdf and pdf by $\Phi(x|\mu,\sigma^{2})$ and
$\phi(x|\mu,\sigma^{2})=[1/(\sqrt{2\pi}\sigma)]\exp[-(x-\mu)^{2}/(2\sigma^{2})]$,
respectively. Note that the mode of $\phi(x|\mu,\sigma^{2})$ is given by
$x_{\rm mod}=\mu$, then (4.3) becomes

$\displaystyle
x^{(t+1)}=x^{(t)}-\sqrt{2\pi}\sigma\Big{[}\Phi(x^{(t)}|\mu,\sigma^{2})-p\Big{]}.$
(4.4)

Example 3 (Quantile of the skew normal distribution). To find the $p$-th
quantile of the skew normal distribution (Azzalini 1985) with pdf

$\displaystyle f(x|\mu,\sigma^{2},\alpha)$ $\displaystyle=$
$\displaystyle\frac{2}{\sigma\sqrt{2\pi}}\exp{\left[-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right]}\Phi\left(\alpha\frac{x-\mu}{\sigma}\right),$
(4.5)

where $x,\mu\in{\mathbb{R}}$, $\sigma^{2}>0$, $\alpha\geqslant 0$ and
$\Phi(\cdot)$ is the cdf of $N(0,1)$, in Appendix B, we proposed an MM
algorithm to iteratively calculate the mode $x_{\rm mod}$ of the skew normal
density as shown in (B.1). From (4.3), the US iteration for calculating the
$p$-th quantile $\xi_{p}$ is

$x^{(t+1)}=x^{(t)}-\frac{F(x^{(t)}|\mu,\sigma^{2},\alpha)-p}{f(x_{\rm
mod}|\mu,\sigma^{2},\alpha)},$

where $F(x|\mu,\sigma^{2},\alpha)$ is the cdf of the skew normal distribution.

####  4.1.2 The mode of the density does not exist

Let $F(x|\alpha,\beta)$ denote the cdf of the beta distribution
$\mbox{Beta}(\alpha,\beta)$ and its pdf is

$\displaystyle f(x|\alpha,\beta)$ $\displaystyle=$
$\displaystyle\frac{1}{B(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1},\quad
x\in(0,1),\;\alpha>0,\;\beta>0.$ (4.6)

Given $\alpha\in(0,1),\beta\in(0,1)$ and $p\in(0,1)$, suppose that we want to
calculate the $p$-th quantile of the distribution as
$\xi_{p}=\mbox{sol}\,\big{\\{}g(x|\alpha,\beta)\triangleq
p-F(x|\alpha,\beta)=0,\;\forall\;x\in(0,1)\big{\\}}$. We know that if at least
one of $\alpha$ and $\beta$ is less than 1, the mode of the beta density
$f(x|\alpha,\beta)$ does not exist. Thus, the US iteration specified by (4.3)
can not be used for such cases.

We employ (4.1)–(4.2) to calculate $\xi_{p}$. First, it is easy to verify the
following inequality:

$\displaystyle\frac{1}{x(1-x)}\leqslant\frac{1}{2x^{2}}+\frac{1}{2(1-x)^{2}},\quad\forall\;x\in(0,1),$
(4.7)

where the equality holds iff $x=1/2$. Next,

$\displaystyle g^{\prime}(x|\alpha,\beta)$ $\displaystyle=$
$\displaystyle-f(x|\alpha,\beta)=-f(x|\alpha+1,\beta+1)\frac{B(\alpha+1,\beta+1)}{B(\alpha,\beta)}\cdot\frac{1}{x(1-x)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive4.7})}}{{\geqslant}}$
$\displaystyle-f(x|\alpha+1,\beta+1)\frac{B(\alpha+1,\beta+1)}{B(\alpha,\beta)}\cdot\left[\frac{1}{2x^{2}}+\frac{1}{2(1-x)^{2}}\right]$
$\displaystyle\geqslant$ $\displaystyle-f(x_{\rm
mod,1}|\alpha+1,\beta+1)\frac{B(\alpha+1,\beta+1)}{2B(\alpha,\beta)}\left[\frac{1}{x^{2}}+\frac{1}{(1-x)^{2}}\right]$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive4.8})}}{{=}}$
$\displaystyle-\frac{a}{x^{2}}-\frac{a}{(1-x)^{2}}\;\triangleq\;b(x|\alpha,\beta),\quad\forall\;x\in(0,1),$

where $x_{\rm mod,1}\triangleq\alpha/(\alpha+\beta)$ denotes the mode of the
beta density $f(x|\alpha+1,\beta+1)$ and

$\displaystyle a\triangleq f(x_{\rm
mod,1}|\alpha+1,\beta+1)\frac{B(\alpha+1,\beta+1)}{2B(\alpha,\beta)}.$ (4.8)

Finally, from (4.2), the US iteration for calculating the $p$-th quantile
$\xi_{p}$ is

$\displaystyle x^{(t+1)}$ $\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(x^{(t)}|\alpha,\beta)+\int_{x^{(t)}}^{x}b(z|\alpha,\beta)\,\mbox{d}z=0,\;\forall\;x,x^{(t)}\in(0,1)\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(x^{(t)}|\alpha,\beta)+\frac{a}{x}-\frac{a}{1-x}-\frac{a}{x^{(t)}}+\frac{a}{1-x^{(t)}}=0,\;\forall\;x,x^{(t)}\in(0,1)\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{a_{1}^{(t)}x^{2}+a_{2}^{(t)}x+a=0,\;\forall\;x,x^{(t)}\in(0,1)\right\\},$

where $a_{1}^{(t)}\triangleq
F(x^{(t)}|\alpha,\beta)-p+a/x^{(t)}-a/(1-x^{(t)})$ and $a_{2}^{(t)}\triangleq-
a_{1}^{(t)}-2a$.

###  4.2 Calculation of exact $p$-values for skew null distributions

####  4.2.1 The chi-squared null distribution

Let $\\{X_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle{\rm
iid}}}{{\sim}}N(\mu,\sigma^{2})$ with unknown $\\{\mu,\sigma^{2}\\}$ and we
want to test $H_{0}$: $\sigma^{2}=\sigma_{0}^{2}$ against $H_{1}$:
$\sigma^{2}\neq\sigma_{0}^{2}$. The test statistic is $\chi=\nu
S^{2}/\sigma_{0}^{2}$ with corresponding observed value $\chi_{\rm
obs}\triangleq\nu s^{2}/\sigma_{0}^{2}$, where
$S^{2}=(1/\nu)\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}$ is the sample variance,
$s^{2}$ is its observed value and $\nu\;\hat{=}\;n-1$. Under $H_{0}$, we have
$\chi\sim\chi^{2}(\nu)$, which is called chi-squared null distribution. We
denote the pdf of $\chi^{2}(\nu)$ by
$h_{\nu}(x)=[2^{\nu/2}\Gamma(\nu/2)]^{-1}x^{\nu/2-1}\mbox{e}^{-x/2}$ for
$x>0$, whose mode is $\nu-2\triangleq m_{0}$.

(a) Case I: $\chi_{\rm obs}<m_{0}$

The exact $p$-value can be calculated as
$p\mbox{-value}=\Pr\\{\chi^{2}(\nu)\leqslant\chi_{\rm
obs}\\}+\Pr\\{\chi^{2}(\nu)\geqslant x_{u}\\}$, where $x_{u}\,(>m_{0})$
satisfies $h_{\nu}(x_{u})=h_{\nu}(\chi_{\rm obs})$; i.e.,
$m_{0}\log(x_{u})-x_{u}=m_{0}\log(\chi_{\rm obs})-\chi_{\rm
obs}\;\triangleq\;c_{0}$. Thus, finding $x_{u}$ is equivalent to finding the
root to the equation $g(x_{u})=c_{0}-m_{0}\log(x_{u})+x_{u}=0$. Since
$g^{\prime}(x_{u})=-m_{0}/x_{u}+1>-m_{0}/x_{u}\;\triangleq\;b(x_{u})$, the US
iteration for calculating $x_{u}$ is

$\displaystyle x_{u}^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(x^{(t)}_{u})+\textstyle\int_{x^{(t)}_{u}}^{x_{u}}b(z)\,\mbox{d}z=0,\;\forall\;x_{u},x^{(t)}_{u}>m_{0}\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\Big{\\{}g(x^{(t)}_{u})-m_{0}\log(x_{u})+m_{0}\log(x^{(t)}_{u})=0,\;\forall\;x_{u},x^{(t)}_{u}>m_{0}\Big{\\}}$
$\displaystyle=$ $\displaystyle
x^{(t)}_{u}\exp\left[g(x^{(t)}_{u})/m_{0}\right].$

(b) Case II: $\chi_{\rm obs}>m_{0}$

The exact $p$-value can be calculated as
$p\mbox{-value}=\Pr\\{\chi^{2}(\nu)\leqslant
x_{l}\\}+\Pr\\{\chi^{2}(\nu)\geqslant\chi_{\rm obs}\\}$, where
$x_{l}\in(0,m_{0})$ satisfies $h_{\nu}(x_{l})=h_{\nu}(\chi_{\rm obs})$; i.e.,
$m_{0}\log(x_{l})-x_{l}=m_{0}\log(\chi_{\rm obs})-\chi_{\rm
obs}\;\triangleq\;c_{0}$. Thus, finding $x_{l}$ is equivalent to finding the
root to the equation $g(x_{l})=m_{0}\log(x_{l})-x_{l}-c_{0}=0$. Since
$g^{\prime}(x_{l})=m_{0}/x_{l}-1>-1\;\triangleq\;b(x_{l})$, the US iteration
for calculating $x_{l}$ is

$\displaystyle x_{l}^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(x^{(t)}_{l})+\textstyle\int_{x^{(t)}_{l}}^{x_{l}}b(z)\,\mbox{d}z=0,\;\forall\;x_{l},x^{(t)}_{l}\in(0,m_{0})\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\Big{\\{}g(x^{(t)}_{l})-x_{l}+x^{(t)}_{l}=0,\;\forall\;x_{l},x^{(t)}_{l}\in(0,m_{0})\Big{\\}}$
$\displaystyle=$ $\displaystyle
g(x^{(t)}_{l})+x^{(t)}_{l}=m_{0}\log\big{(}x^{(t)}_{l}/\chi_{\rm
obs}\big{)}+\chi_{\rm obs}.$

####  4.2.2 The F null distribution

Let $\\{X_{ij}\\}_{j=1}^{n_{i}}\stackrel{{\scriptstyle{\rm
iid}}}{{\sim}}N(\mu_{i},\sigma_{i}^{2})$ for $i=1,2$, and the two random
samples be independent, where $\\{\mu_{i},\sigma_{i}^{2}\\}$ are unknown. We
want to test $H_{0}$: $\sigma_{1}^{2}=\sigma_{2}^{2}$ against $H_{1}$:
$\sigma_{1}^{2}\neq\sigma_{2}^{2}$, and the test statistic is
$F=S_{1}^{2}/S_{2}^{2}$ with corresponding observed value $F_{\rm
obs}\triangleq s_{1}^{2}/s_{2}^{2}$, where $S_{i}^{2}$ is the $i$-th sample
variance, $s_{i}^{2}$ is its observed value and $\nu_{i}\triangleq n_{i}-1$.
Under $H_{0}$, we have $F\sim F(\nu_{1},\nu_{2})$, which is called $F$ null
distribution. We denote the pdf of $F(\nu_{1},\nu_{2})$ by

$h_{\nu_{1},\nu_{2}}(x)=\frac{\Gamma(\nu_{12}/2)\nu^{\nu_{1}/2}}{\Gamma(\nu_{1}/2)\Gamma(\nu_{2}/2)}x^{\nu_{1}/2-1}(1+\nu
x)^{-\nu_{12}/2},\quad
x>0,\;\;\nu\triangleq\frac{\nu_{1}}{\nu_{2}},\;\;\nu_{12}\triangleq\nu_{1}+\nu_{2}.$

The mode of $h_{\nu_{1},\nu_{2}}(x)$ is
$m_{1}\nu_{2}/[\nu_{1}(\nu_{2}+2)]\triangleq m_{0}$, where
$m_{1}\triangleq\nu_{1}-2$.

(a) Case I: $F_{\rm obs}<m_{0}$

The exact $p$-value can be calculated as
$p\mbox{-value}=\Pr\\{F(\nu_{1},\nu_{2})\leqslant F_{\rm
obs}\\}+\Pr\\{F(\nu_{1},\nu_{2})\geqslant x_{u}\\}$, where $x_{u}\,(>m_{0})$
satisfies $h_{\nu_{1},\nu_{2}}(x_{u})=h_{\nu_{1},\nu_{2}}(F_{\rm obs})$; i.e.,
$m_{1}\log(x_{u})-\nu_{12}\log(1+\nu x_{u})=m_{1}\log(F_{\rm
obs})-\nu_{12}\log(1+\nu F_{\rm obs})\triangleq c_{0}$. Thus, finding $x_{u}$
is equivalent to finding the root to the equation
$g(x_{u})=c_{0}-m_{1}\log(x_{u})+\nu_{12}\log(1+\nu x_{u})=0$. Since
$g^{\prime}(x_{u})=-m_{1}/x_{u}+\nu_{12}\nu/(1+\nu
x_{u})>-m_{1}/x_{u}\;\triangleq\;b(x_{u})$, the US iteration for calculating
$x_{u}$ is

$\displaystyle x_{u}^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(x^{(t)}_{u})+\textstyle\int_{x^{(t)}_{u}}^{x_{u}}b(z)\,\mbox{d}z=0,\;\forall\;x_{u},x^{(t)}_{u}>m_{0}\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\Big{\\{}g(x^{(t)}_{u})-m_{1}\log(x_{u})+m_{1}\log(x^{(t)}_{u})=0,\;\forall\;x_{u},x^{(t)}_{u}>m_{0}\Big{\\}}$
$\displaystyle=$ $\displaystyle x^{(t)}_{u}\exp[g(x^{(t)}_{u})/m_{1}].$

(b) Case II: $F_{\rm obs}>m_{0}$

The exact $p$-value can be calculated as
$p\mbox{-value}=\Pr\\{F(\nu_{1},\nu_{2})\leqslant
x_{l}\\}+\Pr\\{F(\nu_{1},\nu_{2})\geqslant F_{\rm obs}\\}$, where
$x_{l}\in(0,m_{0})$ satisfies
$h_{\nu_{1},\nu_{2}}(x_{l})=h_{\nu_{1},\nu_{2}}(F_{\rm obs})$; i.e.,
$m_{1}\log(x_{l})-\nu_{12}\log(1+\nu x_{l})=m_{1}\log(F_{\rm
obs})-\nu_{12}\log(1+\nu F_{\rm obs})\triangleq c_{0}$. Thus, finding $x_{l}$
is equivalent to finding the root to the equation
$g(x_{l})=m_{1}\log(x_{l})-\nu_{12}\log(1+\nu x_{l})-c_{0}=0$. Since
$g^{\prime}(x_{l})=m_{1}/x_{l}-\nu_{12}\nu/(1+\nu
x_{l})>m_{1}/x_{l}-\nu_{12}/x_{l}=-(2+\nu_{2})/x_{l}\;\triangleq\;b(x_{l})$,
the US iteration for calculating $x_{l}$ is

$\displaystyle x_{l}^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(x^{(t)}_{l})+\textstyle\int_{x^{(t)}_{l}}^{x_{l}}b(z)\,\mbox{d}z=0,\;\forall\;x_{l},x^{(t)}_{l}\in(0,m_{0})\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\Big{\\{}g(x^{(t)}_{l})-(2+\nu_{2})\log(x_{l}/x^{(t)}_{l})=0,\;\forall\;x_{l},x^{(t)}_{l}\in(0,m_{0})\Big{\\}}$
$\displaystyle=$ $\displaystyle x^{(t)}_{l}\exp[g(x^{(t)}_{l})/(2+\nu_{2})].$

###  4.3 MLEs of parameters in a class of continuous distributions

####  4.3.1 MLE of $\alpha$ in gamma distribution

Let $\\{X_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle{\rm
iid}}}{{\sim}}\mbox{Gamma}(\alpha,\beta)$ with pdf
$\beta^{\alpha}x^{\alpha-1}\mbox{e}^{-\beta x}/\Gamma(\alpha)$, $x>0$, where
$\\{\alpha,\beta\\}$ are two positive parameters. The log-likelihood function
$\ell(\alpha,\beta)=nG({\boldsymbol{x}})(\alpha-1)-n\bar{x}\beta+n\left[\alpha\log\beta-\log\Gamma(\alpha)\right]$,
where
$G({\boldsymbol{x}})=G(x_{1},\ldots,x_{n})\triangleq(1/n)\sum_{i=1}^{n}\log(x_{i})$
and $\bar{x}=(1/n)\sum_{i=1}^{n}x_{i}$. Let

$\displaystyle 0=\frac{\partial\ell(\alpha,\beta)}{\partial\alpha}$
$\displaystyle=$ $\displaystyle
nG({\boldsymbol{x}})+n\log(\beta)-n\psi(\alpha),$ (4.9) $\displaystyle
0=\frac{\partial\ell(\alpha,\beta)}{\partial\beta}$ $\displaystyle=$
$\displaystyle-n\bar{x}+\frac{n\alpha}{\beta}.$ (4.10)

Given $\alpha$, from (4.10), we have $\beta=\alpha/\bar{x}$. Given $\beta$,
from (4.9), the MLE of $\alpha$ is the root of the equation:
$g(\alpha)=c_{0}-\psi(\alpha)=0$ for $\alpha>0$, where $c_{0}\triangleq
G({\boldsymbol{x}})+\log(\beta)$ is a known constant, $\psi(\cdot)$ is the
digamma function defined as

$\displaystyle\psi(\alpha)\triangleq\frac{\Gamma^{\prime}(\alpha)}{\Gamma(\alpha)}=-\gamma+\sum_{m=0}^{\infty}\left(\frac{1}{m+1}-\frac{1}{m+\alpha}\right)$
(4.11)

here $\gamma=\lim_{n\to\infty}\left(\sum_{m=1}^{n}m^{-1}-\log n\right)\approx
0.5772$ is the Euler–Mascheroni constant. Since

$\displaystyle g^{\prime}(\alpha)$ $\displaystyle=$
$\displaystyle-\psi^{\prime}(\alpha)=-\sum_{m=0}^{\infty}\frac{1}{(m+\alpha)^{2}}=-\frac{1}{\alpha^{2}}-\sum_{m=1}^{\infty}\frac{1}{(m+\alpha)^{2}}$
$\displaystyle>$
$\displaystyle-\frac{1}{\alpha^{2}}-\sum_{m=1}^{\infty}\frac{1}{m^{2}}=-\frac{1}{\alpha^{2}}-\frac{\pi^{2}}{6}\;\triangleq\;b(\alpha),$

the US iteration for calculating the MLE $\hat{\alpha}$ is

$\displaystyle\alpha^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(\alpha^{(t)})+\int_{\alpha^{(t)}}^{\alpha}b(z)\,\mbox{d}z=0,\;\forall\;\alpha,\alpha^{(t)}>0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(\alpha^{(t)})+\frac{1}{\alpha}-\frac{\pi^{2}}{6}\alpha-\frac{1}{\alpha^{(t)}}+\frac{\pi^{2}}{6}\alpha^{(t)}=0,\;\forall\;\alpha,\alpha^{(t)}>0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{(\pi^{2}/6)\alpha^{2}-a_{3}^{(t)}\alpha-1=0,\;\forall\;\alpha,\alpha^{(t)}>0\right\\},$

where
$a_{3}^{(t)}=c_{0}-\psi(\alpha^{(t)})+(\pi^{2}/6)\alpha^{(t)}-1/\alpha^{(t)}$.

####  4.3.2 MLE of $\theta$ in Weibull distribution

Let $\\{X_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle{\rm
iid}}}{{\sim}}\mbox{Weibull}(\theta,\lambda)$ with pdf
$(\theta/\lambda)(x/\lambda)^{\theta-1}\exp[-(x/\lambda)^{\theta}]$, $x>0$,
where $\\{\theta,\lambda\\}$ are two positive parameters (Hallinan 1993, Rinne
2008). The log-likelihood function is
$\ell(\theta,\lambda)=n\log(\theta)-n\theta\log(\lambda)+nG({\boldsymbol{x}})(\theta-1)-\sum_{i=1}^{n}(x_{i}/\lambda)^{\theta}$,
where $G({\boldsymbol{x}})\triangleq(1/n)\sum_{i=1}^{n}\log(x_{i})$. Let
$0=\partial\ell(\theta,\lambda)/\partial\lambda=-n\theta\lambda^{-1}+\theta\lambda^{-\theta-1}\sum_{i=1}^{n}x_{i}^{\theta}$.
Given $\theta$, we obtain
$\lambda=[(1/n)\sum_{i=1}^{n}x_{i}^{\theta}]^{1/\theta}$. Replacing $\lambda$
in $\ell(\theta,\lambda)$ by $[(1/n)\sum_{i=1}^{n}x_{i}^{\theta}]^{1/\theta}$,
we know that $\ell(\theta,\lambda)$ will reduce to $\ell(\theta)=\theta
nG({\boldsymbol{x}})+n\log(\theta)-n\log(\textstyle\sum_{i=1}^{n}x_{i}^{\theta})+\mbox{
constant}$. Thus, the MLE of $\theta$ is the root of $\ell^{\prime}(\theta)=0$
or the root of the equation

$g(\theta)=G({\boldsymbol{x}})+\frac{1}{\theta}-\frac{\sum_{i=1}^{n}x_{i}^{\theta}\log(x_{i})}{\sum_{i=1}^{n}x_{i}^{\theta}}=0,\quad\theta>0.$

Since

$\displaystyle g^{\prime}(\theta)$ $\displaystyle=$
$\displaystyle-\frac{1}{\theta^{2}}-\frac{[\sum_{i=1}^{n}x_{i}^{\theta}(\log
x_{i})^{2}](\sum_{i=1}^{n}x_{i}^{\theta})-(\sum_{i=1}^{n}x_{i}^{\theta}\log
x_{i})^{2}}{(\sum_{i=1}^{n}x_{i}^{\theta})^{2}}$ $\displaystyle\geqslant$
$\displaystyle-\frac{1}{\theta^{2}}-\frac{\sum_{i=1}^{n}x_{i}^{\theta}(\log
x_{i})^{2}}{\sum_{i=1}^{n}x_{i}^{\theta}}$ $\displaystyle\geqslant$
$\displaystyle-\frac{1}{\theta^{2}}-\max_{1\leqslant i\leqslant n}(\log
x_{i})^{2}\;\triangleq\;-\frac{1}{\theta^{2}}-T_{\max}\;\triangleq\;b(\theta),$

the US iteration for calculating the MLE $\hat{\theta}$ is

$\displaystyle\theta^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})+\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z=0,\;\forall\;\theta,\theta^{(t)}>0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})+\frac{1}{\theta}-T_{\max}\theta-\frac{1}{\theta^{(t)}}+T_{\max}\theta^{(t)}=0,\;\forall\;\theta,\theta^{(t)}>0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{T_{\max}\theta^{2}-a_{4}^{(t)}\theta-1=0,\;\forall\;\theta,\theta^{(t)}>0\right\\},$

where $a_{4}^{(t)}\triangleq
g(\theta^{(t)})+T_{\max}\theta^{(t)}-1/\theta^{(t)}$.

###  4.4 MLEs of parameters in a class of discrete distributions

####  4.4.1 MLE of $\theta$ in zeta distribution

The zeta (or Zipf) distribution can be used to model the frequency of
occurrence of a word randomly chosen from a text, or the population rank of a
city randomly chosen from a country. Let
$\\{X_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle{\rm
iid}}}{{\sim}}\mbox{Zeta}(\theta)$ with pmf $Z^{-1}(\theta)x^{-\theta}$,
$x=1,2,\ldots,\infty$, where $\theta\,(>1)$ is a parameter and
$Z(\theta)=\sum_{x=1}^{\infty}x^{-\theta}$ is called Riemann’s zeta function
(Nair et al. 2018, p.131–134). The log-likelihood is
$\ell(\theta)=-nG({\boldsymbol{x}})\theta-n\log Z(\theta)$, where
$G({\boldsymbol{x}})\triangleq(1/n)\sum_{i=1}^{n}\log(x_{i})$. Let
$0=\ell^{\prime}(\theta)=-nG({\boldsymbol{x}})-nZ^{\prime}(\theta)/Z(\theta)$,
then the MLE of $\theta$ is the root of the equation:
$g(\theta)=-G({\boldsymbol{x}})-Z^{\prime}(\theta)/Z(\theta)=0$ for
$\theta>1$. Note that

$g^{\prime}(\theta)=-\;\frac{Z^{\prime\prime}(\theta)Z(\theta)-[Z^{\prime}(\theta)]^{2}}{Z^{2}(\theta)}\geqslant-\;\frac{Z^{\prime\prime}(\theta)}{Z(\theta)}\stackrel{{\scriptstyle(\rm\ref{eqn4.2})\;\&\;(\rm\ref{eqnfiveC.3})}}{{\geqslant}}a_{5}-\frac{a_{6}}{(\theta-1)^{2}}\;\triangleq\;b(\theta),$

where $a_{5}\triangleq-\log 2\cdot(\log 2+2)$ and $a_{6}\triangleq 2\log 2+2$,
then the US iteration for calculating the MLE $\hat{\theta}$ is given by

$\displaystyle\theta^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})+\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z=0,\;\forall\;\theta,\theta^{(t)}>1\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})+a_{5}(\theta-\theta^{(t)})+\frac{a_{6}}{\theta-1}-\frac{a_{6}}{\theta^{(t)}-1}=0,\;\forall\;\theta,\theta^{(t)}>1\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{a_{5}\theta^{2}+a_{7}^{(t)}\theta+a_{8}^{(t)}=0,\;\forall\;\theta,\theta^{(t)}>1\right\\},$

where $a_{7}^{(t)}\triangleq
g(\theta^{(t)})-a_{5}(\theta^{(t)}+1)-a_{6}/(\theta^{(t)}-1)$ and
$a_{8}^{(t)}\triangleq-g(\theta^{(t)})+a_{5}\theta^{(t)}+a_{6}\theta^{(t)}/(\theta^{(t)}-1)$.

####  4.4.2 MLE of $\theta$ in Yule–Simon distribution

The Yule–Simon distribution (Yule 1925, Garcia 2011) can be employed to model
species among genera, words frequencies in texts, numbers of papers published
by researchers, cities by population, income by size (Leisen, Rossini & Villa
2017). Let $\\{X_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle{\rm
iid}}}{{\sim}}\mbox{YS}(\theta)$ with pmf $\theta\,\mbox{Beta}(x,\theta+1)$,
$x=1,2,\ldots,\infty$, where $\theta\,(>0)$ is the shape parameter. The log-
likelihood function is
$\ell(\theta)=n\log(\theta)+n\log\Gamma(\theta+1)-\sum_{i=1}^{n}\log\Gamma(x_{i}+\theta+1)+\mbox{constant}$.
Thus, the MLE of $\theta$ is the root of the following equation:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\ell^{\prime}(\theta)=\frac{n}{\theta}+\sum_{i=1}^{n}[\psi(\theta+1)-\psi(x_{i}+\theta+1)]$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive4.11})}}{{=}}$
$\displaystyle\frac{n}{\theta}+\sum_{i=1}^{n}\left(-\sum_{m=0}^{\infty}\frac{1}{m+\theta+1}+\sum_{m=0}^{\infty}\frac{1}{m+x_{i}+\theta+1}\right)\qquad[\mbox{Let
}m+x_{i}=k]$ $\displaystyle=$
$\displaystyle\frac{n}{\theta}+\sum_{i=1}^{n}\left(-\sum_{m=0}^{\infty}\frac{1}{m+\theta+1}+\sum_{k=x_{i}}^{\infty}\frac{1}{k+\theta+1}\right)$
$\displaystyle=$
$\displaystyle\frac{n}{\theta}+\sum_{i=1}^{n}\left(-\sum_{m=0}^{x_{i}-1}\frac{1}{m+\theta+1}\right)\;\triangleq\;g(\theta).$

Since

$\displaystyle g^{\prime}(\theta)$ $\displaystyle=$
$\displaystyle-\frac{n}{\theta^{2}}+\sum_{i=1}^{n}\left[\sum_{m=0}^{x_{i}-1}\frac{1}{(m+\theta+1)^{2}}\right]$
$\displaystyle=$
$\displaystyle-\frac{n}{\theta^{2}}+\sum_{i=1}^{n}\left[\frac{1}{(\theta+1)^{2}}+I(x_{i}\geqslant
2)\sum_{m=1}^{x_{i}-1}\frac{1}{(m+\theta+1)^{2}}\right]\geqslant-\frac{n}{\theta^{2}}+\frac{n}{(\theta+1)^{2}}\triangleq
b(\theta),$

the US iteration for calculating the MLE $\hat{\theta}$ is

$\displaystyle\theta^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})+\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z=0,\;\forall\;\theta,\theta^{(t)}>0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})+\frac{n}{\theta}-\frac{n}{\theta+1}-\frac{n}{\theta^{(t)}}+\frac{n}{\theta^{(t)}+1}=0,\;\forall\;\theta,\theta^{(t)}>0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{a_{9}^{(t)}\theta^{2}+a_{9}^{(t)}\theta+n=0,\;\forall\;\theta,\theta^{(t)}>0\right\\},$

where $a_{9}^{(t)}\triangleq
g(\theta^{(t)})-n/\theta^{(t)}+n/(\theta^{(t)}+1)$.

####  4.4.3 MLE of $\alpha$ in Gamma–Poisson distribution

Let $\\{X_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle{\rm
iid}}}{{\sim}}\mbox{GPoisson}(\alpha,\beta)$ with pmf
$\Gamma(x+\alpha)\Gamma^{-1}(\alpha)\beta^{\alpha}(\beta+1)^{-x-\alpha}/x!$,
$x=0,1,2,\ldots,\infty$, where $\alpha$ and $\beta$ are two positive
parameters (Arbous & Kerrich 1951, Bates & Neyman 1952). Define
$\bar{x}=(1/n)\sum_{i=1}^{n}x_{i}$, the log-likelihood function is

$\ell(\alpha,\beta)=\sum_{i=1}^{n}\log\Gamma(x_{i}+\alpha)-n\log\Gamma(\alpha)+n\alpha\log(\beta)-n(\alpha+\bar{x})\log(\beta+1)+\mbox{constant}.$

Let
$0=\partial\ell(\alpha,\beta)/\partial\beta=n\alpha/\beta-n(\alpha+\bar{x})/(\beta+1)$,
we obtain $\beta=\alpha/\bar{x}$. Replacing $\beta$ in $\ell(\alpha,\beta)$ by
$\alpha/\bar{x}$, we know that $\ell(\alpha,\beta)$ becomes

$\ell(\alpha)=\sum_{i=1}^{n}\log\Gamma(x_{i}+\alpha)-n\log\Gamma(\alpha)+n\alpha\log\alpha-\
n(\alpha+\bar{x})\log(\alpha+\bar{x})+\mbox{constant}.$

Thus, the MLE of $\alpha$ is the root of the following equation:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\ell^{\prime}(\alpha)=\sum_{i=1}^{n}\left[\psi(x_{i}+\alpha)-\psi(\alpha)\right]+n\log\alpha-n\log(\alpha+\bar{x})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive4.11})}}{{=}}$
$\displaystyle\sum_{i=1}^{n}\left[-\sum_{m=0}^{\infty}\frac{1}{m+x_{i}+\alpha}+\sum_{m=0}^{\infty}\frac{1}{m+\alpha}\right]+n\log\alpha-n\log(\alpha+\bar{x})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\left[-\sum_{m=x_{i}}^{\infty}\frac{1}{m+\alpha}+\sum_{m=0}^{\infty}\frac{1}{m+\alpha}\right]+n\log\alpha-n\log(\alpha+\bar{x})$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{n}\left[I(x_{i}\geqslant
1)\sum_{m=0}^{x_{i}-1}\frac{1}{m+\alpha}\right]+n\log\alpha-n\log(\alpha+\bar{x})\triangleq
g(\alpha).$

Since

$\displaystyle g^{\prime}(\alpha)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\left[-I(x_{i}\geqslant
1)\sum_{m=0}^{x_{i}-1}\frac{1}{(m+\alpha)^{2}}\right]+\frac{n}{\alpha}-\frac{n}{\alpha+\bar{x}}$
$\displaystyle\geqslant$ $\displaystyle-\sum_{i=1}^{n}\left[I(x_{i}\geqslant
2)\sum_{m=1}^{x_{i}-1}\frac{1}{m^{2}}\right]-\frac{1}{\alpha^{2}}\sum_{i=1}^{n}I(x_{i}\geqslant
1)+\frac{n\bar{x}}{\alpha(\alpha+\bar{x})}$ $\displaystyle\geqslant$
$\displaystyle a_{10}+a_{11}/\alpha^{2}\triangleq\;b(\alpha),$

where $a_{10}\triangleq-\sum_{i=1}^{n}\left[I(x_{i}\geqslant
2)\sum_{m=1}^{x_{i}-1}m^{-2}\right]$ and
$a_{11}\triangleq-\sum_{i=1}^{n}I(x_{i}\geqslant 1)$, the US iteration for
calculating the MLE $\hat{\theta}$ is

$\displaystyle\alpha^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(\alpha^{(t)})+\int_{\alpha^{(t)}}^{\alpha}b(z)\,\mbox{d}z=0,\;\forall\;\alpha,\alpha^{(t)}>0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(\alpha^{(t)})+a_{10}(\alpha-\alpha^{(t)})-a_{11}/\alpha+a_{11}/\alpha^{(t)}=0,\;\forall\;\alpha,\alpha^{(t)}>0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{a_{10}\alpha^{2}+a_{12}^{(t)}\alpha-
a_{11}=0,\;\forall\;\alpha,\alpha^{(t)}>0\right\\},$

where $a_{12}^{(t)}\triangleq
g(\alpha^{(t)})-a_{10}\alpha^{(t)}+a_{11}/\alpha^{(t)}$.

####  4.4.4 MLE of $\theta$ in Generalized Poisson distribution

Let $\\{X_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle{\rm
iid}}}{{\sim}}\mbox{GP}(\lambda,\theta)$ with pmf $\lambda(\lambda+\theta
x)^{x-1}\mbox{e}^{-\lambda-\theta x}/x!$, $x=0,1,2,\ldots,\infty$, where
$\lambda>0$ and $\max(-1,-\lambda/r)<\theta\leqslant 1$ are two parameters and
$r\ (\geqslant 4)$ is the largest positive integer for which $\lambda+\theta
r>0$ when $\theta<0$ (Consul & Jain 1973, Consul & Famoye 1989). The log-
likelihood is
$\ell(\lambda,\theta)=\sum_{i=1}^{n}\left[\log\lambda+(x_{i}-1)\log(\lambda+\theta
x_{i})-\lambda-\theta x_{i}\right]+\mbox{constant}$. There is an MM algorithm
for calculating the MLE of $\theta$ for the case of $\theta\in[0,1]$, while,
to the best of our knowledge, there is no efficient algorithm to calculate
$\hat{\theta}$ for the case of $\theta<0$. In this subsection, we only
consider the case of $\theta<0$. Let

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{\partial\ell(\lambda,\theta)}{\partial\lambda}=\frac{n}{\lambda}+\sum_{i=1}^{n}\frac{x_{i}-1}{\lambda+\theta
x_{i}}-n\quad\mbox{and}\quad$ (4.12) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{\partial\ell(\lambda,\theta)}{\partial\theta}=\sum_{i=1}^{n}\frac{(x_{i}-1)x_{i}}{\lambda+\theta
x_{i}}-n\bar{x}.$ (4.13)

By multiplying $\lambda$ and $\theta$ to the two partial differential
equations (4.12) and (4.13), respectively, and then adding them together, we
obtain $\lambda=(1-\theta)\bar{x}$. Replacing $\lambda$ in
$\ell(\lambda,\theta)$ by $(1-\theta)\bar{x}$, we know that
$\ell(\lambda,\theta)$ will reduce to
$\ell(\theta)=n\log(1-\theta)+\sum_{i=1}^{n}\big{\\{}(x_{i}-1)\log[\bar{x}+(x_{i}-\bar{x})\theta]\big{\\}}+\mbox{constant}$.
Thus, the MLE of $\theta$ is the root of the following equation:

$\displaystyle
0=\ell^{\prime}(\theta)=-\frac{n}{1-\theta}+\sum_{i=1}^{n}\frac{(x_{i}-1)(x_{i}-\bar{x})}{\bar{x}+(x_{i}-\bar{x})\theta}\;\triangleq\;g(\theta).$

Since

$\displaystyle g^{\prime}(\theta)$ $\displaystyle=$
$\displaystyle-\;\frac{n}{(1-\theta)^{2}}-\sum_{i=1}^{n}\frac{(x_{i}-1)(x_{i}-\bar{x})^{2}}{[\bar{x}+(x_{i}-\bar{x})\theta]^{2}}\geqslant-\;\frac{n}{(1-\theta)^{2}}-\sum_{i=1}^{n}\frac{x_{i}(x_{i}-\bar{x})^{2}}{[\bar{x}(1-\theta)+x_{i}\theta]^{2}}$
$\displaystyle\geqslant$
$\displaystyle-\;\frac{n}{(1-\theta)^{2}}-\sum_{i=1}^{n}\frac{x_{i}(x_{i}-\bar{x})^{2}}{[\bar{x}(1-\theta)+x_{(n)}\theta]^{2}}=-\;\frac{n}{(1-\theta)^{2}}-\sum_{i=1}^{n}\frac{x_{i}(x_{i}-\bar{x})^{2}}{[\bar{x}+(x_{(n)}-\bar{x})\theta]^{2}}$
$\displaystyle=$
$\displaystyle-\;\frac{n}{(1-\theta)^{2}}-\frac{\sum_{i=1}^{n}[x_{i}(x_{i}-\bar{x})^{2}/(x_{(n)}-\bar{x})^{2}]}{[\bar{x}/(x_{(n)}-\bar{x})+\theta]^{2}}$
$\displaystyle\triangleq$
$\displaystyle-\;\frac{n}{(1-\theta)^{2}}-\frac{a_{13}}{(a_{14}+\theta)^{2}}\triangleq
b(\theta),$

where $x_{(n)}\triangleq\max(x_{1},\ldots,x_{n})$,
$a_{13}\triangleq\sum_{i=1}^{n}[x_{i}(x_{i}-\bar{x})^{2}/(x_{(n)}-\bar{x})^{2}]$,
and $a_{14}\triangleq\bar{x}/(x_{(n)}-\bar{x})$, the US iteration for
calculating the MLE $\hat{\theta}$ is

$\displaystyle\theta^{(t+1)}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.2})}}{{=}}$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})+\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z=0,\;\forall\;\theta,\theta^{(t)}<0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})-\frac{n}{1-\theta}+\frac{a_{13}}{a_{14}+\theta}+\frac{n}{1-\theta^{(t)}}-\frac{a_{13}}{a_{14}+\theta^{(t)}}=0,\;\forall\;\theta,\theta^{(t)}<0\right\\}$
$\displaystyle=$
$\displaystyle\mbox{sol}\left\\{a_{15}^{(t)}\theta^{2}-a_{16}^{(t)}\theta+(n-a_{15}^{(t)})a_{14}-a_{13}=0,\;\forall\;\theta,\theta^{(t)}<0\right\\},$

where $a_{15}^{(t)}\triangleq
g(\theta^{(t)})+n/(1-\theta^{(t)})-a_{13}/(a_{14}+\theta^{(t)})$ and
$a_{16}^{(t)}\triangleq a_{15}^{(t)}(1-a_{14})-n-a_{13}$.

##  5 . Analysis of the convergence rate

Let $\theta^{*}$ be the unique root of the equation $g(\theta)=0$ for all
$\theta\in\Theta$. To get the rate of convergence for the US iteration
$\theta^{(t+1)}\triangleq h(\theta^{(t)})$, we consider the first–order Taylor
expansion of $h(\theta^{(t)})$ around $\theta^{*}$:
$h(\theta^{(t)})=h(\theta^{*})+(\theta^{(t)}-\theta^{*})h^{\prime}(\theta^{*})+0.5(\theta^{(t)}-\theta^{*})^{2}h^{\prime\prime}(\tilde{\theta})$,
where $\tilde{\theta}$ is a point between $\theta^{(t)}$ and $\theta^{*}$. The
convergence rate of the US algorithm is defined by

$\displaystyle\lim_{t\to\infty}\frac{|\theta^{(t+1)}-\theta^{*}|}{|\theta^{(t)}-\theta^{*}|}$
$\displaystyle=$
$\displaystyle\lim_{t\to\infty}\frac{|h(\theta^{(t)})-h(\theta^{*})|}{|\theta^{(t)}-\theta^{*}|}$
$\displaystyle=$
$\displaystyle\lim_{t\to\infty}|h^{\prime}(\theta^{*})+0.5(\theta^{(t)}-\theta^{*})h^{\prime\prime}(\tilde{\theta})|=|h^{\prime}(\theta^{*})|.$

If $|h^{\prime}(\theta^{*})|\in(0,1)$, then the US algorithm has a linear
convergence rate.

In this section, we investigate the convergence rate of the US algorithm for
the three methods of constructing $U$-functions presented in Section 3.
Theorem 2 below states that the US algorithm holds a linear convergence rate,
just like the EM and MM algorithms with monotonic convergence and linear
convergence rate in maximization of an objective function. Theorems 3–4 below
show that the US algorithm holds a quadratic and even a cubic convergence
rates, corresponding to the SLUB constants method and the TLB constant method.

Theorem 2 (Linear convergence rate for the FLB function method). Let the FLB
function $b(\theta)$ be defined by (3.1), and the $U$-function
$U(\theta|\theta^{(t)})$ be given by (3.2). Then the US iteration

$\displaystyle\theta^{(t+1)}=\mbox{sol}\,\left\\{U(\theta|\theta^{(t)})=g(\theta^{(t)})+\int_{\theta^{(t)}}^{\theta}b(z)\,\mbox{d}z=0,\;\forall\;\theta,\theta^{(t)}\in\Theta\right\\}\;\triangleq\;h(\theta^{(t)})$
(5.1)

holds a linear convergence rate given by

$\displaystyle\lim_{t\to\infty}\frac{|\theta^{(t+1)}-\theta^{*}|}{|\theta^{(t)}-\theta^{*}|}=|h^{\prime}(\theta^{*})|=\left|1-\frac{g^{\prime}(\theta^{*})}{b(\theta^{*})}\right|\in(0,1).$
(5.2)

$\|$

Theorem 3 (Quadratic convergence rate for the SLUB constants method). Let two
SLUB constants $\\{b_{21},b_{22}\\}$ be defined by (3.9), and the $U$-function
$U(\theta|\theta^{(t)})$ be given by (3.11). Then the US iteration

$\displaystyle\theta^{(t+1)}$ $\displaystyle=$
$\displaystyle\mbox{sol}\left\\{g(\theta^{(t)})+g^{\prime}(\theta^{(t)})(\theta-\theta^{(t)})+\frac{b_{2}(\theta|\theta^{(t)})}{2}(\theta-\theta^{(t)})^{2}=0,\quad\forall\;\theta,\theta^{(t)}\in\Theta\right\\}$
(5.3) $\displaystyle\triangleq$ $\displaystyle h(\theta^{(t)})$

has a quadratic convergence rate given by

$\displaystyle\lim_{t\to\infty}\frac{|\theta^{(t+1)}-\theta^{*}|}{|\theta^{(t)}-\theta^{*}|^{2}}=\frac{|h^{\prime\prime}(\theta^{*})|}{2}=\left|\frac{b_{22}-g^{\prime\prime}(\theta^{*})}{2g^{\prime}(\theta^{*})}\right|\in(0,\infty).$
(5.4)

$\|$

Theorem 4 (Cubic convergence rate for the TLB constant method). Let the TLB
constant $b_{3}$ be defined by (3.13), and the $U$-function
$U(\theta|\theta^{(t)})$ be given by (3.14). Then the US iteration

$\displaystyle\theta^{(t+1)}$ $\displaystyle=$
$\displaystyle\mbox{sol}\bigg{\\{}g(\theta^{(t)})+g^{\prime}(\theta^{(t)})(\theta-\theta^{(t)})+\frac{g^{\prime\prime}(\theta^{(t)})}{2}(\theta-\theta^{(t)})^{2}+\frac{b_{3}}{6}(\theta-\theta^{(t)})^{3}=0,$
(5.5)
$\displaystyle\qquad\forall\;\theta,\theta^{(t)}\in\Theta\bigg{\\}}\;\triangleq\;h(\theta^{(t)})$

has a cubic convergence rate given by

$\displaystyle\lim_{t\to\infty}\frac{|\theta^{(t+1)}-\theta^{*}|}{|\theta^{(t)}-\theta^{*}|^{3}}=\frac{|h^{\prime\prime\prime}(\theta^{*})|}{6}=\bigg{|}\frac{b_{3}-g^{\prime\prime\prime}(\theta^{*})}{6g^{\prime}(\theta^{*})}\bigg{|}\in(0,\infty).$
(5.6)

$\|$

##  6 . Numerical experiments

To compare the proposed US algorithm with Newton’s method and bisection method
by using numerical experiments, we consider two root-finding problems: Roots
of a high-order polynomial equation and the $p$-th quantile of a normal
distribution.

###  6.1 Root of a high-order polynomial equation

Assume that there exists a unique real root for the following polynomial
equation in $(0,a)$:

$\displaystyle 0=g(\theta)\triangleq
a_{3}\theta^{m}+a_{2}\theta^{2}+a_{1}\theta+a_{0}\triangleq
g_{1}(\theta)+a_{2}\theta^{2}+a_{1}\theta+a_{0},\quad\theta\in(0,a),$ (6.1)

where $m\;(\geqslant 3)$ is a real number (unnecessary positive integer) and
$a_{0}>0$. Based on

$\displaystyle g_{1}^{\prime\prime}(\theta)$ $\displaystyle=$ $\displaystyle
a_{3}m(m-1)\theta^{m-2}\;\in\;(a_{3}m(m-1)a^{m-2},0),\quad\mbox{if }\
a_{3}<0,$ $\displaystyle g_{1}^{\prime\prime\prime}(\theta)$ $\displaystyle=$
$\displaystyle
a_{3}m(m-1)(m-2)\theta^{m-3}\;>\;0,\quad\quad\quad\quad\quad\quad\mbox{if }\
a_{3}>0,$

and Property 1(c), we can construct the following two quadratic $U$-functions:

$\displaystyle U_{2}(\theta|\theta^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.10})}}{{=}}$
$\displaystyle
a_{3}[\theta^{(t)}]^{m}+a_{3}m[\theta^{(t)}]^{m-1}(\theta-\theta^{(t)})+0.5b_{2}(\theta|\theta^{(t)})(\theta-\theta^{(t)})^{2}$
(6.2) $\displaystyle+\;a_{2}\theta^{2}+a_{1}\theta+a_{0}$ $\displaystyle=$
$\displaystyle\Big{[}0.5b_{2}(\theta|\theta^{(t)})+a_{2}\Big{]}\theta^{2}+\Big{\\{}a_{3}m[\theta^{(t)}]^{m-1}-b_{2}(\theta|\theta^{(t)})\theta^{(t)}+a_{1}\Big{\\}}\theta$
$\displaystyle+\;a_{3}(1-m)[\theta^{(t)}]^{m}+0.5b_{2}(\theta|\theta^{(t)})[\theta^{(t)}]^{2}+a_{0}\quad\mbox{and}\quad$
$\displaystyle U_{3}(\theta|\theta^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.12})}}{{=}}$
$\displaystyle
a_{3}[\theta^{(t)}]^{m}+a_{3}m[\theta^{(t)}]^{m-1}(\theta-\theta^{(t)})+0.5a_{3}m(m-1)[\theta^{(t)}]^{m-2}(\theta-\theta^{(t)})^{2}$
(6.3) $\displaystyle+\;a_{2}\theta^{2}+a_{1}\theta+a_{0}$ $\displaystyle=$
$\displaystyle\Big{\\{}0.5a_{3}m(m-1)[\theta^{(t)}]^{m-2}+a_{2}\Big{\\}}\theta^{2}+\Big{\\{}a_{3}m(2-m)[\theta^{(t)}]^{m-1}+a_{1}\Big{\\}}\theta$
$\displaystyle+\;a_{3}(m-1)(0.5m-1)[\theta^{(t)}]^{m}+a_{0},$

where $b_{2}(\theta|\theta^{(t)})=0\cdot
I(\theta\leqslant\theta^{(t)})+a_{3}m(m-1)a^{m-2}\cdot
I(\theta>\theta^{(t)})$.

For the purpose of numerical experiments, we set
$\\{a_{0},a_{1},a_{2},a_{3},m\\}=\\{1,-1,1,-1,3\\}$ in (6.2),
$\\{a_{0},a_{1},a_{2},a_{3},m\\}=\\{1,-1,-3,1,3\\}$ in (6.3), and $a=2$. We
know that the unique root of (6.1) in the interval (0, 2) is
$\theta^{*}=0.4608111$ and $\theta^{*}=1$ for the two cases. We use three
algorithms (i.e., the US algorithm, Newton’s method and bisection method) to
seek the root of (6.1) by running each algorithm for 100,000 repetitions with
initial values randomly chosen from the uniform distribution on the interval
$(0,2)$. The percentage of converged algorithms among 100,000 repetitions
(denoted by Percentage), the average number of iterations for converged
algorithms among 100,000 repetitions (denoted by N), and the total computation
time in second of calculating $\theta^{*}$ for 100,000 repetitions [denoted by
Time(s)] are recorded and reported in Table 2.

From Table 2, we can see that all three algorithms have 100% convergence
performance according to Percentage. The US algorithm and Newton’s methods
have similar performance in terms of N and Time(s). The bisection method
performs worst among the three for the two cases.

Table 2. Comparisons of three algorithms in terms of the percentage of
converged algorithms, the average number of iterations, and the total
computation time in second

Coefficients | $\\{a_{0},a_{1},a_{2},a_{3},m\\}=\\{1,-1,1,-1,3\\}$ | $\\{a_{0},a_{1},a_{2},a_{3},m\\}=\\{1,-1,-3,1,3\\}$
---|---|---
Algorithms | $\mbox{US}_{2}$ | Newton | Bisection | $\mbox{US}_{3}$ | Newton | Bisection
Percentage | 100.0% | 100.0% | 100.0% | 100.0% | 100.0% | 100.0%
N | 7.0000 | 6.3672 | 41.000 | 7.0000 | 6.3630 | 40.000
Time(s) | 1.6931 | 1.5313 | 10.255 | 1.6931 | 1.6206 | 10.184

Note: $\mbox{US}_{2}$ denotes the US algorithm based on
$U_{2}(\theta|\theta^{(t)})$ in (6.2); $\mbox{US}_{3}$ denotes the US
algorithm based on $U_{3}(\theta|\theta^{(t)})$ in (6.3); Percentage is the
percentage of converged algorithms among 100,000 repetitions for different
initial values randomly chosen from $U(0,2)$; N is the average number of
iterations for converged algorithms among 100,000 repetitions; Time(s) is the
total computation time in second of calculating $\theta^{*}$ for 100,000
repetitions.

###  6.2 Calculation of the $p$-th quantile of normal distribution

In Example 2 of subsection 4.1.1, we have applied the US algorithm to find the
$p$-th quantile of the normal distribution $N(\mu,\sigma^{2})$; i.e., to
iteratively calculate the solution to the integral equation $0=g(x)\triangleq
p-\Phi(x|\mu,\sigma^{2})$ with US iteration given by (4.4). Theorem 2 shows
that the US iteration (4.4) only has a linear rate of convergence. To speed
the convergence, we first derive the second–derivative to the
fourth–derivative of $g(x)$ as follows:

$\displaystyle g^{\prime\prime}(x)$ $\displaystyle=$
$\displaystyle\frac{x-\mu}{\sigma^{2}}\phi(x|\mu,\sigma^{2}),$ $\displaystyle
g^{\prime\prime\prime}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma^{2}}\left[1-\frac{(x-\mu)^{2}}{\sigma^{2}}\right]\phi(x|\mu,\sigma^{2})\quad\mbox{and}\quad$
$\displaystyle g^{(4)}(x)$ $\displaystyle=$
$\displaystyle-\frac{x-\mu}{\sigma^{4}}\left[3-\frac{(x-\mu)^{2}}{\sigma^{2}}\right]\phi(x|\mu,\sigma^{2}).$

Note that there are three roots $\\{\mu,\mu\pm\sqrt{3}\sigma\\}$ for the
equation $0=g^{(4)}(x)$, then $g^{\prime\prime\prime}(x)$ achieves its minimum
$-2(\sqrt{2\pi}\sigma^{3}\mbox{e}^{3/2})^{-1}$ when $x=\mu\pm\sqrt{3}\sigma$.
Similarly, since $\mu\pm\sigma$ are two roots of the equation
$0=g^{\prime\prime\prime}(x)$, we know that $g^{\prime\prime}(x)$ achieves its
minimum $-(\sqrt{2\pi}\sigma^{2}\mbox{e}^{1/2})^{-1}$ and maximum
$(\sqrt{2\pi}\sigma^{2}\mbox{e}^{1/2})^{-1}$ when $x=\mu-\sigma$ and
$\mu+\sigma$, respectively. Thus, we obtain

$\displaystyle g^{\prime\prime}(x)$ $\displaystyle\in$
$\displaystyle\left[-(\sqrt{2\pi}\sigma^{2}\mbox{e}^{1/2})^{-1},\;(\sqrt{2\pi}\sigma^{2}\mbox{e}^{1/2})^{-1}\right]\quad\mbox{and}\quad$
$\displaystyle g^{\prime\prime\prime}(x)$ $\displaystyle\geqslant$
$\displaystyle-2(\sqrt{2\pi}\sigma^{3}\mbox{e}^{3/2})^{-1},\quad\forall\;x\in{\mathbb{R}},$

so that we can construct the following two $U$-functions:

$\displaystyle U_{2}(x|x^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.10})}}{{=}}$
$\displaystyle
g(x^{(t)})+g^{\prime}(x^{(t)})(x-x^{(t)})+0.5b_{2}(x|x^{(t)})(x-x^{(t)})^{2}\quad\mbox{and}\quad$
(6.4) $\displaystyle U_{3}(x|x^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.12})}}{{=}}$
$\displaystyle
g(x^{(t)})+g^{\prime}(x^{(t)})(x-x^{(t)})+0.5g^{\prime\prime}(x^{(t)})(x-x^{(t)})^{2}$
(6.5)
$\displaystyle-\;(3\sqrt{2\pi}\sigma^{3}\mbox{e}^{3/2})^{-1}(x-x^{(t)})^{3},$

where
$b_{2}(x|x^{(t)})=(\sqrt{2\pi}\sigma^{2}\mbox{e}^{1/2})^{-1}\times[I(x\leqslant
x^{(t)})-I(x>x^{(t)})]$.

In numerical experiments, we set $p=0.01,\,0.9$, $\mu=-2,\,2$ and $\sigma=1$.
For each combination of the four cases of $\\{p,\mu\\}$, we use four
algorithms (i.e., two US algorithms based on (6.4) and (6.5), Newton’s method
and bisection method) to seek the root of (6.1) by running each algorithm for
100,000 repetitions with initial values randomly chosen from $U(-4,4)$. The
percentage of converged algorithms among 100,000 repetitions (denoted by
Percentage), the average number of iterations for converged algorithms among
100,000 repetitions (denoted by N), and the total computation time in second
of calculating $x^{*}$ for 100,000 repetitions [denoted by Time(s)] are
recorded and displayed in Table 3.

In Table 3, we can see that Newton’s method has the worst convergence
performance according to Percentage. In terms of N, the $\mbox{US}_{3}$
performs the best. In terms of Time(s), although Newton’s method only takes
about 2 seconds for 100,000 repetitions for Case 1, we can see that there are
only about 40% valid computations, hence we prefer $\mbox{US}_{2}$.

Table 3. Comparisons of four algorithms in terms of the percentage of
converged algorithms, the average number of iterations, and the total
computation time in second

Parameters | $\\{p,\mu\\}=\\{0.01,-2\\}$ | $\\{p,\mu\\}=\\{0.01,2\\}$
---|---|---
Algorithms | $\mbox{US}_{2}$ | $\mbox{US}_{3}$ | Newton | Bisection | $\mbox{US}_{2}$ | $\mbox{US}_{3}$ | Newton | Bisection
Percentage | 100.0% | 100.0% | 39.98% | 100.0% | 100.0% | 100.0% | 55.40% | 100.0%
N | 10.544 | 5.3599 | 6.6532 | 33.000 | 10.2855 | 5.0179 | 6.3664 | 36.000
Time(s) | 7.7039 | 14.560 | 2.2438 | 12.678 | 7.62123 | 14.026 | 2.4824 | 13.852
Parameters | $\\{p,\mu\\}=\\{0.9,-2\\}$ | $\\{p,\mu\\}=\\{0.9,2\\}$
Algorithms | $\mbox{US}_{2}$ | $\mbox{US}_{3}$ | Newton | Bisection | $\mbox{US}_{2}$ | $\mbox{US}_{3}$ | Newton | Bisection
Percentage | 100.0% | 100.0% | 42.74% | 100.0% | 100.0% | 100.0% | 38.65% | 100.0%
N | 5.9959 | 4.1400 | 5.0428 | 39.000 | 6.7453 | 4.5215 | 4.9311 | 40.000
Time(s) | 4.7217 | 10.066 | 2.0309 | 14.832 | 5.2619 | 10.515 | 1.9345 | 15.421

Note: $\mbox{US}_{2}$ denotes the US algorithm based on
$U_{2}(\theta|\theta^{(t)})$ in (6.4); $\mbox{US}_{3}$ denotes the US
algorithm based on $U_{3}(\theta|\theta^{(t)})$ in (6.5); Percentage is the
percentage of converged algorithms among 100,000 repetitions for different
initial values randomly chosen from $U(-4,4)$; N is the average number of
iterations for converged algorithms among 100,000 repetitions; Time(s) is the
total computation time in second of calculating $x^{*}$ for 100,000
repetitions.

###  6.3 Solving an equation with multiple roots

Thanks for the strongly stable convergence, the US algorithm could be one of
the powerful tools for solving an equation with multiple roots. For the
purpose of illustration, we consider the following toy equation:

$\displaystyle 0=g(x)\triangleq-0.5x-2\sin(x)+1,\quad x\in{\mathbb{R}}.$ (6.6)

Figure 2(a) shows the equation (6.6) with three roots, denoted by
$x_{1}^{*}<x_{2}^{*}<x_{3}^{*}$, respectively. First, based on the first–order
derivative $g^{\prime}(x)=-0.5-2\cos(x)\geqslant-2.5$, we can construct the
first $U$-function for $g(x)$ at $x=x^{(t)}$ and obtain the corresponding US
iterations as follows:

$\displaystyle U_{1,1}(x|x^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.5})}}{{=}}$ $\displaystyle
g(x^{(t)})-2.5(x-x^{(t)})\quad\mbox{and}\quad$ $\displaystyle x^{(t+1)}$
$\displaystyle=$ $\displaystyle
x^{(t)}+\frac{g(x^{(t)})}{2.5}=x^{(t)}+\frac{-0.5x^{(t)}-2\sin(x^{(t)})+1}{2.5}.$
(6.7)

Given an initial value $x^{(0)}\,(<x_{1}^{*})$, say $x^{(0)}=0$, the US
iterations $\\{x^{(t+1)}\\}_{t=0}^{\infty}$ defined by (6.3) converged to the
first root $x_{1}^{*}=0.4090497$ as shown in Figure 2(b1).

Second, once we obtained $x_{1}^{*}$, we define a new function

$g_{2}(x)\triangleq-g(x)=0.5x+2\sin(x)-1,\quad x\in{\mathbb{R}}.$

Based on the first–order derivative
$g^{\prime}_{2}(x)=0.5+2\cos(x)\geqslant-1.5$, we can construct the second
$U$-function for $g_{2}(x)$ at $x=x^{(t)}$ and obtain the corresponding US
iterations as follows:

$\displaystyle U_{1,2}(x|x^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive3.5})}}{{=}}$ $\displaystyle
g_{2}(x^{(t)})-1.5(x-x^{(t)})\quad\mbox{and}\quad$ $\displaystyle x^{(t+1)}$
$\displaystyle=$ $\displaystyle
x^{(t)}+\frac{g_{2}(x^{(t)})}{1.5}=x^{(t)}+\frac{0.5x^{(t)}+2\sin(x^{(t)})+1}{2.5}.$
(6.8)

Set $x^{(0)}\triangleq x_{1}^{*}+\varepsilon$ as the initial value such that
$g_{2}(x^{(0)})>0$, where $\varepsilon$ is a small positive real number, say
$10^{-6}$. Then the US iterations $\\{x^{(t+1)}\\}_{t=0}^{\infty}$ defined by
(6.3) converged to the second root $x_{2}^{*}=3.535612$ as shown in Figure
2(b2).

Third, once we obtained $x_{2}^{*}$, we set $x^{(0)}\triangleq
x_{2}^{*}+\varepsilon$ as the initial value such that $g(x^{(0)})>0$, then the
US iterations $\\{x^{(t+1)}\\}_{t=0}^{\infty}$ defined by (6.3) converged to
the second root $x_{3}^{*}=5.308993$ as shown in Figure 2(b3).

(a)
(b)

Figure 2. (a) Plots of $0=g(x)=-0.5x-2\sin(x)+1$ with three roots
$x_{1}^{*}=0.4090497$, $x_{2}^{*}=3.535612$ and $x_{3}^{*}=5.308993$. (b) With
the initial values $x^{(0)}=0$, $x^{(0)}=x_{1}^{*}+10^{-6}$ and
$x^{(0)}=x_{2}^{*}+10^{-6}$, the sequence $\\{x^{(t+1)}\\}_{t=0}^{\infty}$
generated respectively by (6.3), (6.3), (6.3) strongly stably converges to
$\\{x_{k}^{*}\\}_{k=1}^{3}$, satisfying
$x^{(0)}<x^{(1)}<\cdots<x^{(t)}<\cdots\leqslant x_{k}^{*}$ for $k=1,2,3$.

##  7 . Discussions

To solve the root of a non-linear equation $g(\theta)=0$, in this paper, we
have established a general framework of a new root–finding method, called as
the US algorithm. Its each iteration consists of a U-step and an S-step, where
the U-step is to construct one $U$-function $U(\theta|\theta^{(t)})$ based on
a new notion of changing direction inequality for the original function
$g(\theta)$, while the S-step solves the $U$-equation
$U(\theta|\theta^{(t)})=0$ with respect to the left argument to obtain the
next iterate $\theta^{(t+1)}$. Unlike the NR algorithm, which is sensitive to
initial values and may diverge if a poor initial value is chosen, from Theorem
1, we know that the US algorithm has the property of strongly stable
convergence; that is, it does not depend on any initial values in $\Theta$ and
strongly stably converges to the root $\theta^{*}$. Especially, because of the
property of strongly stable convergence, the US algorithm could be one of the
powerful tools for solving an equation with multiple roots as shown in Section
6.3.

The critical requirement for applying the US algorithm is to construct one
$U$-function $U(\theta|\theta^{(t)})$ for $g(\theta)$ at $\theta=\theta^{(t)}$
(i.e., to satisfy the CD inequality
$g(\theta)-U(\theta|\theta^{(t)})\stackrel{{\scriptstyle{\rm
sgn}(\theta-\theta^{(t)})}}{{\geqslant}}0$) such that an explicit solution to
the surrogate equation $U(\theta|\theta^{(t)})=0$ is available. By using
Taylor expansion, we have developed three methods for constructing
$U$-functions based on the FLB function, SLUB constants and TLB constant. The
resulting US algorithms have linear/quadratic/cubic convergence rates as shown
in Theorems 2–4, respectively. In the future studies, we will develop other
methods for constructing $U$-functions. Although the US algorithm enjoys
strongly stable convergence, we would like to accelerate the US algorithm with
slow convergence rate if the bound functions/constants are not tight enough,
particularly for the FLB function-based US algorithms. Some accelerating
techniques for the proposed US algorithms may significantly spark the method’s
potential.

In Section 4, we successfully applied the US algorithms to deal with a wide
range of statistical issues that are not perfectly solved by the traditional
numerical algorithms. Beyond those problems, we are trying to use the US
algorithms to handle other attractive statistical models that are still
remained to be conquered.

## Acknowledgements

Guo-Liang TIAN’s research was partially supported by National Natural Science
Foundation of China (No. 12171225).

## References

Arbous AG and Kerrich JE (1951). Accident statistics and the concept of
accident–proneness. Biometrics 7(4): 340–432.

Azzalini A (1985). A class of distributions which includes the normal ones.
Scandinavian Journal of Statistics 12(2): 171–178.

Bates GE and Neyman J (1952). Contributions to the theory of accident
proneness. 1. an optimistic model of the correlation between light and severe
accidents. University of California Publications in Statistics 1(9): 215–254.

Boyd S and Vandenberghe L (2004). Convex Optimization. Cambridge: Cambridge
University Press.

Chabert JL (1999). A History of Algorithms: From the Pebble to the Microchip,
Chapter 3: Methods of false position. Berlin: Springer.

Consul PC and Famoye F (1989). The truncated generalized Poisson distribution
and its estimation. Communications in Statistics — Theory and Methods 18(10):
3635–3648.

Consul PC and Jain GC (1973). A generalization of the Poisson distribution.
Technometrics 15(4): 791–799.

Costabile F, Gualtieri MI and Luceri R (2006). A modification of Muller’s
method. Calcolo, 43(1): 39–50.

Dempster AP, Laird NM and Rubin DB (1977). Maximum likelihood from incomplete
data via the EM algorithm (with discussions). Journal of the Royal Statistical
Society: Series B (Methodological) 39: 1–38.

Garcia JMG (2011). A fixed–point algorithm to estimate the Yule–Simon
distribution parameter. Applied Mathematics and Computation 217(21):
8560–8566.

Hallinan JAJ (1993). A review of the Weibull distribution. Journal of Quality
Technology 25(2): 85–93.

Hunter DR and Lange K (2004). A tutorial on MM algorithms. The American
Statistician 58(1): 30–37.

Lange K, Hunter DR and Yang I (2000). Optimization transfer using surrogate
objective functions (with discussions). Journal of Computational and Graphical
Statistics 9, 1–20.

Leisen F, Rossini L and Villa C (2017). A note on the posterior inference for
the Yule–Simon distribution. Journal of Statistical Computation and Simulation
87(6), 1179–1188.

Nair U, Sankaran PG and Balakrishnan N (2018). Reliability Modelling and
Analysis in Discrete Time Chapter 3: Discrete Lifetime Models. London:
Academic Press.

Rinne H (2008). The Weibull Distribution: A Handbook. Florida: CRC Press.

Wolfe P (1959). The secant method for simultaneous nonlinear equations.
Communications of the ACM 2(12): 12–13.

Wood GR (1992). The bisection method in higher dimensions. Mathematical
Programming 55: 319–337.

Yule GU (1925). A mathematical theory of evolution, based on the conclusions
of Dr. JC Willis, FR S. Philosophical Transactions of the Royal Society of
London. Series B, Containing Papers of A Biological Character 213(402–410):
21–87.

## Appendix A: Proofs of Theorem 1–4

Proof of Theorem 1. (i) We first prove that the US algorithm strongly stably
converges to the unique root $\theta^{*}$. Let $\theta^{(t)}$ be the $t$-th
iteration of $\theta^{*}$, then we only need to prove that

$\displaystyle\theta^{(t)}<\theta^{(t+1)}\leqslant\theta^{*},$
$\displaystyle\mbox{if }\ \theta^{(t)}<\theta^{*},$ (A.1)
$\displaystyle\theta^{*}\leqslant\theta^{(t+1)}<\theta^{(t)},$
$\displaystyle\mbox{if }\ \theta^{(t)}>\theta^{*}.$ (A.2)

For the case of $\theta^{(t)}<\theta^{*}$, note that $U(\theta|\theta^{(t)})$
is a $U$-function for $g(\theta)$ at $\theta=\theta^{(t)}$, then we have

$\displaystyle U(\theta|\theta^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive2.2})\,\&\,(\ref{eqnfive2.3})}}{{\geqslant}}$
$\displaystyle
g(\theta)\stackrel{{\scriptstyle(\ref{eqnfive2.1})}}{{>}}0,\quad\mbox{if }\
\theta\leqslant\theta^{(t)}(<\theta^{*}),$ (A.3) $\displaystyle
U(\theta|\theta^{(t)})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eqnfive2.4})}}{{\leqslant}}$
$\displaystyle g(\theta),\hskip 41.25641pt\mbox{if }\ \theta>\theta^{(t)}.$
(A.4)

Let $\theta^{(t+1)}$ be the $(t+1)$-th US iteration; i.e.,
$U(\theta^{(t+1)}|\theta^{(t)})=0$. Then we can show

$\displaystyle\theta^{(t+1)}>\theta^{(t)}.$ (A.5)

In fact, if (A.5) is not true (i.e., $\theta^{(t+1)}\leqslant\theta^{(t)}$ is
true), by replacing $\theta$ in (A.3) with $\theta^{(t+1)}$, we have
$0=U(\theta^{(t+1)}|\theta^{(t)})\geqslant g(\theta^{(t+1)})>0$, which is
contradictious. So (A.5) is true.

Next, we can show that

$\displaystyle\theta^{(t+1)}\leqslant\theta^{*}.$ (A.6)

In fact, if (A.6) is not true (i.e., $\theta^{(t+1)}>\theta^{*}$ is true), by
replacing $\theta$ in (A.4) with $\theta^{(t+1)}$, we have
$g(\theta^{(t+1)})\geqslant U(\theta^{(t+1)}|\theta^{(t)})=0$, which is
contradictory to the fact that $g(\theta)<0$ for any $\theta>\theta^{*}$ as
shown by (2.1). So (A.6) is true. By combining (A.5) with (A.6), we obtain
(A.1). For the case of $\theta^{(t)}>\theta^{*}$, similarly, we can prove
(A.2).

(ii) From (A.1) and (A.2), we can see that the US algorithm does not depend on
any initial values in $\Theta$. $\Box$

Proof of Theorem 2. From (5.1), we have

$0=U(h(\theta^{(t)})|\theta^{(t)})=g(\theta^{(t)})+\int_{\theta^{(t)}}^{h(\theta^{(t)})}b(z)\,\mbox{d}z.$

By taking derivative with respect to $\theta^{(t)}$ on the both sides of the
above equation, we obtain
$0=g^{\prime}(\theta^{(t)})+b(h(\theta^{(t)}))h^{\prime}(\theta^{(t)})-b(\theta^{(t)})$,
so that

$h^{\prime}(\theta^{(t)})=\frac{b(\theta^{(t)})-g^{\prime}(\theta^{(t)})}{b(h(\theta^{(t)}))}.$

Note that the rate of convergence for the US algorithm is

$\displaystyle|h^{\prime}(\theta^{*})|$ $\displaystyle=$
$\displaystyle\lim_{t\to\infty}|h^{\prime}(\theta^{(t)})|=\lim_{t\to\infty}\left|\frac{b(\theta^{(t)})-g^{\prime}(\theta^{(t)})}{b(h(\theta^{(t)}))}\right|$
$\displaystyle=$
$\displaystyle\left|\frac{b(\theta^{*})-g^{\prime}(\theta^{*})}{b(h(\theta^{*}))}\right|=\left|\frac{b(\theta^{*})-g^{\prime}(\theta^{*})}{b(\theta^{*})}\right|=\left|1-\frac{g^{\prime}(\theta^{*})}{b(\theta^{*})}\right|,$

we only need to show that

$\displaystyle|1-g^{\prime}(\theta^{*})/b(\theta^{*})|<1.$ (A.7)

If the condition, $b(\theta^{*})\leqslant g^{\prime}(\theta^{*})<0$, holds,
then it is easy to check the correctness of (A.7). From (3.1), we of course
have $b(\theta^{*})\leqslant g^{\prime}(\theta^{*})$.

Now, we only need to prove $g^{\prime}(\theta^{*})<0$. For any
$\varepsilon>0$, from the assumption (2.1), we have
$g(\theta^{*}+\varepsilon)<0$ and $g(\theta^{*}-\varepsilon)>0$ so that we
obtain

$\displaystyle g_{+}^{\prime}(\theta^{*})$ $\displaystyle\triangleq$
$\displaystyle\lim_{\varepsilon\to
0^{+}}\frac{g(\theta^{*}+\varepsilon)-g(\theta^{*})}{\varepsilon}=\lim_{\varepsilon\to
0^{+}}\frac{g(\theta^{*}+\varepsilon)}{\varepsilon}<0\quad\mbox{and}\quad$
$\displaystyle g_{-}^{\prime}(\theta^{*})$ $\displaystyle\triangleq$
$\displaystyle\lim_{\varepsilon\to
0^{+}}\frac{g(\theta^{*}-\varepsilon)-g(\theta^{*})}{-\varepsilon}=\lim_{\varepsilon\to
0^{+}}\frac{g(\theta^{*}-\varepsilon)}{-\varepsilon}<0.$

Since the first–order derivative $g^{\prime}(\cdot)$ exists, we have
$g^{\prime}(\theta^{*})=g_{+}^{\prime}(\theta^{*})=g_{-}^{\prime}(\theta^{*})<0$.
$\Box$

Proof of Theorem 3. From (5.3), we have

$\displaystyle
g(\theta^{(t)})+g^{\prime}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}]+\frac{b_{2}(h(\theta^{(t)}),\theta^{(t)})}{2}[h(\theta^{(t)})-\theta^{(t)}]^{2}=0.$
(A.8)

By taking the first-derivative about $\theta^{(t)}$ on the both sides of
(A.8), we obtain

$g^{\prime}(\theta^{(t)})h^{\prime}(\theta^{(t)})+\Big{\\{}g^{\prime\prime}(\theta^{(t)})+b_{2}(h(\theta^{(t)}),\theta^{(t)})[h^{\prime}(\theta^{(t)})-1]\Big{\\}}[h(\theta^{(t)})-\theta^{(t)}]=0.$

By letting $t\to\infty$, we have
$g^{\prime}(\theta^{*})h^{\prime}(\theta^{*})=0$ because

$\displaystyle\lim_{t\to\infty}[h(\theta^{(t)})-\theta^{(t)}]=h(\theta^{*})-\theta^{*}=0.$
(A.9)

Usually, we know that $g^{\prime}(\theta^{*})\neq 0$ so that

$\displaystyle h^{\prime}(\theta^{*})=0.$ (A.10)

By taking the second-derivative with respect to $\theta^{(t)}$ on the both
sides of (A.8), we have

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
g^{\prime\prime}(\theta^{(t)})h^{\prime}(\theta^{(t)})+g^{\prime}(\theta^{(t)})h^{\prime\prime}(\theta^{(t)})+\Big{\\{}g^{\prime\prime\prime}(\theta^{(t)})+b_{2}(h(\theta^{(t)}),\theta^{(t)})h^{\prime\prime}(\theta^{(t)})\Big{\\}}[h(\theta^{(t)})-\theta^{(t)}]$
$\displaystyle+\;\Big{\\{}g^{\prime\prime}(\theta^{(t)})+b_{2}(h(\theta^{(t)}),\theta^{(t)})[h^{\prime}(\theta^{(t)})-1]\Big{\\}}[h^{\prime}(\theta^{(t)})-1].$

By letting $t\to\infty$ and using (A.9) & (A.10), we obtain

$|h^{\prime\prime}(\theta^{*})|=\lim_{t\to\infty}|h^{\prime\prime}(\theta^{(t)})|=\left|\frac{b_{2}(\theta^{*},\theta^{*})-g^{\prime\prime}(\theta^{*})}{g^{\prime}(\theta^{*})}\right|\;\stackrel{{\scriptstyle(\ref{eqnfive3.9})}}{{=}}\;\left|\frac{b_{22}-g^{\prime\prime}(\theta^{*})}{g^{\prime}(\theta^{*})}\right|,$

indicating that the US iteration (5.3) has a quadratic convergence rate given
by (5.4). $\Box$

Proof of Theorem 4. From (5.5), we have

$\displaystyle
0=g(\theta^{(t)})+g^{\prime}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}]+\frac{g^{\prime\prime}(\theta^{(t)})}{2}[h(\theta^{(t)})-\theta^{(t)}]^{2}+\frac{b_{3}}{6}[h(\theta^{(t)})-\theta^{(t)}]^{3}.\quad$
(A.11)

By taking the first-derivative about $\theta^{(t)}$ on the both sides of
(A.11), we obtain

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
g^{\prime}(\theta^{(t)})h^{\prime}(\theta^{(t)})+g^{\prime\prime}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}]+\frac{g^{\prime\prime\prime}(\theta^{(t)})}{2}[h(\theta^{(t)})-\theta^{(t)}]^{2}$
$\displaystyle+\;g^{\prime\prime}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}][h^{\prime}(\theta^{(t)})-1]+\frac{b_{3}}{2}[h(\theta^{(t)})-\theta^{(t)}]^{2}[h^{\prime}(\theta^{(t)})-1].$

By letting $t\to\infty$, we have $h^{\prime}(\theta^{*})=0$ when
$g^{\prime}(\theta^{*})\neq 0$, because

$\displaystyle\lim_{t\to\infty}[h(\theta^{(t)})-\theta^{(t)}]=h(\theta^{*})-\theta^{*}=0.$
(A.12)

By taking the second-derivative about $\theta^{(t)}$ on the both sides of
(A.11), we obtain

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
g^{\prime\prime}(\theta^{(t)})h^{\prime}(\theta^{(t)})+g^{\prime}(\theta^{(t)})h^{\prime\prime}(\theta^{(t)})+g^{\prime\prime\prime}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}]+g^{\prime\prime}(\theta^{(t)})[h^{\prime}(\theta^{(t)})-1]$
$\displaystyle+\;\frac{g^{(4)}(\theta^{(t)})}{2}[h(\theta^{(t)})-\theta^{(t)}]^{2}+g^{\prime\prime\prime}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}][h^{\prime}(\theta^{(t)})-1]$
$\displaystyle+\;\left\\{g^{\prime\prime\prime}(\theta^{(t)})+\frac{b_{3}}{2}[h^{\prime}(\theta^{(t)})-1]\right\\}\cdot[h(\theta^{(t)})-\theta^{(t)}][h^{\prime}(\theta^{(t)})-1]$
$\displaystyle+\;\left\\{g^{\prime\prime}(\theta^{(t)})+\frac{b_{3}}{2}[h(\theta^{(t)})-\theta^{(t)}]\right\\}\cdot\Big{\\{}[h^{\prime}(\theta^{(t)})-1]^{2}+[h(\theta^{(t)})-\theta^{(t)}]h^{\prime\prime}(\theta^{(t)})\Big{\\}}.$

By letting $t\to\infty$, we have $h^{\prime\prime}(\theta^{*})=0$, because
$h^{\prime}(\theta^{*})=0$ and (A.12). By taking the third-derivative about
$\theta^{(t)}$ on the both sides of equation (A.11), we obtain

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
g^{\prime\prime\prime}(\theta^{(t)})h^{\prime}(\theta^{(t)})+g^{\prime\prime}(\theta^{(t)})h^{\prime\prime}(\theta^{(t)})+g^{\prime\prime}(\theta^{(t)})h^{\prime\prime}(\theta^{(t)})+g^{\prime}(\theta^{(t)})h^{\prime\prime\prime}(\theta^{(t)})$
$\displaystyle+\;g^{(4)}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}]+g^{\prime\prime\prime}(\theta^{(t)})[h^{\prime}(\theta^{(t)})-1]+g^{\prime\prime\prime}(\theta^{(t)})[h^{\prime}(\theta^{(t)})-1]$
$\displaystyle+\;g^{\prime\prime}(\theta^{(t)})h^{\prime\prime}(\theta^{(t)})+\frac{1}{2}g^{(5)}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}]^{2}+g^{(4)}(\theta^{(t)})[h(\theta^{(t)})-\theta^{(t)}][h^{\prime}(\theta^{(t)})-1]$
$\displaystyle+\;\left\\{2g^{(4)}(\theta^{(t)})+\frac{b_{3}}{2}h^{\prime\prime}(\theta^{(t)})\right\\}\cdot[h(\theta^{(t)})-\theta^{(t)}][h^{\prime}(\theta^{(t)})-1]$
$\displaystyle+\;\left\\{2g^{\prime\prime\prime}(\theta^{(t)})+\frac{b_{3}}{2}[h^{\prime}(\theta^{(t)})-1]\right\\}\cdot\left\\{[h^{\prime}(\theta^{(t)})-1]^{2}+[h(\theta^{(t)})-\theta^{(t)}]h^{\prime\prime}(\theta^{(t)})\right\\}$
$\displaystyle+\;\left\\{g^{\prime\prime\prime}(\theta^{(t)})+\frac{b_{3}}{2}[h^{\prime}(\theta^{(t)})-1]\right\\}\cdot\left\\{[h^{\prime}(\theta^{(t)})-1]^{2}+[h(\theta^{(t)})-\theta^{(t)}]h^{\prime\prime}(\theta^{(t)})\right\\}$
$\displaystyle+\;\left\\{g^{\prime\prime}(\theta^{(t)})+\frac{b_{3}}{2}[h(\theta^{(t)})-\theta^{(t)}]\right\\}\cdot\left\\{3[h^{\prime}(\theta^{(t)})-1]h^{\prime\prime}(\theta^{(t)})+[h(\theta^{(t)})-\theta^{(t)}]h^{\prime\prime\prime}(\theta^{(t)})\right\\}.$

By letting $t\to\infty$ and using $h^{\prime}(\theta^{*})=0$, (A.12) and
$h^{\prime\prime}(\theta^{*})=0$, we obtain

$|h^{\prime\prime\prime}(\theta^{*})|=\bigg{|}\frac{b_{3}-g^{\prime\prime\prime}(\theta^{*})}{g^{\prime}(\theta^{*})}\bigg{|},$

indicating that the US iteration (5.5) has a cubic convergence rate given by
(5.6). $\Box$

## Appendix B: Mode of the skew normal distribution

Let $x_{\rm mod}$ denote the mode of the skew normal distribution with pdf
given by (4.5), in this appendix, we will propose an MM algorithm to
iteratively calculate

$x_{\rm
mod}=\arg\;\max_{x\in{\mathbb{R}}}\log[f(x|\mu,\sigma^{2},\alpha)]=\arg\;\max_{x\in{\mathbb{R}}}\left\\{-\frac{(x-\mu)^{2}}{2\sigma^{2}}+\log\left[\Phi\left(\frac{x-\mu}{\sigma_{*}}\right)\right]\right\\},$

where $\sigma_{*}\triangleq\sigma/\alpha$. Let
${\boldsymbol{\theta}}=(\mu,\sigma_{*}^{2})^{\\!\top\\!}$ and define the left-
truncated normal density as

$h(x|x^{(s)},{\boldsymbol{\theta}})\triangleq\phi\left(x+x^{(s)}-\mu|{\boldsymbol{\theta}}\right)\cdot\Phi^{-1}\left(\frac{x^{(s)}-\mu}{\sigma_{*}}\right)\cdot
I(x<\mu),$

where $\phi(\cdot|{\boldsymbol{\theta}})$ denotes the pdf of
$N(\mu,\sigma_{*}^{2})$. By applying the integral version of Jensen’s
inequality, we have

$\displaystyle\log\left[\Phi\left(\frac{x-\mu}{\sigma_{*}}\right)\right]$
$\displaystyle=$
$\displaystyle\log\left[\int_{-\infty}^{x}\phi(y|{\boldsymbol{\theta}})\,\mbox{d}y\right]=\log\left[\int_{-\infty}^{\mu}\phi(z+x-\mu|{\boldsymbol{\theta}})\,\mbox{d}z\right]$
$\displaystyle=$
$\displaystyle\log\left[\int_{-\infty}^{\mu}\frac{\phi(z+x-\mu|{\boldsymbol{\theta}})}{h(z|x^{(s)},{\boldsymbol{\theta}})}\cdot
h(z|x^{(s)},{\boldsymbol{\theta}})\,\mbox{d}z\right]$ $\displaystyle\geqslant$
$\displaystyle\int_{-\infty}^{\mu}\log\left[\frac{\phi(z+x-\mu|{\boldsymbol{\theta}})}{h(z|x^{(s)},{\boldsymbol{\theta}})}\right]\cdot
h(z|x^{(s)},{\boldsymbol{\theta}})\,\mbox{d}z,$

where the equality holds iff $x=x^{(s)}$. Thus, the minorizing function can be
constructed as

$\displaystyle Q(x|x^{(s)})$ $\displaystyle\triangleq$
$\displaystyle-\frac{(x-\mu)^{2}}{2\sigma^{2}}+\int_{-\infty}^{\mu}\log\left[\frac{\phi(z+x-\mu|{\boldsymbol{\theta}})}{h(z|x^{(s)},{\boldsymbol{\theta}})}\right]\cdot
h(z|x^{(s)},{\boldsymbol{\theta}})\,\mbox{d}z$ $\displaystyle=$
$\displaystyle-\;\frac{(x-\mu)^{2}}{2\sigma^{2}}-\int_{-\infty}^{\mu}\frac{(z+x-2\mu)^{2}}{2\sigma_{*}^{2}}\cdot
h(z|x^{(s)},{\boldsymbol{\theta}})\,\mbox{d}z+\mbox{constant}$
$\displaystyle=$
$\displaystyle-\;\left(\frac{1}{2\sigma^{2}}+\frac{1}{2\sigma_{*}^{2}}\right)x^{2}+\left[\frac{\mu}{\sigma^{2}}+\frac{2\mu}{\sigma_{*}^{2}}-\frac{1}{\sigma_{*}^{2}}\int_{-\infty}^{\mu}z\cdot
h(z|x^{(s)},{\boldsymbol{\theta}})\,\mbox{d}z\right]x+\mbox{constant}.$

Let $\,\mbox{d}Q(x|x^{(s)})/\,\mbox{d}x=0$, we obtain the $(s+1)$-th MM
iteration:

$\displaystyle
x^{(s+1)}=\frac{1}{1+\alpha^{2}}\left[\mu(1+2\alpha^{2})-\alpha^{2}\int_{-\infty}^{\mu}z\cdot
h(z|x^{(s)},{\boldsymbol{\theta}})\,\mbox{d}z\right].$ (B.1)

## Appendix C: Inequalities on Riemann’s zeta function

First, for any $\theta>1$, we have

$\displaystyle\left(\theta-\frac{3}{2}\right)^{2}+\frac{3}{4}>0\Leftrightarrow\theta^{2}-3\theta+3>0\Leftrightarrow\frac{1}{\theta-1}<\frac{1}{(\theta-1)^{2}}+1.$
(C.1)

Next, we prove the following key inequality:

$\displaystyle\int_{1}^{\infty}\frac{\log(z+1)}{z^{\theta}}\,\mbox{d}z$
$\displaystyle=$
$\displaystyle\frac{1}{1-\theta}\int_{1}^{\infty}\log(z+1)\,\mbox{d}z^{1-\theta}$
(C.2) $\displaystyle=$
$\displaystyle\frac{\log(2)}{\theta-1}+\frac{1}{\theta-1}\int_{1}^{\infty}\frac{z^{1-\theta}}{z+1}\,\mbox{d}z$
$\displaystyle<$
$\displaystyle\frac{\log(2)}{\theta-1}+\frac{1}{\theta-1}\int_{1}^{\infty}\frac{z^{1-\theta}}{z}\,\mbox{d}z$
$\displaystyle=$
$\displaystyle\frac{\log(2)}{\theta-1}+\frac{1}{(\theta-1)^{2}}\;\stackrel{{\scriptstyle\rm(\ref{eqnfiveC.1})}}{{<}}\;\frac{\log(2)+1}{(\theta-1)^{2}}+\log(2).$

Finally, we obtain

$\displaystyle Z(\theta)$ $\displaystyle=$
$\displaystyle\sum_{x=1}^{\infty}\frac{1}{x^{\theta}}\geqslant\int_{1}^{\infty}\frac{1}{x^{\theta}}\,\mbox{d}x=\frac{1}{\theta-1},$
(C.3) $\displaystyle Z^{\prime}(\theta)$ $\displaystyle=$
$\displaystyle-\sum_{x=2}^{\infty}\frac{\log(x)}{x^{\theta}}\geqslant-\int_{1}^{\infty}\frac{\log(z+1)}{z^{\theta}}\,\mbox{d}z\;\stackrel{{\scriptstyle\rm(\ref{eqnfiveC.2})}}{{>}}\;-\frac{\log(2)+1}{(\theta-1)^{2}}-\log(2),$
$\displaystyle Z^{\prime\prime}(\theta)$ $\displaystyle=$
$\displaystyle\sum_{x=2}^{\infty}\frac{\log^{2}(x)}{x^{\theta}}\leqslant\int_{1}^{\infty}\frac{\log^{2}(z+1)}{z^{\theta}}\,\mbox{d}z=\int_{1}^{\infty}\frac{\log^{2}(z+1)}{1-\theta}\,\mbox{d}z^{1-\theta}$
(C.4) $\displaystyle=$
$\displaystyle\frac{\log^{2}(2)}{\theta-1}+\frac{2}{\theta-1}\int_{1}^{\infty}\frac{z}{z+1}\cdot
z^{-\theta}\log(z+1)\,\mbox{d}z$ $\displaystyle\leqslant$
$\displaystyle\frac{\log^{2}(2)}{\theta-1}+\frac{2}{\theta-1}\int_{1}^{\infty}\frac{\log(z+1)}{z^{\theta}}\,\mbox{d}z$
$\displaystyle\stackrel{{\scriptstyle\rm(\ref{eqnfiveC.2})}}{{<}}$
$\displaystyle\frac{\log^{2}(2)}{\theta-1}+\frac{2}{\theta-1}\left[\frac{\log(2)+1}{(\theta-1)^{2}}+\log(2)\right].$

## Appendix D: R program for solving the root
of high–order polynomial equation

    root.HPE = function(a3, a2, a1, a0, a.max = 10, m, x0 = 0)
    { # Function name: root.HPE(a3, a2, a1, a0, a.max = 10, m, x0 = 0)
      # ---------------------------Aim------------------------------------------
      # Solve the equation: 0 = a3*x^m + a2*x^2 + a1*x + a0, 0 < x < a.max
      # ---------------------------Input----------------------------------------
      #    a3: The coefficient of x^m
      #    a2: The coefficient of x^2
      #    a1: The coefficient of x
      #    a0: The constant, that is a positive real number
      # a.max: The maximum of x
      #     m: The highest order of the equation
      #    x0: An initial value and the default value is 0
      # ---------------------------Output---------------------------------------
      # The unique root of the polynomial equation in (0, a.max)
      ##########################################################################
      xt.save = c(); yt.save = c();
      xt = x0; y1t = a3*xt^m + a2*xt^2 + a1*xt + a0
      xt.save[1] = xt; yt.save[1] = y1t
      k = 1; esp = 1
      while(esp>1e-8){
        if(a3<0){                ### When a3<0, we apply the SLUF method
          if(y1t>0){ b2 = a3*m*(m-1)*a.max^(m-2) }else{b2 = 0}
          A2 = b2/2+a2
          A1 = a3*m*xt^(m-1)-b2*xt + a1
          A0 = a3*(1-m)*xt^m+b2/2*xt^2 + a0
        }
        else{                  ### When a3>0, we apply the TLF method
          A2 = a3*m*(m-1)/2*xt^(m-2)+a2
          A1 = a3*m*(2-m)*xt^(m-1)+a1
          A0 = a3*(m-1)*(m/2-1)*xt^m+a0
        }
        Delta.det = A1^2-4*A2*A0
        if(Delta.det >= 0){
          xt.save[k+1] = -(A1+sqrt(Delta.det))/(2*A2)
          yt.save[k+1] = a3*xt.save[k+1]^m + a2*xt.save[k+1]^2 +
            a1*xt.save[k+1]+a0
          esp = abs(yt.save[k+1])
          xt = xt.save[k+1]
          y1t = yt.save[k+1]
          k = k+1
        }
        else{
          stop(paste(’There exists no root in (0,’,a.max,’)’))
        }
      }
      return(xt)
    }
    # --------------------------------------------------------------------------
    # --------Example: x^3-3x^2+x+1 = 0-----------------------------------------
    > a3 = 1; a2 = -3; a1 = 1; a0 = 1
    > root.HPE(a3, a2, a1, a0, a.max=2, m=3, x0=0)
            [,1]      [,2]       [,3]         [,4] [,5]
    xt.save    0 0.7675919 0.99418291 9.999999e-01    1
    yt.save    1 0.4522631 0.01163398 1.968315e-07    0