# Bayesian Estimation of the Degrees of Freedom Parameter of the Student-$t$
Distribution—A Beneficial Re-parameterization

Darjus Hosszejni

###### Abstract

In this paper, conditional data augmentation (DA) is investigated for the
degrees of freedom parameter $\nu$ of a Student-$t$ distribution. Based on a
restricted version of the expected augmented Fisher information, it is
conjectured that the ancillarity DA is progressively more efficient for MCMC
estimation than the sufficiency DA as $\nu$ increases; with the break even
point lying at as low as $\nu\approx 4$. The claim is examined further and
generalized through a large simulation study and a application to U.S.
macroeconomic time series. Finally, the ancillarity-sufficiency interweaving
strategy is empirically shown to combine the benefits of both DAs. The
proposed algorithm may set a new standard for estimating $\nu$ as part of any
model.

##  1 Introduction

The Student-$t$ distribution is a standard item in the statistician’s toolbox.
Utilized as the regression error distribution, it robustifies parameter
estimation against extreme observations in the data set [9]. Furthermore, the
properties of Student-$t$ are classical building blocks of hierarchical models
for financial Bollerslev1987,Harvey1994 and polychotomous data [1], among
others.

Maximum likelihood (ML) estimation of the degrees of freedom ($\nu$) parameter
can be cumbersome due to the complexity of the likelihood function; [23] give
an overview and develop a quasi-Newton, an expectation-maximization (EM), and
a score function-based approach to ML. Liu1997 reviews extensions to the EM
algorithm for the multivariate Student-$t$ distribution. In the non-asymptotic
branch of frequentist literature, [8] discusses structural inference for the
parameters, and Sutradhar1986,Singh1988 provide method of moments estimators
for $\nu$.

Bayesian treatments of robust regression date back to [7], who explores scale-
mixtures of normal distributions as a tool for handling outliers and presents
the Cauchy distribution as an example, which corresponds to $\nu=1$. Still
with a fixed $\nu$, [44] discusses posterior inference for the parameters of a
multivariate linear Student-$t$ regression using conjugate prior
distributions. Ramsay1980,West1984 further extend de Finetti’s framework and
provide more examples including the Student-$t$ distribution. Chib1991
demonstrate conditions, when this type of robust regression can not work due
to a non-informative likelihood for $\nu$.

In terms of Bayesian inference for $\nu$ using Markov chain Monte Carlo
methods [[, MCMC;]]Tierney1994, the Gibbs sampler [12] by [15] has been
established as the standard method and it comprises the main focus in this
paper. The Gibbs sampler relies on the scale-mixture representation of the
Student-$t$ distribution to draw a posterior sample using a model with data
augmentation [[, DA;]]VanDyk2001. There exist many options for DA; however, to
our knowledge, only the one outlined by [15], termed sufficient augmentation
(SA), has been in use for Student-$t$.

Multiple examples demonstrate that SA may not work well for sampling $\nu$;
see, for instance, [14, efficiency values in Appendix C], or [26], who blame
the hierarchical structure and the lack of marginalization for the low
sampling efficiency. However, inspired by the general framework of [30], an
alternative way to higher efficiency is found. They demonstrate through
multiple example models that the transformation of SA to the ancillary
augmentation (AA) often boosts the efficiency already quite drastically;
unfortunately, whether AA or SA is superior depends on the data set. [29]
propose as the way out of this data-dependence a hybrid Gibbs-sampler, which
randomly picks AA or SA in every pass; they conclude that a combination of
sufficiency and ancillarity may be desirable in any setting. Finally, [43]
introduce the ancillarity-sufficiency interweaving strategy (ASIS), which
exploits the complementary roles of AA and SA in any hierarchical model and
integrates them into one generally effective recipe. The applicability and
utility of ASIS is documented in the literature [[, see,
e.g.,]]Kastner2014,Bitto2019,Monterrubio-Gomez2020.

The contribution of this paper is two-fold: the ancillary parameterization for
the Student-$t$ distribution is introduced, and it is shown through
simulations that the AA sampler explores the posterior space more efficiently
than the usual SA sampler for a large set of data generating processes;
furthermore, an implementation of ASIS is proposed that can be applied as a
plug-and-play substitute in any MCMC procedure for models that include the
Student-$t$ distribution with unknown degrees of freedom.

The remaining part of the paper is structured as follows: Section 2 introduces
the model and the notation. Section 3 provides details on measurement of
efficiency, the Gibbs sampling algorithm, and the corresponding conditional
distributions. Section 4 showcases the efficiency of the introduced method on
a large grid of simulated data sets, and Section 5 is an application inspired
by [16]. Section 6 concludes.

##  2 Simple Student-$t$-Model

The simplest setup with Student-$t$ errors is investigated;

$y\sim t_{\nu}(0,1),$ (1)

where $t_{\nu}(\mu,\sigma^{2})$ denotes the Student-$t$ distribution with
$\nu$ degrees of freedom, location $\mu$, and variance $\sigma^{2}$; the
kernel of its density function is $(1+x^{2}/\nu)^{-(\nu+1)/2}$. Taking this
simple model enables a pure analysis of sampling $\nu$, decoupled from other
effects present in larger models; furthermore, Equation (1) is an ingredient
of all models employing the Student-$t$ distribution. At the same time, the
simplicity of Equation (1) is not restrictive either due to the modular nature
of Gibbs-sampling: in case an econometric model is extended by Student-$t$
errors, its corresponding MCMC procedure shall merely include the additional
steps for drawing $\nu$. This modularity is demonstrated in Section 5 through
a more complex model.

As [44] points out, the density function of the Student-$t$ distribution can
be formulated as a scale mixture of normal distributions, which inspires the
introduction of missing data $\tau$ and the data augmented representation

$\begin{split}y\sim&\;\mathcal{N}({0},{\tau}),\\\
\tau^{-1}\sim&\;\mathcal{G}({\nu/2},{\nu/2}),\end{split}$ (2)

where $\mathcal{N}({\mu},{\sigma^{2}})$ denotes the univariate normal
distribution with mean $\mu$ and variance $\sigma^{2}$, and
$\mathcal{G}({\alpha},{\beta})$ is the gamma distribution with shape $\alpha$
and rate $\beta$. [30] call Equation (2) the centered parameterization, while,
following [43], we prefer the term sufficient parameterization. This
terminology derives from the fact that $\tau$ is a sufficient statistic for
$\nu$.

Based on Equation (2), it is possible to derive other data augmentation
schemes, where the auxiliary variable and $\nu$ are a priori independent;
these are called ancillary parameterization. As [30] point out, one such
option is to replace $\tau$ by $u=F(\tau;\nu)$, where $F(x;\nu)$ is the
cumulative distribution function (CDF) for the prior distribution of $\tau$.
Note that $u$ follows a standard uniform distribution a priori; therefore, the
ancillary data augmentation for Equation (1) is

$\begin{split}y\sim&\;\mathcal{N}({0},{F^{-1}(u;\nu)}),\\\
u\sim&\;\mathcal{U}(0,1),\end{split}$ (3)

where $F^{-1}(u;\nu)$ is the inverse CDF for the prior distribution of $\tau$.

As demonstrated by [43], for improved sampling efficiency, it is beneficial to
consider both a sufficient and an ancillary data augmentation scheme in the
design of an MCMC sampling algorithm. The added value in using both stems in
the differences between the joint distributions $p(y,\nu,\tau)$ and
$p(y,\nu,u)$; two slices for each distribution are depicted in Figures 1 and
2, respectively. The figures show contour plots of unnormalized111There is no
closed form for these distributions, and normalization is done over the shown
window, which contains most of the probability mass. posterior densities
$p(\nu,\tau\mid y=y_{0})$ and $p(\nu,u\mid y=y_{0})$ for small ($y_{0}=0$) and
large ($y_{0}=4$) observation error. The top panel of Figure 1 depicts
rotating shapes, which signals that the dependence between $\nu$ and $\tau$ is
different at different parts of the space. In comparison, this rotation is
non-visible in the bottom panel, and, thus, dependence between $\nu$ and $u$
seems less complex; moreover, for the shown range, the shapes of isolines
resemble concentric standard ellipses, which is a sign for approximate
posterior independence between $\nu$ and $u$; a Gibbs-sampler can explore such
a space very efficiently. Figure 2 shows an example, where a suspect outlier
is modeled using Student-$t$ errors and a small value for $\nu$ is expected;
here, the dependence between $\nu$ and $\tau$ varies less in space than the
dependence between $\nu$ and $u$. In conclusion, SA seems beneficial for fat
tailed observations; on the other hand, AA may excel for lighter tailed
observations.

Figure 1: Contour plot of the bivariate density functions $p(\nu,\tau\mid
y=0)$ (top) and $p(\nu,u\mid y=0)$ (bottom), where $\nu\in[1,30]$,
$\tau\in(0,20]$, and $u\in(0,1)$.

Figure 2: Contour plot of the bivariate density functions $p(\nu,\tau\mid
y=4)$ (top) and $p(\nu,u\mid y=4)$ (bottom), where $\nu\in[1,30]$,
$\tau\in(0,4]$, and $u\in(0,1)$.

The Bayesian description of the model is completed by the prior distribution
$p(\nu)$. Following [15], an exponential distribution is assumed;

$\nu\sim\mathcal{E}(\lambda),$ (4)

where $\lambda$ denotes the rate parameter. This choice is motivated by the
availability of a suitable sampling procedure for SA using this prior [14]. It
is to be noted, that the adaptation of the proposed methods to alternative
priors, such as the uniform distribution over an interval $[a,b]$ [5] or a
gamma distribution [27], is uncomplicated. This is discussed in Section 3.2.

##  3 Bayesian Estimation

Two Gibbs-samplers are implemented that correspond to AA and SA, respectively;
both of them consist of sequentially drawing random values from two
conditional posterior distributions. Above all, primary interest lies in how
quickly the two resulting Markov chains mix, given different sets of
observations; i.e., how close the posterior sample for $\nu$ is to a set of
independent draws. In particular, fast mixing is a desirable requirement for
any MCMC procedure. [40] discuss several measures for mixing and propose to
minimize the expected augmented Fisher information matrix

$\begin{split}I_{\ast}(y,\nu)=\text{E}\left[-\frac{\partial^{2}\log p(\nu\mid
y,\ast)}{\partial\nu\cdot\partial\nu}\;\middle|\;y,\nu\right]\end{split}$

for choosing between DAs, where $\ast\equiv\tau$ for SA and $\ast\equiv u$ for
AA. Although a simple derivation gives
$I_{\tau}(y,\nu)=n(\Psi(\nu/2)-1/\nu)/2$, where $n$ denotes the number of
observations and $\Psi(z)=\text{d}\log(\Gamma(x))/\text{d}z$ is the digamma
function, we have been unsuccessful at computing $I_{u}$ analytically.
Fortunately, Monte Carlo (MC) estimation is possible. Figure 3 exemplifies the
difference between the estimated $I_{u}$ and $I_{\tau}$ for data sets of size
1. The break even point (BEP) between AA and SA is at $\nu\approx 4$ almost
independently from $y$; that is where AA overtakes SA at this benchmark as
$\nu$ increases.222For the evaluation of $I_{u}(y_{0},\nu_{0})$, a MC sample
$\vec{u}_{\text{MC}}=(u^{(1)},\dots,u^{(L)})^{\top}$ of length $L=10000$ was
simulated from $p(u\mid y=y_{0},\nu=\nu_{0})$ using the algorithm in Section
3.2; next, $\partial^{2}\log p(\nu\mid
y=y_{0},u=u^{(i)})/\partial\nu\cdot\partial\nu$ is evaluated at points
$i=1,\dots,L$; finally, the average of the values is taken to be the estimator
for $I_{u}(y_{0},\nu_{0})$. The numerical second derivative of $\log p(\nu\mid
y,u)$ is approximated using Richardson’s extrapolation with 6 steps, which is
implemented in the hessian function from the software package numDeriv [18].

Figure 3: Contour plot of the estimated difference
$I_{u}(y,\nu)-I_{\tau}(y,\nu)$ in the region
$(y,\nu)^{\top}\in[0,10]\times[0.1,100]$. The depicted function is symmetric
about $y=0$. Note that the horizontal axis is on a logarithmic scale.

The comparison of $I_{\tau}$ and $I_{u}$ through MC quickly becomes
impractical for larger data sets. Therefore, an extensive computer simulation
is conducted and empirical measurement of mixing is computed. Ideally, the
latter is based on exact sampling from the four conditional posterior
distributions.

###  3.1 Conditional Posterior Distributions

Given data points $\vec{y}=(y_{1},\dots,y_{n})^{\top}$, the specification for
$\vec{\tau}=(\tau_{1},\dots,\tau_{n})^{\top}$ is conditionally conjugate [15],
and, therefore, sampling from $p(\vec{\tau}\mid\vec{y},\nu)$ is
straightforward. Next, as [30] mention, sampling from
$p(\vec{u}\mid\vec{y},\nu)$ can be reduced to sampling $\vec{\tau}$.

The rejection sampling (RS) method by [15] is employed for
$p(\nu\mid\vec{y},\vec{\tau})$: denote by $\xi^{\ast}$ the unique solution to

$(\log(\xi/2)+1-\Psi(\xi/2))n/2+\xi^{-1}-\eta=0,$

where $\eta=0.5\sum_{i=1}^{n}(\log(\tau_{i})+1/\tau_{i})-\lambda$ is a
sufficient statistic. By finding $\xi^{\ast}$, the normalizing constant in
$p(\nu\mid\vec{y},\vec{\tau})$ is estimated indirectly, which is necessary for
RS. Then, a proposal $\nu_{\text{p}}$ is drawn from an exponential
distribution with rate $1/\xi^{\ast}$, and it is retained with probability

$\exp(\eta\xi^{\ast}-1)\left(\frac{(\xi^{\ast}/2)^{\xi^{\ast}/2}}{\Gamma(\xi^{\ast}/2)}\right)^{-n}\exp(\nu_{\text{p}}(1/\xi^{\ast}-\eta))\left(\frac{(\nu_{\text{p}}/2)^{\nu_{\text{p}}/2}}{\Gamma(\nu_{\text{p}}/2)}\right)^{n}.$

Last, we are unaware of exact sampling procedures for
$p(\nu\mid\vec{y},\vec{u})$. A modern and easy-to-implement substitute is to
employ adaptive Metropolis [[, AM;]]Haario2001 within Gibbs, which can be used
to sample from any unnormalized density that can be evaluated. However,
contrary to RS, AM produces an auto-correlated Markov chain, and therefore it
is inferior to RS in our setup. As a simple improvement, independent sampling
is approximated by repeating AM in one step of the Gibbs-sampler. Our
implementation of AM proposes $\log(\nu_{\text{p}})$ for the logarithm of
$\nu$ through a Gaussian random walk with variance $\sigma^{2}$. Next,
$\nu_{\text{p}}$ is accepted or rejected according to the Metropolis ratio
$\alpha$. The proposal distribution is tuned to achieve the famous
$\alpha=0.44$ [11, 10, 34] by increasing or decreasing $\sigma$ after batches
of several hundred draws; this is similar to the examples by [35].

###  3.2 MCMC Algorithm

Following [15], Equation (2) naturally translates into the following MCMC
scheme, henceforth termed algorithm SA, which simulates from the posterior
distribution $p(\nu,\vec{\tau}\mid\vec{y})$:

  1. 1.

Initialize $\nu$;

  2. 2.

Draw
$\tau_{i}^{-1}\sim\mathcal{G}\left({\frac{\nu+1}{2}},{\frac{\nu+y_{i}^{2}}{2}}\right)$
independently for each $i=1,\dots,n$;

  3. 3.

Draw $\nu\sim p(\nu\mid\vec{y},\vec{\tau})$, where
$\vec{\tau}=(\tau_{1},\dots,\tau_{n})^{\top}$ using RS;

  4. 4.

Repeat Steps 2 and 3.

For alternative prior distributions, Step 3 can be replaced by AM as long as
it is possible to evaluate the prior density function $p(\nu)$. This change
may affect the mixing speed for two reasons: autocorrelation is introduced
into the Markov chain, and $I_{\tau}$ might be altered.

Next, the two novel algorithms are introduced. Corresponding to Equation (3),
algorithm AA samples from the posterior distribution $p(\nu,u\mid\vec{y})$;

  1. 1.

Initialize $\nu$;

  2. 2.

Draw $u_{i}\sim p(u_{i}\mid y_{i},\nu)$ independently for each $i=1,\dots,n$
by first drawing
$\tau_{i}^{-1}\sim\mathcal{G}\left({\frac{\nu+1}{2}},{\frac{\nu+y_{i}^{2}}{2}}\right)$
and then computing $u_{i}=F(\tau_{i};\nu)$ [[, inspired
by]]Papaspiliopoulos2007;

  3. 3.

Draw $\nu\sim p(\nu\mid\vec{y},\vec{u})$, where
$\vec{u}=(u_{1},\dots,u_{n})^{\top}$ using AM; this step includes the
evaluation of $F^{-1}(u_{i};\nu)$; repeat this step $k_{\text{AA}}$ times;

  4. 4.

Repeat Steps 2 and 3.

This is algorithm is trivially modified for other prior distributions: the
Metropolis ratio $\alpha$ needs to be adapted to the new prior density.

Finally, algorithm ASIS combines SA and AA:

  1. 1.

Initialize $\nu$;

  2. 2.

Do Steps 2 and 3 of SA;

  3. 3.

Compute $u_{i}=F(\tau_{i};\nu)$ for $i=1,\dots,n$;

  4. 4.

Do Step 3 of AA;

  5. 5.

Repeat Steps 2 through 4.

Notice that ASIS is not a sequence of SA and AA; $\vec{u}$ is not sampled
separately. Rather, focus lies on the differences between
$p(\nu\mid\vec{y},\vec{\tau})$ and $p(\nu\mid\vec{y},\vec{u})$, which are
exploited in a minimal setup.

Correctness of the implementations for AA, SA, and ASIS is confirmed using a
variant on Geweke’s test [17]. Specifically, if a new vector of observations
is simulated at the end of every round of the particular algorithm, then the
sampling distribution of $\nu$ is equivalent to its prior; and this is
verified visually with the help of quantile-quantile-plots.

##  4 Simulation Study

In order to make a general comparison between the mixing of AA, SA, and ASIS,
an extensive simulation study is conducted. To this end, the algorithms in
Section 3.2 are implemented in the R language [32, version 4.1.0]. Then,
multiple independent posterior Markov chains are generated for data simulated
from Equation (1); for AA, $k_{\text{AA}}=20$ is set. Next, convergence of the
Markov chains is automatically validated by ensuring that their $\hat{R}$
value [41] is smaller than the recommended $1.1$; this is done using the rhat
function of the R package posterior rposterior. Finally, relative numerical
efficiencies [[, RNE;]]Geweke1989 of the chains are reported. RNE is defined
for an identically distributed but serially correlated Monte Carlo sample
$\\{\theta^{(j)}\\}$ and a function of interest $g$ as
$\hat{\text{var}}(g)/\hat{S}(0)$, where $\hat{\text{var}}(g)$ is a consistent
estimator for the variance of the distribution $p(\theta)$ and
$\hat{S}(\omega)$ is a consistent estimator for the spectral density
$S(\omega)$ of $\\{g(\theta^{(j)})\\}$, $j=1,\dots,M$.333RNE is the reciprocal
of the inefficiency factor reported in other works Nakajima2009a. This
computation is implemented in function effectiveSize from the coda package
[31].

###  4.1 Setup of Data Generating Process

A grid of true values $\nu_{\text{true}}=1,1.5,2,2.5,3,4,5,10,20,50,100$ and
data set lengths $n=1,3,10,30,100,300,1000,3000,10000$ is considered for the
comparison. This covers “realistic” values [[,
$\nu_{\text{true}}>2$;]]Gelman2013a and both problematic regions of low and
high degrees of freedom. Moreover, dependence on the size of the data set can
be evaluated. Then, for each of the 99 setups, five data sets are simulated
from Equation (1). Furthermore, to examine prior sensitivity, values
$\lambda=0.05,0.1,0.2,0.5,1$ are considered, which correspond to prior
expectations of $1$ through $20$ for $\nu$. In addition, following the
recommendations for the computation of $\hat{R}$, four independent Markov
chains are simulated with initial values $\nu=0.5,2,10,100$ (these receive the
same $\hat{R}$ during post-processing), for each of the three algorithms and
each setup. Overall, $29700$ samples of length $M=10000$ are drawn after a
burn-in phase of length $1000$, and samples with $\hat{R}>1.1$ are removed.

###  4.2 Posterior Estimates

A sample of interval estimators for $\nu$ is presented in Table 1 for $n=1000$
and $\lambda=0.2$. With the exception of two cells, the estimates by AA, SA,
and ASIS are similar for all values of $\nu_{\text{true}}$. The twin of these
exceptions is an example of a general caveat about AA: for a long vector of
heavy-tailed observations, the numerical evaluation of $F^{-1}$ may fail
within AA; especially for values $\nu$ at the outskirts of the posterior
distribution. The issue materializes as a constant estimator $(100,100)$ in
Table 1.444 Our implementation for $F^{-1}$ is based on the built-in R
function qgamma. Typically, the Markov chain is constant in these settings
because the MH step always rejects; in all other settings, AA is reliable.
Consequently, by choosing (our implementation of) AA, one assumes a complex
restriction on the prior distribution $p(\nu,\vec{u})$. Moreover, the three
rows for $\nu_{\text{true}}=100$ provide two further insights. First,
obviously, even at sample size $n=1000$, most of the posterior mass is closer
to the prior expectation of $1/\lambda=5$ than to 100; indeed, the prior
choice is quite influential compared to the flat likelihood for light-tailed
data. Second, the variation in the estimates from SA is the largest in this
row both in relative and in absolute terms; the $90\%$ quantile ranges between
$37.9$ and $40.7$. This sheds light on the outstandingly deteriorating
effectiveness of SA as $\nu_{\text{true}}$ increases; this is discussed in
Section 4.3. As a result, its finite sample properties make SA less reliable
for estimating a large $\nu$. Importantly, ASIS does not suffer from any of
the aforementioned issues; indeed, it combines the best of both worlds.

Table 1: Estimated 10th and 90th percentiles for $\nu$. The same data of length $n=1000$ is used inside each group of $\nu_{\text{true}}$; further, $\lambda=0.2$ for each entry. Sets of chains with $\hat{R}\geq 1.1$ are marked. $\nu_{\text{true}}$ |  | Initial $\nu=0.5$ | 2 | 10 | 100
---|---|---|---|---|---
| AA | $(0.894,1.01)^{\ast}$ | $(0.895,1.01)^{\ast}$ | $(0.893,1.01)^{\ast}$ | $(100,100)^{\ast}$
| SA | $(0.893,1.01)$ | $(0.894,1.01)$ | $(0.893,1.01)$ | $(0.895,1.01)$
1 | ASIS | $(0.894,1.01)$ | $(0.893,1.01)$ | $(0.894,1.01)$ | $(0.894,1.01)$
| AA | $(1.74,2.05)^{\ast}$ | $(1.75,2.05)^{\ast}$ | $(1.75,2.05)^{\ast}$ | $(100,100)^{\ast}$
| SA | $(1.75,2.05)$ | $(1.74,2.05)$ | $(1.75,2.05)$ | $(1.74,2.04)$
2 | ASIS | $(1.75,2.05)$ | $(1.75,2.05)$ | $(1.75,2.05)$ | $(1.74,2.04)$
| AA | $(3.99,5.25)$ | $(3.98,5.23)$ | $(3.99,5.26)$ | $(3.98,5.23)$
| SA | $(4,5.31)$ | $(3.94,5.21)$ | $(3.98,5.25)$ | $(3.99,5.26)$
5 | ASIS | $(3.98,5.23)$ | $(3.98,5.24)$ | $(3.98,5.24)$ | $(3.99,5.24)$
| AA | $(10,19.4)$ | $(10.2,19.1)$ | $(10.2,19.3)$ | $(10.2,19.4)$
| SA | $(10.2,19.4)$ | $(10.2,19.4)$ | $(10.4,19.8)$ | $(10.2,18.6)$
20 | ASIS | $(10.2,19.3)$ | $(10.1,19.1)$ | $(10.1,19.3)$ | $(10.2,19.1)$
| AA | $(16.6,39.9)$ | $(16.5,40.1)$ | $(16.7,39.9)$ | $(16.7,40.2)$
| SA | $(16.7,40.7)$ | $(16.8,37.9)$ | $(17.1,38.4)$ | $(17,38.9)$
100 | ASIS | $(16.6,39.9)$ | $(16.7,40.1)$ | $(16.7,39.8)$ | $(16.7,40.2)$

###  4.3 Sampling Efficiency

Detailed results are presented in Table 2. In most cases of practical
interest, i.e., where the observations have finite variance by assumption, the
RNE of AA is at least comparable to that of the classical SA; and, for
lighter-tailed inputs, AA is incomparably more efficient. Nevertheless, as
discussed in Section 4.2, AA may fail; these cases are marked as missing
values in Table 2. In particular, AA breaks down for the combination of a data
set of length $n\gtrapprox 100$ simulated with $\nu_{\text{true}}\lessapprox
3$, and a chain with initialization $\nu=100$; however, apparent from the
results, lighter-tailed data sets may be increasingly affected as their size
grows. On the other hand, AA is reliable in all other settings. Luckily, ASIS
proves to be a solution in this regard, as SA does not suffer from the issue;
therefore, the chain can move out of this initial bad region. Moreover, ASIS
outperforms both AA and SA in all setups. Even an RNE of 100% is reached for
the small $n=10$, which is measurably on par with an i.i.d. sample.

Table 2: Mean RNE (%) for AA, SA, and ASIS, and different values of $n$ and $\nu_{\text{true}}$.  |  | $\nu_{\text{true}}=1$ | 1.5 | 2 | 2.5 | 3 | 4 | 5 | 10 | 20 | 50 | 100
---|---|---|---|---|---|---|---|---|---|---|---|---
| $n=10$ | 23.2 | 39.6 | 53.1 | 63.9 | 70.2 | 77.6 | 81.7 | 87.7 | 90.1 | 90.9 | 91.1
| 100 | 4.3 | 15.4 | 22.0 | 30.5 | 38.2 | 45.5 | 51.6 | 60.4 | 68.0 | 70.9 | 71.3
| 1000 | - | - | 23.3 | 29.9 | 36.3 | 43.7 | 48.1 | 45.1 | 39.3 | 44.4 | 46.6
AA | 10000 | - | - | - | - | 37.3 | 43.2 | 48.9 | 45.3 | 28.4 | 21.8 | 24.5
| 10 | 45.5 | 31.8 | 25.7 | 22.5 | 21.2 | 19.8 | 19.5 | 19.2 | 19.7 | 19.8 | 19.8
| 100 | 45.4 | 32.7 | 24.7 | 19.2 | 15.4 | 10.0 | 7.2 | 4.4 | 4.0 | 3.9 | 3.9
| 1000 | 44.4 | 31.3 | 23.2 | 17.8 | 14.2 | 9.5 | 6.8 | 2.2 | 1.0 | 0.9 | 0.9
SA | 10000 | 44.6 | 31.5 | 23.5 | 18.2 | 14.5 | 9.9 | 7.1 | 2.3 | 0.7 | 0.3 | 0.2
| 10 | 76.7 | 80.8 | 86.0 | 90.4 | 93.1 | 95.3 | 96.8 | 99.4 | 99.8 | 100.6 | 100.6
| 100 | 63.1 | 59.5 | 58.9 | 60.4 | 62.1 | 65.5 | 65.4 | 68.3 | 73.5 | 74.3 | 76.1
| 1000 | 63.0 | 60.3 | 60.6 | 61.6 | 63.1 | 64.1 | 64.1 | 49.6 | 40.9 | 45.3 | 48.0
ASIS | 10000 | 62.5 | 59.7 | 60.5 | 62.0 | 63.3 | 65.3 | 64.4 | 50.0 | 29.3 | 21.9 | 25.1

Some further results are highlighted using visualizations. First, Figure 4 is
a demonstration of the differences between AA, SA, and ASIS using trace plots:
these are shown for hand-picked but representative examples. The top three
panels correspond to very heavy-tailed input data ($\nu_{\text{true}}=1$) and
a good initialization for AA; SA and ASIS mix quickly, but AA seems to
slightly lag behind. In contrast, the bottom three panels represent the case
of a light-tailed input vector ($\nu_{\text{true}}=100$); now, AA and ASIS
show good performance, but SA is highly ineffective.

Figure 4: Example trace plots for a posterior sample from $p(\nu\mid\vec{y})$,
where $n=1000$, $\lambda=0.2$, and the data $\vec{y}$ is simulated using
$\nu=1$ (top) and $\nu=100$ (bottom). For better visibility, only $5000$ draws
are shown.

Next, Figure 5 depicts a sample of RNE values conditional on the algorithm,
$\nu_{\text{true}}$, $\lambda$, and the simulated data set; for improved
readability, the plot is restricted to $n=100$,
$\nu_{\text{true}}\in\\{1,2,3,5,10\\}$, and four out of five data sets.
Clearly, the sampling efficiency strongly depends on the exact properties of
the observations. For instance, the data set with $DataID=1$ is an outlier in
the panel for $\nu=5$; and the reason for the comparably poor performance of
AA is the extreme outlier $\simeq 15$ contained in the data. Similar behavior
can be observed in the rightmost panel for $DataID=1$, which is due to the
same random seed being used for generating data sets with the same $DataID$.
Further, at this sample size, $\lambda$ does not affect RNE. Indeed, there is
no clustering of colors in any of the 15 columns. Since $\lambda$ cancels out
in the Fisher information, this speaks for the similarity between RNE and
$I_{\tau}$ or $I_{u}$. In summary, part of the conclusion of Figure 5 is
consistent with that of Figure 3: SA outperforms AA only for particularly
heavy-tailed input data. However, the data do affect the BEP, which, in
addition, comes at a lower $\nu_{\text{true}}$ than the aforementioned value
4: between 2 and 3. Finally and most importantly, ASIS is immune to both the
numerical issues of AA and the quickly deteriorating efficiency of SA.

Figure 5: RNE in dependence on four variables: the algorithm (AA, SA, or
ASIS), $\nu$, $\lambda$, and the input data set; $n=100$ is fixed. Points
correspond to the individual Markov chains; therefore, four points belong to
each setting. For increased visibility, only four out of the five simulated
data sets are included, and horizontal jitter is applied to each point.
Algorithm AA has only one chain for $\nu=1$ with an $\hat{R}<1.1$.

The significance of ASIS becomes even more apparent through Figure 6. A
summary of RNE is displayed for AA, SA, and ASIS, as it decreases with
increasing $n$ and $\nu_{\text{true}}$. Above all, the larger
$\nu_{\text{true}}$ is, the greater the difference between the slopes for SA
and ASIS. Based on the bottom panel, one may hypothesize that the RNE of SA is
linearly proportional to $1/n$ for light-tailed input; this should be
connected to the multiplicative term $n$ in $I_{\tau}$. Finally, for the last
time, the dominance of ASIS is conspicuous.

Figure 6: RNE in dependence on $n$, $\nu_{\text{true}}$, and the algorithm.
Quantiles are computed over the $100$ Markov chains that span over the five
values for $\lambda$, four parallel chains, and five simulated data sets; as
described in Section 4.1. Both axes are on a logarithmic scale. Note that
algorithm AA does not have chains for $n\geq 1000$, whose $\hat{R}<1.1$.

##  5 Application

The efficiency of ASIS is evaluated on an extension of the 14 U.S.
macroeconomic time series by [28] as part of the model specified by
Geweke1992,Geweke1993,Geweke1994a; these three closely connected publications
are followed in this section. The data set can be downloaded with the R
package tseries rlang2021,tseries2020.555The data set was formerly stored in
the data archive of the Journal of Business and Economic Statistics, currently
at http://korora.econ.yale.edu/phillips/data/np&enp.dat. The observation
series are annual, they start between 1860 and 1909, and all end in 1988.

###  5.1 Model Specification

The model is a re-parameterized auto-regressive linear model of lag five with
a time trend and Student-$t$ increments, and, as highlighted by [15], it is an
alternative to the specifications employed by [28] and others. To underline
the time-series nature of the data set, index $t$ is used in this section
instead of $i$. Formally, observations $\\{y_{t}\\}$ for $t=1,\dots,T$ are
assumed to be explained through

$\begin{split}y_{t}&=\gamma+\delta t+u_{t},\\\
A^{\star}(L)u_{t}&=\varepsilon_{t},\\\ A^{\star}(L)&=(1-\rho
L)+(1-L)\sum_{j=1}^{4}a_{j}L^{j},\\\ \\{\varepsilon_{t}\\}&\sim\text{i.i.d.
}t_{\nu}(0,\sigma^{2}),\end{split}$ (5)

where $L$ is the lag operator, $\gamma\in\mathbb{R}$ is the intercept,
$\delta\in\mathbb{R}$ is the slope of the main time trend, and the series
$\\{u_{t}\\}$ stands for the cyclic deviations from the trend. The cycle is
assumed to vary according to the fifth degree polynomial specification for
$A^{\star}$ with parameters $0\leq\rho<1$ and
$a_{1},a_{2},a_{3},a_{4}\in\mathbb{R}$.

It is possible to eliminate the lag operator [15] and write down Equation (5)
as

$\displaystyle
y_{t}=\gamma(1-\rho)+\delta\left(\rho-\sum_{j=1}^{4}a_{j}\right)+\delta(1-\rho)t+\rho
y_{t-1}+\sum_{j=1}^{4}a_{j}(y_{t-j}-y_{t-j-1})+\varepsilon_{t}.$

Following [15], the following prior distributions are assumed;

$\begin{split}p(\rho)&=5\rho^{4}I_{[0,1]}(\rho),\\\
\delta&\sim\mathcal{N}({0},{0.05^{2}}),\\\
a_{j}&\sim\mathcal{N}({0},{0.731\cdot 0.342^{j-1}}),\\\
\gamma&\sim\mathcal{N}({y_{0}},{10^{2}}),\\\
p(\sigma)&\propto\sigma^{-1}I_{(0,\infty)}(\sigma),\end{split}$ (6)

where $y_{0}$ is the first observation and it is not included in (5), and
$I_{A}(x)$ is the indicator function that takes $1$ for $x\in A$ and $0$
otherwise. These prior distributions are in line with the hypotheses in the
original work, where the model is employed with Gaussian increments
($\nu=\infty$) to show that business cycles deviate from a linear time trend
in a non-stationary manner. Finally, Equation (4) is adopted for $p(\nu)$ with
$\lambda=0.333$.

A variant on the algorithm by [15] is implemented;

  1. 1.

Initialize $(\gamma,\delta,\rho,a_{1},a_{2},a_{3},a_{4})^{\top}$ as the
ordinary least squares (OLS) estimate (enforcing the bounds for $\rho$ if
needed), $\sigma^{2}$ as the sample variance given the OLS estimates,
$\tau\equiv 1$, and $\nu=4$;

  2. 2.

Draw $(\gamma,\delta)^{\top}$ given all other variables from a bivariate
normal distribution;

  3. 3.

Draw $(a_{1},a_{2},a_{3},a_{4})^{\top}$ given all other variables from a
multivariate normal distribution;

  4. 4.

Draw $\rho$ given all other variables from a non-standard distribution using
rejection sampling;

  5. 5.

Draw $\nu$ and $\vec{\tau}$ (for SA), $\vec{u}$ (for AA), or both (for ASIS)
given all other variables as in Section 3.2; note that $\\{\varepsilon_{t}\\}$
is a sufficient statistic;

  6. 6.

Draw $\sigma^{2}$ given all other variables from a scaled inverse-chi-square
distribution with $T$ degrees of freedom.

Our modification only affects Step 5, where the method described in Section 3
is applied. The exact distributions and methods of the remaining steps can be
found in Geweke1993 and in the software code as part of the Supplementary
Material.

Similarly to Section 3.2, correctness of our implementations for Steps 1-5 is
ensured through Geweke’s test [17]; the prior distributions are recovered by
simulating data from Equation (5) at the end of each pass. This does not work
for Step 6 due to $p(\sigma)$ being improper. There, the empirical quantile
function of the conditional posterior distribution is compared to the
theoretical quantile function.

###  5.2 Posterior Estimates

For all 14 variables, a posterior sample of size $10000$ is drawn after a
burnin of $1000$ passes using AA, SA, and ASIS, resulting in 42 Markov chains;
after the verification of convergence through trace plots, these numbers are
considered to be adequate. There were no numerical issues with the computation
of $F^{-1}$ as part of AA.

Estimated posterior distributions are summarized in Table 3. The three
algorithms largely agree on the results, and there do not seem to be
systematic biases present. The figures suggest highly heavy-tailed posterior
distributions for $\\{\varepsilon_{t}\\}$ in all cases, given the prior
information. Median values range between approximately $1.2$ and $6$, and,
therefore, they fall on both sides of the previously measured BEP’s.

Table 3: Summary of posterior distributions for $\nu$. Credible intervals show the 10th and the 90th posterior percentiles after 10000 draws. | Median | 80% Credible Interval
---|---|---
| AA | SA | ASIS | AA | SA | ASIS
CPI | $2.21$ | $2.15$ | $2.17$ | $(1.57,3.31)$ | $(1.55,3.27)$ | $(1.54,3.21)$
Employment | $2.3$ | $2.34$ | $2.29$ | $(1.47,4.03)$ | $(1.5,4)$ | $(1.49,4.09)$
GNP Deflator | $2.23$ | $2.21$ | $2.2$ | $(1.58,3.27)$ | $(1.54,3.25)$ | $(1.55,3.26)$
Industrial Production | $2.88$ | $2.81$ | $2.89$ | $(1.96,4.71)$ | $(1.92,4.49)$ | $(1.92,4.74)$
Interest Rate | $1.23$ | $1.22$ | $1.22$ | $(0.889,1.72)$ | $(0.891,1.73)$ | $(0.888,1.75)$
Money Stock | $3.55$ | $3.63$ | $3.56$ | $(2.26,6.09)$ | $(2.27,6.46)$ | $(2.27,6.19)$
Nominal GNP | $2.31$ | $2.32$ | $2.32$ | $(1.52,3.79)$ | $(1.53,3.75)$ | $(1.53,3.77)$
Real GNP | $3$ | $3.09$ | $3.05$ | $(1.8,5.77)$ | $(1.82,6.05)$ | $(1.8,5.85)$
Real per Capita GNP | $2.79$ | $2.84$ | $2.87$ | $(1.68,5.48)$ | $(1.65,5.68)$ | $(1.69,5.59)$
Real Wages | $5.41$ | $5.41$ | $5.38$ | $(2.97,10.4)$ | $(2.97,10.3)$ | $(2.95,10.2)$
Stock Prices | $6.11$ | $5.99$ | $5.99$ | $(3.74,10.3)$ | $(3.78,10.1)$ | $(3.72,10.2)$
Unemployment Rate | $3.49$ | $3.35$ | $3.45$ | $(2.05,6.56)$ | $(2,6.26)$ | $(2.03,6.55)$
Velocity | $3.62$ | $3.7$ | $3.73$ | $(2.26,6.94)$ | $(2.28,7.22)$ | $(2.29,7.06)$
Wages | $1.84$ | $1.84$ | $1.81$ | $(1.25,2.91)$ | $(1.25,2.95)$ | $(1.24,2.88)$

###  5.3 Sampling Efficiency

RNE of all chains is reported in Table 4. Parallels can be found to Table 2:
AA is most efficient compared to SA for large posterior medians. Although
other parts of the Gibbs sampler also affect the figures, the median $\nu$
seems to be a strong factor in the ratio between the RNE’s of AA and SA. In
both columns, the lowest value belongs to Interest Rate, and one of the
largest to Stock Prices, respectively. Moreover, the two columns predominantly
stay close in magnitude.

Table 4: RNE for estimating $\nu$ using AA, SA, and ASIS, along with the ratio between AA and SA and the mean of the three corresponding posterior median estimates. | RNE (%) |
---|---|---
| AA | SA | ASIS | AA / SA | Median
CPI | 11.1 | 4.7 | 13.2 | 2.4 | $\sim 2.18$
Employment | 10.6 | 2.8 | 10.3 | 3.7 | $\sim 2.31$
GNP Deflator | 14.6 | 8.0 | 19.9 | 1.8 | $\sim 2.21$
Industrial Production | 12.7 | 3.2 | 12.6 | 3.9 | $\sim 2.86$
Interest Rate | 7.7 | 6.6 | 7.9 | 1.2 | $\sim 1.22$
Money Stock | 20.3 | 4.4 | 23.9 | 4.6 | $\sim 3.58$
Nominal GNP | 13.8 | 7.3 | 18.0 | 1.9 | $\sim 2.32$
Real GNP | 13.8 | 4.3 | 13.1 | 3.2 | $\sim 3.05$
Real per Capita GNP | 10.7 | 4.0 | 11.6 | 2.7 | $\sim 2.84$
Real Wages | 27.6 | 4.5 | 31.5 | 6.1 | $\sim 5.4$
Stock Prices | 35.4 | 5.8 | 40.4 | 6.1 | $\sim 6.03$
Unemployment Rate | 14.7 | 4.2 | 14.0 | 3.5 | $\sim 3.43$
Velocity | 14.9 | 3.0 | 14.4 | 5.0 | $\sim 3.69$
Wages | 9.3 | 4.6 | 11.9 | 2.0 | $\sim 1.83$

##  6 Discussion

In this paper, efficient Bayesian estimation of the degrees of freedom
parameter for the Student-$t$ distribution is investigated. Specifically, the
mixing speed of AA and SA is compared through various methods: an estimate for
the expected augmented Fisher information, an extensive simulation study, and
a suitable application. AA is found to be beneficial both for real-world data
and in most synthetic setups. However, for extreme starting values, AA is
found to be completely unreliable.

As a solution, ASIS is proposed to combine the good parts of AA and SA. It is
found to be as efficient as the better of AA and SA and also invulnerable
compared to AA. Moreover, ASIS is easy to implement as soon as AA and SA have
been derived.

Lastly, it is to be highlighted that the prime interest of the present work is
the simulation-based description of the mixing speed of algorithms AA, SA, and
ASIS. In particular, the results are helpful for understanding the effect of
re-parameterization on sampling efficiency, which is independent from
execution time. Therefore, optimization for effective sampling rate, which is
a runtime-adjusted variant of RNE, and, granted, which is of most relevance to
end users but is highly dependent on many factors [22], is not considered in
the paper at hand. One would almost certainly turn away from rejection
sampling nowadays and rather experiment with other modern methods. This
remains for future research.

## Acknowledgement

We would like to thank Gregor Kastner for useful suggestions and helpful
comments.

## References

  * [1] James H. Albert and Siddhartha Chib  “Bayesian Analysis of Binary and Polychotomous Response Data”  In _Journal of the American Statistical Association_ 88.422, 1993, pp. 669–679  DOI: 10.1080/01621459.1993.10476321
  * [2] Angela Bitto and Sylvia Frühwirth-Schnatter  “Achieving Shrinkage in a Time-Varying Parameter Model Framework”  In _Journal of Econometrics_ 210.1 Elsevier B.V., 2019, pp. 75–97  DOI: 10.1016/j.jeconom.2018.11.006
  * [3] Tim Bollerslev  “A Conditionally Heteroskedastic Time Series Model for Speculative Prices and Rates of Return”  In _The Review of Economics and Statistics_ 69.3, 1987, pp. 542  DOI: 10.2307/1925546
  * [4] Paul-Christian Bürkner, Jonah Gabry, Matthew Kay and Aki Vehtari  “posterior: Tools for Working with Posterior Distributions”, 2021  URL: https://mc-stan.org/posterior/
  * [5] Siddhartha Chib, Federico Nardari and Neil Shephard  “Markov Chain Monte Carlo Methods for Stochastic Volatility Models”  In _Journal of Econometrics_ 108.2, 2002, pp. 281–316  DOI: 10.1016/S0304-4076(01)00137-3
  * [6] Siddhartha Chib, Jacek Osiewalski and Mark F.J. Steel  “Posterior Inference on the Degrees of Freedom Parameter in Multivariate $t$ Regression Models”  In _Economics Letters_ 37.4, 1991, pp. 391–397  DOI: 10.1016/0165-1765(91)90076-W
  * [7] Bruno de Finetti  “The Bayesian Approach to the Rejection of Outliers”  In _Proceedings of the Berkeley Symposium on Mathematical Statistics and Probability_ 4.1, 1961, pp. 199–210
  * [8] Donald A.. Fraser  “Necessary Analysis and Adaptive Inference”  In _Journal of the American Statistical Association_ 71.353, 1976, pp. 99–110  DOI: 10.1080/01621459.1976.10481486
  * [9] Andrew Gelman et al.  “Models for Robust Inference”  In _Bayesian Data Analyais_ New York: ChapmanHall/CRC, 2013, pp. 435–447  DOI: 10.1201/b16018
  * [10] Andrew Gelman, Walter R. Gilks and Gareth Owen Roberts  “Weak Convergence and Optimal Scaling of Random Walk Metropolis algorithms”  In _The Annals of Applied Probability_ 7.1, 1997  DOI: 10.1214/aoap/1034625254
  * [11] Andrew Gelman, Gareth Owen Roberts and Walter R. Gilks  “Efficient Metropolis Jumping Rules”  In _Bayesian Statistics 5 – Proceedings of the Fifth Valencia International Meeting_ Clarendon Press, Oxford, UK, 1994, pp. 599–607
  * [12] Stuart Geman and Donald Geman  “Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 6.6, 1984, pp. 721–741  DOI: 10.1109/TPAMI.1984.4767596
  * [13] John F. Geweke  “Bayesian Inference in Econometric Models Using Monte Carlo Integration”  In _Econometrica_ 57.6, 1989, pp. 1317  DOI: 10.2307/1913710
  * [14] John F. Geweke  “Priors for Macroeconomic Time Series and Their Application”, 1992, pp. 1–52
  * [15] John F. Geweke  “Bayesian Treatment of the Independent Student-$t$ Linear Model”  In _Journal of Applied Econometrics_ 8.S1, 1993, pp. S19–S40  DOI: 10.1002/jae.3950080504
  * [16] John F. Geweke  “Priors for Macroeconomic Time Series and Their Application”  In _Econometric Theory_ 10.3-4, 1994, pp. 609–632
  * [17] John F. Geweke  “Getting It Right: Joint Distribution Tests of Posterior Simulators”  In _Journal of the American Statistical Association_ 99.467, 2004, pp. 799–804  DOI: 10.1198/016214504000001132
  * [18] Paul Gilbert and Ravi Varadhan  “numDeriv: Accurate Numerical Derivatives”, 2019  URL: https://cran.r-project.org/package=numDeriv
  * [19] Heikki Haario, Eero Saksman and Johanna Tamminen  “An Adaptive Metropolis Algorithm”  In _Bernoulli_ 7.2, 2001, pp. 223–242  DOI: 10.2307/3318737
  * [20] Andrew C. Harvey, Esther Ruiz and Neil Shephard  “Multivariate Stochastic Variance Models”  In _The Review of Economic Studies_ 61.2, 1994, pp. 247–264  DOI: 10.2307/2297980
  * [21] Gregor Kastner and Sylvia Frühwirth-Schnatter  “Ancillarity-Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC Estimation of Stochastic Volatility Models”  In _Computational Statistics and Data Analysis_ 76 Elsevier B.V., 2014, pp. 408–423  DOI: 10.1016/j.csda.2013.01.002
  * [22] Hans-Peter Kriegel, Erich Schubert and Arthur Zimek  “The (Black) Art of Runtime Evaluation: Are We Comparing Algorithms or Implementations?”  In _Knowledge and Information Systems_ 52.2, 2017, pp. 341–378  DOI: 10.1007/s10115-016-1004-2
  * [23] Kenneth L. Lange, Roderick J.A. Little and Jeremy M.G. Taylor  “Robust Statistical Modeling Using the $t$ Distribution”  In _Journal of the American Statistical Association_ 84.408, 1989, pp. 881–896  DOI: 10.1080/01621459.1989.10478852
  * [24] Chuanhai Liu  “ML Estimation of the Multivariate $t$ Distribution and the EM Algorithm”  In _Journal of Multivariate Analysis_ 63.2, 1997, pp. 296–312  DOI: 10.1006/jmva.1997.1703
  * [25] Karla Monterrubio-Gómez et al.  “Posterior Inference for Sparse Hierarchical Non-stationary Models”  In _Computational Statistics and Data Analysis_ 148 Elsevier B.V., 2020, pp. 106954  DOI: 10.1016/j.csda.2020.106954
  * [26] Jouchi Nakajima and Yasuhiro Omori  “Leverage, Heavy-Tails and Correlated Jumps in Stochastic Volatility Models”  In _Computational Statistics and Data Analysis_ 53.6 Elsevier B.V., 2009, pp. 2335–2353  DOI: 10.1016/j.csda.2008.03.015
  * [27] Jouchi Nakajima and Yasuhiro Omori  “Stochastic Volatility Model with Leverage and Asymmetrically Heavy-Tailed Error Using GH Skew Student’s $t$ Distribution”  In _Computational Statistics and Data Analysis_ 56.11, 2012, pp. 3690–3704  DOI: 10.1016/j.csda.2010.07.012
  * [28] Charles R. Nelson and Charles I. Plosser  “Trends and Random Walks in Macroeconmic Time Series”  In _Journal of Monetary Economics_ 10.2, 1982, pp. 139–162  DOI: 10.1016/0304-3932(82)90012-5
  * [29] Omiros Papaspiliopoulos and Gareth Owen Roberts  “Stability of the Gibbs Sampler for Bayesian Hierarchical Models”  In _Annals of Statistics_ 36.1, 2008, pp. 95–117  DOI: 10.1214/009053607000000749
  * [30] Omiros Papaspiliopoulos, Gareth Owen Roberts and Martin Sköld  “A General Framework for the Parametrization of Hierarchical Models”  In _Statistical Science_ , 2007  DOI: 10.1214/088342307000000014
  * [31] Martyn Plummer, Nicky Best, Kate Cowles and Karen Vines  “CODA: Convergence Diagnosis and Output Analysis for MCMC”  In _R News_ 6.1, 2006, pp. 7–11
  * [32] R Core Team  “R: A Language and Environment for Statistical Computing”, 2021  R Foundation for Statistical Computing  URL: https://www.r-project.org/
  * [33] James O. Ramsay and Melvin R. Novick  “PLU Robust Bayesian Decision Theory: Point Estimation”  In _Journal of the American Statistical Association_ 75.372, 1980, pp. 901–907  DOI: 10.1080/01621459.1980.10477570
  * [34] Gareth Owen Roberts and Jeffrey S. Rosenthal  “Optimal Scaling for Various Metropolis-Hastings Algorithms”  In _Statistical Science_ 16.4, 2001  DOI: 10.1214/ss/1015346320
  * [35] Gareth Owen Roberts and Jeffrey S. Rosenthal  “Examples of adaptive MCMC”  In _Journal of Computational and Graphical Statistics_ 18.2, 2009, pp. 349–367  DOI: 10.1198/jcgs.2009.06134
  * [36] Radhey S. Singh  “Estimation of Error Variance in Linear Regression Models with Errors Having Multivariate Student-$t$ Distribution with Unknown Degrees of Freedom”  In _Economics Letters_ 27.1, 1988, pp. 47–53  DOI: 10.1016/0165-1765(88)90218-2
  * [37] Brajendra C. Sutradhar and Mir M. Ali  “Estimation of the Parameters of a Regression Model with a Multivariate $t$ Error Variable”  In _Communications in Statistics - Theory and Methods_ 15.2, 1986, pp. 429–450  DOI: 10.1080/03610928608829130
  * [38] Luke Tierney  “Markov Chains for Exploring Posterior Distributions”  In _The Annals of Statistics_ 22.4, 1994, pp. 1701–1728  DOI: 10.1214/aos/1176325750
  * [39] Adrian Trapletti and Kurt Hornik  “tseries: Time Series Analysis and Computational Finance”, 2020  URL: https://cran.r-project.org/package=tseries
  * [40] David A van Dyk and Xiao-Li Meng  “The Art of Data Augmentation”  In _Journal of Computational and Graphical Statistics_ 10.1, 2001, pp. 1–50  DOI: 10.1198/10618600152418584
  * [41] Aki Vehtari et al.  “Rank-Normalization, Folding, and Localization: An Improved $\hat{R}$ for Assessing Convergence of MCMC”  In _Bayesian Analysis_ Advance Publication, 2021, pp. 1–38  DOI: 10.1214/20-BA1221
  * [42] Mike West  “Outlier Models and Prior Distributions in Bayesian Linear Regression”  In _Journal of the Royal Statistical Society: Series B (Methodological)_ 46.3, 1984, pp. 431–439  DOI: 10.1111/j.2517-6161.1984.tb01317.x
  * [43] Yaming Yu and Xiao-Li Meng  “To Center or Not to Center: That Is Not the Question—An Ancillarity-Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC Efficiency”  In _Journal of Computational and Graphical Statistics_ 20.3, 2011, pp. 531–570  DOI: 10.1198/jcgs.2011.203main
  * [44] Arnold Zellner  “Bayesian and Non-Bayesian Analysis of the Regression Model with Multivariate Student-$t$ Error Terms”  In _Journal of the American Statistical Association_ 71.354, 1976, pp. 400–405  DOI: 10.1080/01621459.1976.10480357


# Improved Coefficients for the Karagiannidis–Lioumpas
Approximations and Bounds to the Gaussian $Q$-Function

Islam M. Tanash *  — and Taneli Riihonen *  — Manuscript received November 19,
2020; revised December 18, 2020 and January 8, 2021; accepted January 9, 2021.
Date of publication January 18, 2021; date of current version . This work was
partially supported by the Academy of Finland under Grant 326448. The
associate editor coordinating the review of this letter and approving it for
publication was A.-A. A. Boulogeorgos. (Corresponding author: Islam M.
Tanash.)The authors are with Tampere University, Tampere 33720, Finland
(e-mail<EMAIL_ADDRESS>taneli.riihonen@tuni.fi).Digital Object
Identiﬁer 10.1109/LCOMM.2021.3052257

###### Abstract

We revisit the Karagiannidis–Lioumpas (KL) approximation of the $Q$-function
by optimizing its coefficients in terms of absolute error, relative error and
total error. For minimizing the maximum absolute/relative error, we describe
the targeted uniform error functions by sets of nonlinear equations so that
the optimized coefficients are the solutions thereof. The total error is
minimized with numerical search. We also introduce an extra coefficient in the
KL approximation to achieve significantly tighter absolute and total error at
the expense of unbounded relative error. Furthermore, we extend the KL
expression to lower and upper bounds with optimized coefficients that minimize
the error measures in the same way as for the approximations.

###### Index Terms:

Communication theory, error probability.

##  I Introduction

Karagiannidis and Lioumpas presented in [1] a relatively tight, yet
analytically tractable, approximation for the Gaussian $Q$-function[2] as
follows:

$\displaystyle\begin{split}Q(x)&\triangleq\frac{1}{\sqrt{2\pi}}\int_{x}^{\infty}\exp\left(-{\textstyle\frac{1}{2}}t^{2}\right)\,\mathrm{d}t\\\
&\approx
a\exp\left(-b\,x^{2}\right)\cdot\frac{1-\exp\left(-c\,x\right)}{x}\triangleq\tilde{Q}(x)\end{split}$
(1)

for which their original study sets
$(a,b,c)=(\frac{1}{B\sqrt{2\pi}},\frac{1}{2},\frac{A}{\sqrt{2}})$ and proposes
for error minimization example coefficient values $A=1.98$ and $B=1.135$
rendering $(a,c)\approx(0.3515,1.4001)$.

Despite drawing some criticism [3] shortly after publication, the
‘Karagiannidis–Lioumpas (KL) approximation’ has gradually established itself
as one of the most usable substitutes for the Gaussian $Q$-function in
communication theory problems and the paper [1] has received a large number of
citations; it is only fitting to begin calling the expression after its
inventors.

A diverse set of applications for the KL approximation can be found in [4, 5,
6, 7, 8] to name but a few prominent articles. In general, the approximation
is often used in the calculation of average bit or symbol error probability as
a tractable replacement for the Gaussian $Q$-function such that analysis can
be carried out and completed in a closed form at the cost of making results
tight approximations instead of exact ones. This usually involves integrating
something like, e.g., $f(\tilde{Q}(x(\gamma)))$, where even simple functions
$f(q)$ and $x(\gamma)$, which are derived from the communication system under
study, may forbid exact analysis using the actual $Q$-function[9, 10, 11]. One
should note especially that the approximation is always used in an
intermediate step of analytical derivations and it is not meant for numerical
probability computations per se — instead, rational Chebyshev functions [12]
are perfect to that end.

This Letter is inspired by the fact that the original study [1] presents
explicit values of $a$ and $c$ for only one approximation (which has low
integrated total error when $b=\frac{1}{2}$, to be exact). However, the KL
approximation family is actually much more versatile, whereby new coefficients
can be acquired in terms of other criteria for better accuracy depending on
the application. The KL expression can be also repurposed to achieve lower and
upper bounds (that are also tight approximations) and, in certain cases,
coefficients admit explicit values. Furthermore, by introducing the extra
coefficient $b$ in (1) that originally was $b=\frac{1}{2}$ and permitting
$b<\frac{1}{2}$, we achieve significantly improved accuracy in terms of
absolute and total error.

The objective of this Letter is to apply the KL expression of the $Q$-function
to derive improved approximations and bounds which are global and tight over
$x\geq 0$ by optimizing the coefficients $(a,b,c)$ in respect to their minimum
global absolute or relative error or minimum integrated total error. Like[11]
for another popular expression[9], we present new formulation that minimizes
the maximum global error of (1) by constructing a set of equations, which
describes the corresponding error function, and solve them numerically to find
the optimized coefficients. The total error is optimized with exhaustive
search for reference. In general, when optimizing one of the three criteria,
better performance will be achieved at the expense of decreased accuracy in
terms of the others.

The new coefficients solved herein are applicable as one-to-one replacements
for the original ones of [1] adopted into the analysis of [4, 5, 6, 7, 8] and
many other studies. Literature is rich in approximations/bounds for the
$Q$-function and, typically, the application’s mathematics define, which ones
are tractable for it. Whenever (1) is preferred, our coefficients offer
variety to tailor accuracy for the application or to use bounds. The tractable
series expansion of the KL expression proposed in [13] can likewise be used
with these substitutes (but without guarantee that bounds remain global).
Consequently, they are useful also in the contexts of, e.g., [14, 15, 16, 17,
18, 19, 20] and many others.

##  II Preliminaries

The case $x\geq 0$ is presumed throughout this Letter with little loss of
generality because the relation $Q(x)=1-Q(-x)$ extends all the considered
functions to the negative real axis. In fact, this is the main motive for
optimizing approximations and bounds also subject to an additional constraint
$\tilde{Q}(0)=\frac{1}{2}$ that makes their extensions continuous at the
origin like $Q(x)$.

This study solves optimized approximations and bounds for three criteria and
for combinations thereof, viz. $\min_{(a,b,c)}d_{\mathrm{max}}$ (‘minimax
absolute error’), $\min_{(a,b,c)}r_{\mathrm{max}}$ (‘minimax relative error’)
and $\min_{(a,b,c)}d_{\mathrm{tot}}$ (‘integrated total error’ [1]), where
$\displaystyle d_{\mathrm{max}}\triangleq\max_{x\geq
0}|d(x)|,\,r_{\mathrm{max}}\triangleq\max_{x\geq
0}|r(x)|,\,d_{\mathrm{tot}}\triangleq\int_{0}^{\infty}|d(x)|\,\mathrm{d}x,$
and the error functions are defined as

$\displaystyle d(x)$ $\displaystyle\triangleq\tilde{Q}(x)-Q(x),$ (2)
$\displaystyle r(x)$
$\displaystyle\triangleq\frac{d(x)}{Q(x)}=\frac{\tilde{Q}(x)}{Q(x)}-1.$ (3)

For baseline reference, the coefficients originally given in [1] render
$d_{\mathrm{max}}\approx 0.00789$, $r_{\mathrm{max}}\approx 0.119$, and
$d_{\mathrm{tot}}\approx 0.00385$.

As implied above, the presented approximations and bounds will be global ones,
i.e., tight over the whole non-negative real axis (for all $x\geq 0$). The
error functions converge to explicit values, which may be local extrema, at
both ends of this range:

$\displaystyle\lim_{\mathclap{x\to 0}}d(x)$
$\displaystyle=ac-{\textstyle\frac{1}{2}},$
$\displaystyle\>\>\>\>\lim_{\mathclap{x\to 0}}r(x)$ $\displaystyle=2ac-1,$ (4)
$\displaystyle\lim_{\mathclap{x\to\infty}}d(x)$ $\displaystyle=0,$
$\displaystyle\>\>\>\>\lim_{\mathclap{x\to\infty}}r(x)$
$\displaystyle=\begin{cases}\infty,&\text{if }b<{\textstyle\frac{1}{2}},\\\
a\sqrt{2\pi}-1,&\text{if }b={\textstyle\frac{1}{2}},\\\ -1,&\text{if
}b>{\textstyle\frac{1}{2}}.\end{cases}$

The last limit shows especially that global approximations and bounds in terms
of relative error exist if and only if we set $b=\frac{1}{2}$. However, as a
novel fact, our study demonstrates that absolute error and total error can be
instead significantly reduced by permitting $b<\frac{1}{2}$. Therefore, two
scenarios of approximations for the absolute and total error are considered in
this Letter, i.e., approximations with $b=\frac{1}{2}$ or $b<\frac{1}{2}$.

Local error extrema may occur also at critical points, where the derivatives
of the continuous error functions vanish. Denoting differentiation with an
apostrophe, they are given by
$\displaystyle
d^{\prime}(x)=\tilde{Q}^{\prime}(x)-Q^{\prime}(x),\,\,\,r^{\prime}(x)=\frac{\tilde{Q}^{\prime}(x)Q(x)-\tilde{Q}(x)Q^{\prime}(x)}{[Q(x)]^{2}},$
where

$\displaystyle\tilde{Q}^{\prime}(x)$
$\displaystyle=-\dfrac{a\left(\left(2bx^{2}+1\right)\left(\mathrm{e}^{cx}-1\right)-cx\right)\mathrm{e}^{-bx^{2}-cx}}{x^{2}},$
(5) $\displaystyle Q^{\prime}(x)$
$\displaystyle=-\frac{1}{\sqrt{2\pi}}\exp\left(-{\textstyle\frac{1}{2}}x^{2}\right).$
(6)

Two variations of approximations are considered herein: $d(0)=r(0)=0$ and
$d(0)=-d_{\mathrm{max}}$ (resp. $r(0)=-r_{\mathrm{max}}$). The former case
maintains the continuity of the $Q$-function when extending to $x<0$ and
results in $c=\frac{1}{2\,a}$, when substituted in $\lim_{{x\to 0}}d(x)$
(resp. $\lim_{{x\to 0}}r(x)$) that is given in (4). The latter case provides
slightly better accuracy at the cost of discontinuity occurring at $x=0$ and
results in $c=\sqrt{\frac{\pi}{2}}$ in the cases of relative error, by solving
$\lim_{{x\to 0}}r(x)=-r_{\mathrm{max}}$ with
$\lim_{{x\to\infty}}r(x)=-r_{\mathrm{max}}$ that are defined in (4).

##  III Alternative Improved Coefficients for (1)

In this section, we describe the methodologies to solve the new coefficients
$(a,b,c)$ for the KL expression. They are optimized either in the minimax
sense or in terms of the integrated total error to yield an approximation, an
upper bound or a lower bound. All the 17 thus-obtained improved/alternative
coefficient sets and accuracy thereof are listed in Table I.

###  III-A Global Uniform Approximations and Bounds

The minimax optimization problems are solved in terms of both absolute and
relative errors defined in (2) and (3), respectively, by constructing a set of
nonlinear equations. This set describes the resulting error function, which
should be uniform with equal values for all the extrema points. Each extremum
point yields two equations, where one expresses its value and the other sets
the derivative of the error function to zero at that point. In addition, one
equation (for $d(x)$) or two equations (for $r(x)$) is/are obtained from
evaluating the limits at the two endpoints of the considered range,
$[0,\infty]$, per (4).

The resulting sets of equations, which have equal number of equations and
unknowns, can be solved straightforwardly by any numerical tool for the
considered variations to find the optimized sets of coefficients that satisfy
$\min_{(a,b,c)}d_{\mathrm{max}}$ for the absolute error and
$\min_{(a,b,c)}r_{\mathrm{max}}$ for the relative error. We used iteratively
random initial guesses for the unknowns in this approach, namely $(a,b,c)$,
$d_{\mathrm{max}}$ or $r_{\mathrm{max}}$, and the location of the extrema
($x_{k}$), until fsolve in Matlab converged to the solution, which is
confirmed by substitution. The formulations for the minimax
approximations/bounds are described below.

####  III-A1 Approximations in Terms of Absolute Error

The coefficients $(a,b,c)$ are optimized for approximations in terms of the
absolute error by formulating a set of equations as

$\displaystyle\begin{cases}d^{\prime}(x_{k})=0,&\text{for }k=1,2\text{ or
}1,2,3,\\\ d(x_{k})=(-1)^{k+1}\,d_{\mathrm{max}},&\text{for }k=1,2\text{ or
}1,2,3,\\\ \begin{cases}a\,c=\frac{1}{2},&\text{when }d(0)=0,\\\
a\,c=\frac{1}{2}-d_{\mathrm{max}},&\text{when
}d(0)=-d_{\mathrm{max}},\end{cases}\end{cases}$ (7)

where $x_{k}$ is an extremum point. The number of the error function’s extrema
depends on the value of $b$; if $b$ is fixed to $\frac{1}{2}$, then we have
two extrema, whereas if $b$ is allowed to be any positive value, then we need
three separate extrema.

####  III-A2 Lower Bounds in Terms of Absolute Error

For the lower bounds, we need to find the optimized coefficients which
minimize the global absolute error for $d(x)\leq 0$ when $x\geq 0$. The value
of $b$ must always equal to $\frac{1}{2}$. The tightest resulting uniform
error function will start from $d(0)=-d_{\mathrm{max}}$, with its maximum
equal to zero and its minimum equal to $-d_{\mathrm{max}}$ so that we can
formulate a set of equations as

$\displaystyle\begin{cases}d^{\prime}(x_{1})=d^{\prime}(x_{2})=0,\\\
d(x_{1})=0,\,d(x_{2})=-d_{\mathrm{max}},\\\
a\,c=\frac{1}{2}-d_{\mathrm{max}}.\end{cases}$ (8)

When $d(0)=0$, we get $a=\sqrt{\frac{\pi}{32}}$ and $c=\sqrt{\frac{8}{\pi}}$
by imposing $d^{\prime}(0)=0$ (only in this case), which produces
$a\,c^{2}=\sqrt{2/\pi}$, and solving with $c=\frac{1}{2\,a}$ that results from
setting $d(0)=0$.

####  III-A3 Upper Bounds in Terms of Absolute Error

The set of equations becomes

$\displaystyle\begin{cases}d^{\prime}(x_{1})=d^{\prime}(x_{2})=d^{\prime}(x_{3})=0,\\\
d(x_{1})=d(x_{3})=d_{\mathrm{max}},\,d(x_{2})=0,\\\
a\,c=\frac{1}{2}.\end{cases}$ (9)

In particular, we shape the uniform error function to have three extrema with
$d(x)\geq 0$ when $x\geq 0$ in which its maxima are equal to
$d_{\mathrm{max}}$ and its minimum is equal to zero. The corresponding error
function must always start from $d(0)=0$.

####  III-A4 Approximations in Terms of Relative Error

The targeted uniform error function in terms of the relative error consists of
only one maximum point and converges to $-r_{\mathrm{max}}$ as $x$ tends to
infinity, which results in $-r_{\mathrm{max}}=a\sqrt{2\,\pi}-1$ according to
(4). Therefore, we can formulate the set of equations as

$\displaystyle\begin{cases}r^{\prime}(x_{1})=0,\,r(x_{1})=r_{\mathrm{max}},\\\
\begin{cases}a\,c=\frac{1}{2},&\text{when }r(0)=0,\\\
a\,c=\frac{1-r_{\mathrm{max}}}{2},&\text{when
}r(0)=-r_{\mathrm{max}},\end{cases}\\\
a=\frac{1-r_{\mathrm{max}}}{\sqrt{2\,\pi}}.\end{cases}$ (10)

####  III-A5 Lower Bounds in Terms of Relative Error

We need to find the optimized coefficients, $a$ and $c$, in the minimax sense
for $r(x)\leq 0$ when $x\geq 0$ which converges to $-r_{\mathrm{max}}$ as $x$
tends to infinity. The resulting error function can either start from
$r(0)=-r_{\mathrm{max}}$ to formulate a set of equations as

$\displaystyle\begin{cases}r^{\prime}(x_{1})=r(x_{1})=0,\\\
a\,c=\frac{1-r_{\mathrm{max}}}{2},\,a=\frac{1-r_{\mathrm{max}}}{\sqrt{2\,\pi}},\end{cases}$
(11)

or from $r(0)=0$ yielding $a=\sqrt{\frac{\pi}{32}}$ and
$c=\sqrt{\frac{8}{\pi}}$ like with the corresponding lower bound in terms of
absolute error.

####  III-A6 Upper Bound in Terms of Relative Error

We must ensure that $r(x)\geq 0$ when $x\geq 0$ for the uniform error
function. The resulting error function has only one maximum point and
converges to zero as $x$ tends to infinity. Therefore,
$a=\frac{1}{\sqrt{2\pi}}$ and $c=\sqrt{\frac{\pi}{2}}$ as proposed earlier in
[21] and $b$ is known to be equal to $\frac{1}{2}$. The optimized upper bound
in terms of relative error is also optimal in terms of absolute error and
integrated total error for the case where $b=\frac{1}{2}$.

###  III-B Numerical Optimization in Terms of Total Error

Instead of defining
$d_{\mathrm{tot}}\triangleq\int_{0}^{R}|d(x)|\,\mathrm{d}x$ like in [1] and so
making optimized coefficients specific to the value chosen for $R$ and limited
to the range $[0,R]$, we measure total error with $R\to\infty$ and obtain
globally optimized approximations and bounds. In particular, we optimized the
coefficients for the two variations of the approximations with or without
setting $b=\frac{1}{2}$ by performing an extensive search, where we evaluated
the target metric ($d_{\mathrm{tot}}$) over wide one/two/three-dimensional
grids for the unknowns $a$, $(a,b)$, $(a,c)$, or $(a,b,c)$ with granularity of
$0.000001$ and selected the grid point with the minimum total error for each
variation. This renders four sets of optimized coefficients. Extra constraint
checks guarantee $d(x)<0$ for the lower bound and $d(x)>0$ for the upper
bound.

##  IV Numerical Results and Conclusions

TABLE I: New coefficients for (1) and approximation error thereof Figure 1:
The improved approximations and bounds compared to the KL approximation with
the original coefficients [1] and to expressions from [9] and [10].

We summarize the improved coefficients for the minimax approximations and
bounds and for the total absolute error in Table I and illustrate their error
functions in Fig. 1, together with the original KL approximation from[1] and
reference approximations and bounds from [9] and [10].111The labels having the
form $Xy$-$n$ in the results refer to the approximations and bounds as
follows: $X$ is U for upper bounds, A for approximations, and L for lower
bounds; whereas $y$ is d for absolute error, r for relative error, and t for
total error; in addition, $n$ refers to rank of the coefficients according to
the accuracy of the absolute error of each variation in an ascending order.

The numerical results show that the improved coefficients of the proposed KL
approximations and bounds are optimal subject to their optimization targets,
yet expressed precisely in implicit form as solutions to systems of nonlinear
equations as opposed to relying on numerical search to minimize error
measures. In some specific cases, a part or even all of the three coefficients
can be expressed as explicit constants. The best approximation/bound from
Table I for a specific application is chosen by contrasting requirements
against Fig. 1, provided that the KL expression (1) is suitable for it to
begin with.

As an ultimate conclusion, the presented data suggests good alternatives to
the original coefficients given in [1] for the case of $b=\frac{1}{2}$: In
some applications, the accuracy of the KL approximation might be improved by
choosing instead $A=1.95$, $B=1.113$ (a compromise between all A$y$-$n$) for
decreasing both absolute and relative error by round $15$% at the cost of
increasing total error by round $65$%; or $A=2.03$, $B=1.162$ (At-4) for
decreasing absolute error and total error by round $10$% and $25$%,
respectively, at the cost of increasing relative error by round $15$%.
Sometimes it may also be useful to choose $A=B\sqrt{\pi}\approx 1.88$,
$B=1.061$ (Ar-5) for minimizing relative error (with round $50$% reduction)
subject to zero error at the origin. In contrast, when primarily minimizing
absolute error, accuracy can be improved significantly by generalizing the KL
approximation to allow any positive $b$: Namely, the choice $a=0.32$,
$b=0.4703$, $c=1.5625$ (Ad-2) guarantees zero error at the origin while
decreasing absolute error and total error as much as round $90$% and $65$%,
respectively, at the cost of making relative error unbounded for large
arguments.

## References

  * [1] G. K. Karagiannidis and A. S. Lioumpas, “An improved approximation for the Gaussian $Q$-function,” _IEEE Commun. Lett._ , vol. 11, no. 8, pp. 644–646, Aug. 2007.
  * [2] S. Aggarwal, “A survey-cum-tutorial on approximations to Gaussian ${Q}$ function for symbol error probability analysis over Nakagami-${m}$ fading channels,” _IEEE Commun. Surveys Tuts._ , vol. 21, no. 3, pp. 2195–2223, Jul.–Sep. 2019.
  * [3] J. Dyer and S. Dyer, “Corrections to, and comments on, “An improved approximation for the Gaussian $Q$-function”,” _IEEE Commun. Lett._ , vol. 12, no. 4, p. 231, Apr. 2008.
  * [4] C. Potter, G. Venayagamoorthy, and K. Kosbar, “RNN based MIMO channel prediction,” _Signal Process._ , vol. 90, no. 2, pp. 440–450, Feb. 2010\.
  * [5] L. Tan and L. Le, “Distributed MAC protocol for cognitive radio networks: Design, analysis, and optimization,” _IEEE Trans. Veh. Technol._ , vol. 60, no. 8, pp. 3990–4003, Oct. 2011.
  * [6] J. Wu _et al._ , “Unified spectral efficiency analysis of cellular systems with channel-aware schedulers,” _IEEE Trans. Commun._ , vol. 59, no. 12, pp. 3463–3474, Dec. 2011.
  * [7] D. Malak, M. Al-Shalash, and J. Andrews, “Optimizing content caching to maximize the density of successful receptions in device-to-device networking,” _IEEE Trans. Commun._ , vol. 64, no. 10, pp. 4365–4380, Oct. 2016.
  * [8] S. Lin _et al._ , “Rayleigh fading suppression in one-dimensional optical scatters,” _IEEE Access_ , vol. 7, pp. 17 125–17 132, Jan. 2019\.
  * [9] M. Chiani, D. Dardari, and M. K. Simon, “New exponential bounds and approximations for the computation of error probability in fading channels,” _IEEE Trans. Wireless Commun._ , vol. 2, no. 4, pp. 840–845, Jul. 2003.
  * [10] M. López-Benítez and F. Casadevall, “Versatile, accurate, and analytically tractable approximation for the Gaussian $Q$-function,” _IEEE Trans. Commun._ , vol. 59, no. 4, pp. 917–922, Apr. 2011.
  * [11] I. M. Tanash and T. Riihonen, “Global minimax approximations and bounds for the Gaussian $Q$-function by sums of exponentials,” _IEEE Trans. Commun._ , vol. 68, no. 10, pp. 6514–6524, Oct. 2020.
  * [12] W. Cody, “Rational Chebyshev approximations for the error function,” _Math. Comp._ , vol. 23, no. 107, pp. 631–637, Jul. 1969.
  * [13] Y. Isukapalli and B. Rao, “An analytically tractable approximation for the Gaussian $Q$-function,” _IEEE Commun. Lett._ , vol. 12, no. 9, pp. 669–671, Sep. 2008.
  * [14] Y. Isukapalli and B. Rao, “Packet error probability of a transmit beamforming system with imperfect feedback,” _IEEE Trans. Signal Process._ , vol. 58, no. 4, pp. 2298–2314, Apr. 2010.
  * [15] Q. Zhou _et al._ , “Decode-and-forward two-way relaying with network coding and opportunistic relay selection,” _IEEE Trans. Commun._ , vol. 58, no. 11, pp. 3070–3076, Nov. 2010.
  * [16] M. Seyfi, S. Muhaidat, and J. Liang, “Performance analysis of relay selection with feedback delay and channel estimation errors,” _IEEE Signal Process. Lett._ , vol. 18, no. 1, pp. 67–70, Jan. 2011.
  * [17] M. Seyfi _et al._ , “Effect of feedback delay on the performance of cooperative networks with relay selection,” _IEEE Trans. Wireless Commun._ , vol. 10, no. 12, pp. 4161–4171, Dec. 2011.
  * [18] S. Banani and R. Vaughan, “Capacity maximisation in eigen-multiple-input multiple-output using adaptive modulation and Reed–Solomon coding,” _IET Commun._ , vol. 6, no. 15, pp. 2413–2424, Oct. 2012.
  * [19] B. Choi and L. Hanzo, “Adaptive WHT aided QAM for fading channels subjected to impulsive noise,” _IEEE Commun. Lett._ , vol. 17, no. 7, pp. 1317–1320, Jul. 2013.
  * [20] D. Basnayaka and H. Haas, “A new degree of freedom for energy efficiency of digital communication systems,” _IEEE Trans. Commun._ , vol. 65, no. 7, pp. 3023–3036, Jul. 2017.
  * [21] W. M. Jang, “A simple upper bound of the Gaussian $Q$-function with closed-form error bound,” _IEEE Commun. Lett._ , vol. 15, no. 2, pp. 157–159, Feb. 2011.


# On a class of Lebesgue-Ramanujan-Nagell equations

Azizul Hoque<EMAIL_ADDRESS>
###### Abstract.

We deeply investigate the Diophantine equation $cx^{2}+d^{2m+1}=2y^{n}$ in
integers $x,y\geq 1,m\geq 0$ and $n\geq 3$, where $c$ and $d$ are given
coprime positive integers such that $cd\not\equiv 3\pmod{4}$. We first solve
this equation for prime $n$, under the condition $n\nmid h(-cd)$, where
$h(-cd)$ denotes the class number of the quadratic field
$\mathbb{Q}(\sqrt{-cd})$. We then completely solve this equation for both $c$
and $d$ primes under the assumption that $\gcd(n,h(-cd))=1$. We also
completely solve this equation for $c=1$ and $d\equiv 1\pmod{4}$, under the
condition $\gcd(n,h(-d))=1$. For some fixed values of $c$ and $d$, we derive
some results concerning the solvability of this equation.

###### Key words and phrases:

Diophantine equation, Fibonacci number, Lehmer number, Lucas number, Primitive
divisor

###### 2010 Mathematics Subject Classification:

Primary: 11D61, 11D41; Secondary: 11B39, 11Y50

##  1\. Introduction

Many special cases of the Lebesgue-Ramanujan-Nagell type equation

$x^{2}+d=2y^{n},~{}~{}x,y\geq 1,\gcd(x,y)=1,n\geq 3,$ (1.1)

where $d\geq 1$ is an integer, have been considered over the years. In 1895,
Störmer [7] proved using an elementary factorization argument that for odd $n$
and $d=1$, (1.1) has no solution with $y>1$. Later, Cohn extended this result
in [4], and proved that $(x,y,n)=(239,13,4)$ is the only solution of (1.1)
with $y>1$ when $d=1$. In [6], Pink and Tengely considered (1.1) for $d=a^{2}$
and gave an upper bound for the exponent $n$ depending only on $a$. Moreover,
they completely solved it when $1<a<1000$ and $3<n<80$. For $d=a^{2}$ with odd
$a$ and $3\leq a\leq 501$, Tengely [8] extended the work of [6] and completely
solved (1.1) when $n$ is prime. On the other hand, for $d=q^{2m}$ with $m\geq
1$ and prime $q\geq 3$, Tengely [9] established some bounds for prime $n$.
Applying these bounds, he completely solved (1.1) for $d=q^{2m}$ and $y=17$
(resp. $d=3^{2m}$) with $m\geq 1$ and prime $q\geq 3$. Later in [12], Zhu et
al. extended the work of Tengely [9], and described all possible solutions of
(1.1) for $d=q^{2m}$ with $q$ an odd prime and $m\geq 1$. Muriefah et al.
investigated (1.1) in [1], and they completely solved it when $d\equiv
1\pmod{4}$ and less than $100$. On the other hand, Ljunggren [5] considered a
general Lebesgue-Ramanujan-Nagell equation

$cx^{2}+d=2y^{n},~{}~{}x,y\geq 1,n\geq 3,\gcd(cx,y)=1,\gcd(h(-cd),n)=1,$ (1.2)

where $c$ and $d$ are fixed positive integers, and $h(\Delta)$ denotes the
class number of the quadratic field $\mathbb{Q}(\sqrt{\Delta})$. He described
the solutions of some very particular cases of (1.2), under certain conditions
on $c,d$ and $n$.

The purpose of this paper is to perform a deep investigation on the Labesgue-
Ramanujan-Nagell type equation

$cx^{2}+d^{2m+1}=2y^{n},~{}~{}x,y\geq 1,m\geq 0,n\geq 3,$ (1.3)

where $c$ and $d$ are given positive integers with square-free $cd$. This work
generalizes some of the earlier works; for instance, the work of Ljunggren
[5]. In particular, we prove three theorems concerning the solvability of
(1.3), and then we derive some consequences for certain fixed values of $c,d$
and $n$. Note that the elegant results of Bilu et al. [2] and Voutier [10] on
the existence of primitive divisors of Lehmer numbers have turned out to be
extremely powerful tools in the proofs.

We first solve (1.3) for odd prime $n$ under the conditions $\gcd(h(-cd),n)=1$
and $cd\not\equiv 3\pmod{4}$. Precisely, we prove:

######  Theorem 1.1.

Let $c$ and $d$ be positive integers such that $cd$ is square-free and
$cd\not\equiv 3\pmod{4}$. Assume that $p$ is an odd prime such that $p\nmid
h(-cd)$.

  * (I)

The Diophantine equation

$cx^{2}+d^{2m+1}=2y^{p}$ (1.4)

has no positive integer solution $(c,d,x,y,m,p)$ for $p>7$, except for
$(x,y)=(1,1)$.

  * (II)

If $p=7$, then (1.4) has only one solution, viz. $(c,d,x,y,m)=(7,11,1169,9,0)$
provided $(x,y)\neq(1,1)$.

  * (III)

If $p=5$ and $cd\neq 1$, then for any $m\geq 0$ the solutions of (1.4) are
given by

$\displaystyle(c,d,x,y)\in$
$\displaystyle\Big{\\{}\left(F_{k-2\varepsilon}/2u^{2},L_{k+\varepsilon}/2v^{2},u\left((2F_{k}-F_{k-2\varepsilon})^{2}+F_{k}L_{k+\varepsilon}\right),F_{k}\right),$
$\displaystyle\left(L_{t-2\varepsilon}/2u^{2},5F_{t+\varepsilon}/2v^{2},u\left((2L_{t}-L_{t-2\varepsilon})^{2}+5L_{t}F_{t+\varepsilon}\right),L_{t}\right)\Big{\\}},$

where $u$ and $v$ are suitable coprime odd positive integers, $k\geq 3$ and
$t\geq 2$ are integers, $F_{k}$ (resp. $L_{k}$) is $k$-th Fibonacci (resp.
Lucas) number and $\varepsilon=\pm 1$.

  * (IV)

If $p=3$ and $cd\neq 1$, then the solutions of (1.4) are given by

$(c,d,x,y,m)=\left(c,\frac{3u^{2}c-2}{v^{2}},u(4u^{2}c-3),2u^{2}c-1,m\right),$

where $v=d^{m}$, $u\geq 1$ is an odd integers and $k\geq 2$ is an integer.

### Remark

There are many cases where both $F_{k+2}$ and $F_{k-2}$ are odd. In fact, both
$F_{k+2}$ and $F_{k-2}$ are odd when $3\mid k$. In such cases,
$c=F_{k-2\varepsilon}/2u^{2}$ can not be an integer. Therefore the solutions
mentioned in (III) may not always be possible.

Theorem 1.1 yields the following straightforward corollary.

######  Corollary 1.1.

Let $c\equiv 1\pmod{4}$ be a square-free positive integer and $p$ an odd prime
such that $p\nmid h(-c)$. Then the Diophantine equation $cx^{2}+1=2y^{p}$ has
no positive integer solution $(c,x,y,p)$, except for $(c,x,y)=(1,1,1)$.

The second result is about the positive integer solutions of (1.3), when both
$c$ and $d$ are distinct odd primes satisfying $cd\not\equiv 3\pmod{4}$ and
$n\geq 3$ is any integer other than powers of $2$. Particularly, we prove:

######  Theorem 1.2.

Let $\ell$ and $q$ be two distinct odd primes such that $\ell q\not\equiv
3\pmod{4}$. Suppose that $n>2$ is an integer which is not a power of $2$ and
$\gcd(n,h(-\ell q))=1$. Then

  * (i)

the Diophantine equation

$\ell x^{2}+q^{2m+1}=2y^{n}$ (1.5)

has no positive integer solution $(\ell,q,x,y,m,n)$ when $n$ has a prime
divisor bigger than $7$;

  * (ii)

if $7$ is the largest prime divisor of $n$, then the solutions of (1.5) are:

$(\ell,q,x,y,m,n)\in\\{(7,11,1169,9,0,7),(7,11,1169,3,0,14)\\};$

  * (iii)

if $5$ is the largest prime divisor of $n$, then the solutions (1.5) are given
by

$\hskip
36.98857pt(x,y,n)=\left(\frac{1}{4}\left|\sqrt{\frac{F_{k-2\varepsilon}}{2\ell}}\left(F_{k-2\varepsilon}^{2}-\frac{5}{2}F_{k-2\varepsilon}L_{k+\varepsilon}+\frac{5}{4}L_{k+\varepsilon}^{2}\right)\right|,\left(\frac{F_{k-2\varepsilon}+L_{k+\varepsilon}}{4}\right)^{\frac{1}{2^{r}3^{s}5^{t}}},2^{r}3^{s}5^{t+1}\right),$

where $F_{k}$ and $L_{k}$ are respectively Fibonacci and Lucas numbers
satisfying the negative Pell equation (3.13) and $\ell,q,m$ are given by
(3.11). Here $r,s$ and $t$ are suitable non-negative integers.

  * (iv)

if $3$ is the largest prime divisor of $n$, then the solutions of (1.4) are
given by

$\hskip
36.98857pt(\ell,q,x,y,m,n)=\left(\frac{q^{2m+1}+2}{3u^{2}},q,\frac{u(4q^{2m+1}-1)}{3},\left(\frac{2q^{2m+1}+1}{3}\right)^{1/2^{r}3^{t}},m,2^{r}3^{t+1}\right),$

where $u\geq 1$ is a suitable odd integer, and $r\geq 0$ and $t\geq 0$ are any
suitable integers.

We now consider $\ell\in\\{3,7,11,19\\}$ and $q\in\\{43,79,87,119,239\\}$.
Then $\ell q\not\equiv 3\pmod{4}$ and we see that $2$ and $3$ are the only
prime factors of $h(-\ell q)$. Thus as a consequence of Theorem 1.2, one gets
the following straightforward corollary.

######  Corollary 1.2.

Fix $\ell\in\\{3,7,11,19\\}$ and $q\in\\{43,79,87,119,239\\}$. If $n$ is a
positive integer such that $\gcd(n,30)=1$, then (1.5) has no positive integer
solution.

If we fix $\ell,q\in\\{3,7,11,23\\}$, then $\ell q\not\equiv 3\pmod{4}$ and
$2$ is the only prime factor of $h(-\ell d)$. Therefore by Theorem 1.2, we get
the following:

######  Corollary 1.3.

Let $\ell,q\in\\{3,7,11,23\\}$ such that $\ell\neq q$. Then for an odd integer
$n>1$, the positive integer solutions of (1.5) are given by

$(\ell,q,x,y,m,n)\in\left\\{(7,11,1169,9,0,7),\left(3,7,9,5,0,3\right)\right\\}.$

In the next result, we fix $c=1$ and completely solve (1.3) under the
restriction $\gcd(n,h(-d))=1$. This result gives a particular case of Theorem
1.1 when $n$ is a prime.

######  Theorem 1.3.

Let $d\equiv 1$ be a square-free integer with $d>1$. Assume that $n>2$ is an
integer which is not a power of $2$ such that $\gcd(n,h(-d))=1$. Then the
Diophantine equation

$x^{2}+d^{2m+1}=2y^{n},$ (1.6)

has no positive integer solution when $n$ has a prime divisor other than $2$
and $3$. When $2$ and $3$ are the only prime divisors of $n$, all positive
integer solutions are given by

$(d,x,y,m,n)=\left((3u^{2}-2)/v^{2},4u^{3}-3u,(2u^{2}-1)^{1/2^{r}3^{t}},m,2^{r}3^{t+1}\right),$

where $v=d^{m}$, $u\geq 1$ is an odd integer coprime to $v$, and $r$ and $t$
are suitable non-negative integers.

### Remarks

We make the following remarks about Theorem 1.3.

  * (i)

For $d\equiv 2\pmod{4}$, (1.6) has no solution.

  * (ii)

For $m=0$, (1.6) reduces to $x^{2}+d=2y^{n}$, which was deeply investigated by
Abu Muriefah et al. [1] for $1\leq d<100$, under the condition $\gcd(x,y)=1$.
However their method (or they) had failed to determine the following
solutions,

$(d,x,y,n)\in\\{(73,485,49,3),(73,485,7,6)\\},$

which are obtained by our Theorem 1.3

  * (iii)

Since in [1] the authors were allowed $d$ to be square, they had also failed
to determine $(d,x,y,n)=(25,99,17,3)$. Although it is out the assumptions of
Theorem 1.3, we obtain this solution from our theorem.

  * (iv)

For $(m,n)=(0,3)$, (1.6) has infinitely many positive integer solutions
$(d,x,y)$, which are given by

$(d,x,y)=(3u^{2}-2,4u^{3}-3u,2u^{2}-1),$

where $u\geq 1$ is an odd integer. It is noted that each fixed $d$ gives at
most one positive integer solution $(x,y)$.

  * (v)

From Theorem 1.3, we see that $x^{2}+d=2y^{6}$ has an infinite family of
positive integer solutions, corresponding to the solutions of the negative
Pell equation $y^{2}-2u^{2}=-1$. More precisely, these solutions are given by

$(d,x,y)=(3u_{t}^{2}-2,4u_{t}^{3}-3u_{t},y_{t}),$

where $(u_{t},y_{t})=(2y_{t-1}+3u_{t-1},3y_{t-1}+4u_{t-1})$ with $t\geq 0$ an
integer and $(u_{0},y_{0})=(1,1)$.

##  2\. Lehmer sequences and some lemmas

We need some definitions and notations to state some crucial results on the
existence of primitive divisors of Lehmer numbers. Suppose that $\alpha$ and
$\beta$ are algebraic integers. Then $(\alpha,\beta)$ is a Lehmer pair if
$(\alpha+\beta)^{2}$ and $\alpha\beta$ are two non-zero coprime rational
integers, and $\alpha/\beta$ is not a root of unity. For a positive integer
$\ell$ and a Lehmer pair $(\alpha,\beta)$, the corresponding sequence of
Lehmer numbers is defined by

$\mathcal{L}_{\ell}(\alpha,\beta)=\begin{cases}\dfrac{\alpha^{\ell}-\beta^{\ell}}{\alpha-\beta}&\text{
if }\ell\text{ is odd},\vspace{1mm}\\\
\dfrac{\alpha^{\ell}-\beta^{\ell}}{\alpha^{2}-\beta^{2}}&\text{ if }\ell\text{
is even}.\end{cases}$

Note that all Lehmer numbers are non-zero rational integers. Recall that two
Lehmer pairs $(\alpha_{1},\beta_{1})$ and $(\alpha_{2},\beta_{2})$ are
equivalent if $\alpha_{1}/\alpha_{2}=\beta_{1}/\beta_{2}\in\\{\pm
1,\pm\sqrt{-1}\\}$. For such pairs,
$\mathcal{L}_{\ell}(\alpha_{1},\beta_{1})=\pm\mathcal{L}_{\ell}(\alpha_{2},\beta_{2})$
for any positive integer $\ell$. A prime number $p$ is a primitive divisor of
$\mathcal{L}_{\ell}(\alpha,\beta)$ if $p\mid\mathcal{L}_{\ell}(\alpha,\beta)$
and
$p\nmid(\alpha^{2}-\beta^{2})^{2}\mathcal{L}_{1}(\alpha,\beta)\mathcal{L}_{2}(\alpha,\beta)\cdots\mathcal{L}_{\ell-1}(\alpha,\beta)$.

We now define $a:=(\alpha+\beta)^{2}$ and $b:=a-4\alpha\beta$. Then
$\alpha=(\sqrt{a}\pm\sqrt{b})/2$ and $\beta=(\sqrt{a}\mp\sqrt{b})/2$. This
pair $(a,b)$ is known as the parameter of the Lehmer pair $(\alpha,\beta)$. We
combine a classical result of Voutier [10, Theorem 1] and a result of Bilu et
al. [2, Theorem 1.3]; precisely Tables 2 and 4 in [2], to get the following:

######  Theorem A.

Let $\ell$ be a prime number such that $3\leq\ell\leq 30$. If the Lehmer
numbers $\mathcal{L}_{\ell}(\alpha,\beta)$ have no primitive divisor, then up
to equivalence, the parameters $(a,b)$ of the corresponding Lehmer pair
$(\alpha,\beta)$ are given by the following:

  * (i)

For $\ell=3,(a,b)=\begin{cases}(1+t,1-3t)\text{ with }t\neq 1,\\\
(3^{k}+t,3^{k}-3t)\text{ with }t\not\equiv
0\pmod{3},(k,t)\neq(1,1);\end{cases}$

  * (ii)

For
$\ell=5,(a,b)=\begin{cases}(F_{k-2\varepsilon},F_{k-2\varepsilon}-4F_{k})\text{
with }k\geq 3,\\\ (L_{k-2\varepsilon},L_{k-2\varepsilon}-4L_{k})\text{ with
}k\neq 1;\end{cases}$

  * (iii)

For $\ell=7,(a,b)=(1,-7),(1,-19),(3,-5),(5,-7),(13,-3),(14,-22)$;

  * (iv)

For $\ell=13,(a,b)=(1,-7)$;

where $t$ is any non-zero integer, $k$ is non-negative integer,
$\varepsilon=\pm 1$, and $F_{k}$ (resp. $L_{k}$) denotes the $k$-th Fibonacci
(resp. Lucas) number.

We also need the following result of Bilu et al. [2, Theorem 1.4] in the
proofs.

######  Theorem B.

For every integer $\ell>30$, the Lehmer numbers
$\mathcal{L}_{\ell}(\alpha,\beta)$ have primitive divisors.

Let $F_{k}$ (resp. $L_{k}$) denote the $k$-th term in the Fibonacci (resp.
Lucas) sequence defined by $F_{0}=0,F_{1}=1$, and $F_{k+2}=F_{k}+F_{k+1}$
(resp. $L_{0}=2,L_{1}=1$, and $L_{k+2}=L_{k}+L_{k+1}$), where $k\geq 0$ is an
integer. Then $L_{k}=F_{k-1}+F_{k+1}$ and $5F_{k}=L_{k-1}+L_{k+1}$ for all
$k\geq 1$. Applying these relations, we can derive the following:

######  Lemma 2.1.

For an integer $k\geq 0$, let $F_{k}$ (resp. $L_{k}$) denote the $k$-th
Fibonacci (resp. Lucas) number. Then for $\varepsilon=\pm 1$,

  * (i)

$4F_{k}-F_{k-2\varepsilon}=L_{k+\varepsilon}$,

  * (ii)

$4L_{k}-L_{k-2\varepsilon}=5F_{k+\varepsilon}$

######  Theorem C ([3, Theorems 2 and 4]).

Let $F_{k}$ and $L_{k}$ be as in Lemma 2.1. Then

  * (i)

if $L_{k}=2x^{2}$, then $k=0,6$;

  * (ii)

if $F_{k}=2x^{2}$, then $k=0,3,6$.

The following classical lemma is a special case of [11, Corollary 3.1].

######  Lemma 2.2.

Let $c$ and $d$ be positive integers such that $cd$ is square-free and
$cd\not\equiv 3\pmod{4}$. Assume that $n$ is a positive odd integer such that
$\gcd(n,h(-cd))=1$. Then all the positive integer solutions $(X,Y,Z)$ of the
equation

$cX^{2}+dY^{2}=2Z^{n},~{}~{}\gcd(cX,dY)=1,$ (2.1)

can be expressed as

$\frac{X\sqrt{c}+Y\sqrt{-d}}{\sqrt{2}}=\delta\left(\frac{u\sqrt{c}\pm
v\sqrt{-d}}{\sqrt{2}}\right)^{n},$

where $u$ and $v$ are positive integers satisfying
$Z=\dfrac{u^{2}c+v^{2}d}{2}$ and $\gcd(uc,vd)=1$, and $\delta\in\\{-1,1\\}$.

##  3\. Proofs

###### Proof of Theorem 1.1.

Suppose that $(c,d,x,y,m,p)$ is a positive integer solution of (1.4). We first
observe that $2\nmid cdx$ since $cd$ is square-free. Also reading (1.4) modulo
$4$ and then applying $cd\not\equiv 3\pmod{4}$, we see that $y$ is odd.

Since $cd$ is square-free, so that $\gcd(cx,d)=1$, and thus by Lemma 2.2, we
can write:

$\frac{x\sqrt{c}+d^{m}\sqrt{-d}}{\sqrt{2}}=\delta\left(\frac{u\sqrt{c}\pm
v\sqrt{-d}}{\sqrt{2}}\right)^{p},$ (3.1)

where $u$ and $v$ are positive integers satisfying

$y=\dfrac{u^{2}c+v^{2}d}{2}$ (3.2)

and $\gcd(uc,vd)=1$, and $\delta\in\\{-1,1\\}$. Note that $2\nmid uv$ since
$2\nmid cdxy$.

Assume that $\alpha=\dfrac{u\sqrt{c}+v\sqrt{-d}}{\delta\sqrt{2}}$ and
$\bar{\alpha}=\dfrac{u\sqrt{c}-v\sqrt{-d}}{\delta\sqrt{2}}$. Then using (3.2),
we check that both $\alpha$ and $\bar{\alpha}$ are algebraic integers. Also,
$(\alpha+\bar{\alpha})^{2}=2u^{2}c$ and
$\alpha\bar{\alpha}=\dfrac{u^{2}c+v^{d}}{2}=y$ are coprime positive integers
since $\gcd(2u^{2}c,y)=1$. We see that $\alpha/\bar{\alpha}$ is a root of

$yZ^{2}-(u^{2}c-v^{2}d)Z+y=0.$

As $\gcd(y,u^{2}c-v^{2}d)=\gcd(uc,vd)=1$, so that $\alpha/\bar{\alpha}$ is not
an algebraic integer and thus it not a root of unity. Therefore
$(\alpha,\bar{\alpha})$ is a Lehmer pair with the corresponding parameters
$(2u^{2}c,-2v^{2}d)$.

Suppose that $\mathcal{L}_{\ell}(\alpha,\bar{\alpha}),\ell\in\mathbb{N}$, is
the Lehmer number corresponding to the Lehmer pair $(\alpha,\bar{\alpha})$.
Then using (3.1), we get

$\mathcal{L}_{p}(\alpha,\bar{\alpha})=\frac{\alpha^{p}-\bar{\alpha}^{p}}{\alpha-\bar{\alpha}}=\frac{\delta
d^{m}}{v}.$ (3.3)

This shows that all the prime divisors of
$\mathcal{L}_{p}(\alpha,\bar{\alpha})$ are also the divisors of $d$. On the
other hand $(\alpha^{2}-\bar{\alpha}^{2})^{2}=-2u^{2}v^{2}cd$, which shows
that all the divisors of $d$ are also divisors of
$(\alpha^{2}-\bar{\alpha}^{2})^{2}$. Therefore
$\mathcal{L}_{p}(\alpha,\bar{\alpha})$ has no primitive divisor.

Since $(2u^{2}c,-2v^{2}d)$ is the parameter corresponding to the Lehmer pair
$(\alpha,\bar{\alpha})$, so that by Theorems A and B there is no Lehmer number
without primitive divisor for $p>7$. This shows that (1.4) has no positive
integer solution when $p>7$.

We now treat the remaining cases , viz. $p=3,5,7$, individually. For $p=7$, by
Theorem A the only possibility is $(2u^{2}c,-2v^{2}d)=(14,-22)$, which gives
$(c,d,u,v)=(7,11,1,1)$. Using this in (3.2), we get $y=9$. Thus (1.4) reduces
to $7x^{2}+11^{2m+1}=2\times 9^{7}$. We use MAGMA to solve it, which gives
$(x,m)=(1169,0)$. Therefore the corresponding solution of (1.4) is
$(c,d,x,y,m,p)=(7,11,1169,9,0,7)$.

We now consider the case $p=5$. In this case, Theorem A gives

$\begin{cases}F_{k-2\varepsilon}=2u^{2}c,\\\
F_{k-2\varepsilon}-4F_{k}=-2v^{2}d,\end{cases}$ (3.4)

and

$\begin{cases}L_{t-2\varepsilon}=2u^{2}c,\\\
L_{t-2\varepsilon}-4L_{t}=-2v^{2}d,\end{cases}$ (3.5)

where $k\geq 3$ and $t\geq 2$ are integers, and $\varepsilon=\pm 1$. Now (3.4)
gives $2F_{k}=u^{2}c+v^{2}d$, and thus by (3.2), we get $F_{k}=y$ with $k\geq
3$. Also (3.4) and Lemma 2.1 together give $c=F_{k-2\varepsilon}/2u^{2}$ and
$d=L_{k+\varepsilon}/2v^{2}$.

We now compare the real parts in (3.1) for $p=5$ to get

$x=\frac{\delta u}{4}(u^{4}c^{2}-10u^{2}v^{2}cd+5v^{4}d^{2}).$

Applying (3.4) and then simplifying, we have

$x=\delta
u\left((2F_{k}-F_{k-2\varepsilon})^{2}+4F_{k}^{2}-F_{k}F_{k-2\varepsilon}\right).$

We apply Lemma 2.1 (precisely (i)) to get

$x=u\left((2F_{k}-F_{k-2\varepsilon})^{2}+F_{k}L_{k+\varepsilon}\right).$

As in the previous case, (3.1), (3.2), (3.5) and (ii) of Lemma 2.1 all
together give

$(c,d,x,y)=\left(L_{t-2\varepsilon}/2u^{2},5F_{t+\varepsilon}/2v^{2},u\left((2L_{t}-L_{t-2\varepsilon})^{2}+5L_{t}F_{t+\varepsilon}\right),L_{t}\right),$

where $t\geq 2$.

We now consider the remaining case, $p=3$. In this case, we get by Theorem A:

$\begin{cases}2u^{2}c=a+1,\\\ 2v^{2}d=3a-1,\end{cases}$ (3.6)

and

$\begin{cases}2u^{2}c=b+3^{k},\\\ 2v^{2}d=3b-3^{k},\end{cases}$ (3.7)

where $a>1,b>1$ and $k>1$ are integers with $b\not\equiv 0\pmod{3}$. The
equations (3.2) and (3.6) together give $d=(3u^{2}c-2)/v^{2}$ and
$y=2u^{2}c-1$. Also comparing the real and imaginary parts of (3.1), and then
applying (3.6), we get $x=u(4u^{2}c-3)$ and $2d^{m}=\delta(3u^{2}cv-dv^{3})$.
Since $v$ is a divisor of $d^{m}$, so that $v=d^{r}$ for $0\leq r\leq m$. Thus
we can write the last equation as

$2d^{m-r}=\delta(3u^{2}c-d^{2r+1}).$ (3.8)

This shows that $d=1,3$, except for $r=m$. If $d=1$, then (3.6) implies
$3a-1=2$, which gives $a=1$ and thus $c=1$. Therefore, we have
$(c,d,x,y)=(1,1,1,1)$.

Again for $d=3$, (3.8) can be written as $2\times
3^{m-r-1}=\delta(u^{2}c-3^{2r})$, which implies either $r=0$ or $r=m-1$. If
$r=0$, then $v=1$ and thus as before, we get $(c,d,x,y)=(1,1,1,1)$. For
$r=m-1$, (3.6) implies $2\times 3^{2m-1}=3a-1$, which is not possible. Finally
for $r=m$, we have

$(c,d,x,y)=\left(c,(3u^{2}c-2)/v^{2},u(4u^{2}c-3),2u^{2}c-1\right),$

where $v=d^{m}$ and $u\geq 1$ is any odd integer.

We now compare the imaginary parts in (3.1), and then use (3.7) to get $d=3$
and $m=k>1$. Thus using (3.2) and (3.7), we get $y=2v^{2}+3^{k-1}$ and
$c=(v^{2}+2\times 3^{k-1})/u^{2}$. Finally, equating the real parts in (3.1)
and using these values, we have $x=u(4v^{2}-2\times 3^{k-1})$. This is not
possible since $x$ is odd. ∎

###### Proof of Theorem 1.2.

Let $p$ be a prime divisor of $n$. Then (1.5) can be written as

$\ell x^{2}+q^{2m+1}=2Y^{p},$ (3.9)

where $Y=y^{n/p}$. Thus by Theorem 1.1, (3.9) has no solution when $p>7$, and
therefore (1.5) has no solution when $n$ has a prime divisor $>7$.

For $p=7$, Theorem 1.1 gives $(\ell,q,x,y,m)=(7,11,1169,9,0)$. Since
$Y=y^{n/7}=9$, so that $(y,n)=(9,7),(3,14)$.

We now consider the case $p=5$. As in (3.1)), we get

$\frac{x\sqrt{\ell}+q^{m}\sqrt{-q}}{\sqrt{2}}=\delta\left(\frac{u\sqrt{\ell}\pm
v\sqrt{-q}}{\sqrt{2}}\right)^{5},$ (3.10)

where $u$ and $v$ are positive integers satisfying
$Y=y^{n/5}=\dfrac{u^{2}\ell+v^{2}q}{2}$ with $\gcd(uc,vd)=1$, and
$\delta\in\\{-1,1\\}$. Also analogous to (3.3), one gets
$\mathcal{L}_{5}(\alpha,\bar{}\alpha)=q^{m}/v$, which implies $v=q^{r}$ for
some integer $0\leq r\leq m$. We first assume that $r\leq m-1$. Then the
imaginary parts of (3.10) gives $q=5$ since $\gcd(q,u\ell)=1$. Upon
simplifying these imaginary parts, we get

$u^{4}\ell^{2}-5^{2r}\times 2u^{2}\ell+5^{4r+1}=5^{m-r-1}\times 4\delta,$

which implies that $r=m-1$, and thus $v=5^{m-1}$. Therefore the above equation
can be written as

$u^{4}\ell^{2}-5^{2m-2}\times 2u^{2}\ell+5^{4m-3}=4\delta.$

This can be simplified to

$\left(\frac{u^{2}\ell-5^{2m-2}}{2}\right)^{2}+5^{4m-4}=\pm 1.$

Reading this modulo $4$, we can further reduce it to

$\left(\frac{u^{2}\ell-5^{2m-2}}{2}\right)^{2}+5^{4m-4}=1,\text{ with
}\ell\equiv 1\pmod{4}.$

This can be written as

$\left(\frac{2+5^{2m-2}-u^{2}\ell}{2}\right)\left(\frac{2+u^{2}\ell-5^{2m-2}}{2}\right)=5^{4m-4}.$

Clearly, $(2+5^{2m-2}-u^{2}\ell)/2$ and $(2+u^{2}\ell-5^{2m-2})/2$ are
coprime, so that either $(2+5^{2m-2}-u^{2}\ell)/2=1$ or
$(2+u^{2}\ell-5^{2m-2})/2=1$. Both of them lead to $u^{2}\ell=5^{2m-2}$, which
is not possible as $\gcd(u\ell,5)=1$.

When $r=m$, we have $|\mathcal{L}_{5}(\alpha,\bar{}\alpha)|=1$ and thus it has
no primitive divisor. Therefore by Theorem A and Lemma 2.1, we get

$(F_{k-2\varepsilon},L_{k+\varepsilon})=(2u^{2}\ell,2q^{2m+1}),$ (3.11)

or

$(L_{k-2\varepsilon},5F_{k+\varepsilon})=(2u^{2}\ell,2q^{2m+1}).$ (3.12)

Now (3.12) gives $q=5$, and hence $F_{k+\varepsilon}=2\times 5^{2m}$. Thus by
Theorem C, $(k,\varepsilon,m)=(4,-1,0),$ $(2,1,0)$, and hence using (3.12), we
get $(u,\ell)=(1,3)$. Therefore $Y=(u^{2}\ell+v^{2}q)/2$ implies that $Y=4$,
and thus (3.9) becomes

$3x^{2}+5=2\times 4^{5}.$

This equation has no integer solution.

Now the imaginary parts of (3.10) gives

$5u^{4}\ell^{2}-10u^{2}\ell q^{2m+1}+q^{4m+2}=\pm 4.$

This can be written as

$\left(\frac{q^{2m+1}-5u^{2}\ell}{2}\right)^{2}-5u^{4}\ell^{2}=\pm 1.$

By reading this equation modulo $4$, it can be reduced to

$\left(\frac{q^{2m+1}-5u^{2}\ell}{2}\right)^{2}-5u^{4}\ell^{2}=-1.$

Employing (3.11), we get

$\left(\frac{L_{k+\varepsilon}-5F_{k-2\varepsilon}}{4}\right)^{2}-\frac{5}{4}F_{k-2\varepsilon}^{2}=-1.$
(3.13)

Again equating the real parts of (3.10), we get

$4x=|u(u^{4}\ell^{2}-10u^{2}\ell v^{2}q+5v^{4}q^{2})|.$

Applying (3.11), this can be written as

$x=\frac{1}{4}\left|\sqrt{\frac{F_{k-2\varepsilon}}{2\ell}}\left(F_{k-2\varepsilon}^{2}-\frac{5}{2}F_{k-2\varepsilon}L_{k+\varepsilon}+\frac{5}{4}L_{k+\varepsilon}^{2}\right)\right|.$

Now $y^{n/5}=(u^{2}\ell+v^{2}q)/2$ gives

$y=\left(\frac{F_{k-2\varepsilon}+L_{k+\varepsilon}}{4}\right)^{5/n}.$

Finally for $p=3$, we follow the previous approach to get

$(\ell,q,x,Y)=\left(\frac{q^{2m+1}+2}{3u^{2}},q,\frac{u(4q^{2m+1}-1)}{3},\frac{2q^{2m+1}+1}{3}\right).$

This implies

$(\ell,q,x,y)=\left(\frac{q^{2m+1}+2}{3u^{2}},q,\frac{u(4q^{2m+1}-1)}{3},\left(\frac{2q^{2m+1}+1}{3}\right)^{1/2^{r}3^{t}}\right),$

where $r\geq 0$ and $t\geq 0$ are suitable integers. ∎

###### Proof of Theorem 1.3.

The proof of this theorem is mutatis mutandis the same as that of Theorem 1.2
when $n$ has a prime divisor $>5$. For the remaining cases, we rewrite (1.6)
as

$x^{2}+d^{2m+1}=2Y^{p},$ (3.14)

where $p\leq 5$ is a prime divisor of $n$ and $Y=y^{n/p}$. Note that both $x$
and $y$ are odd as $d\equiv 1\pmod{4}$.

Assume that $p=5$. Then by Theorem 1.1, $F_{k-2\varepsilon}=2u^{2}$ and
$L_{{t-2\varepsilon}}=2u^{2}$. Since $u\geq 1$ is odd, $k\geq 3$ and $t\geq
2$, so that by Theorem C we get:

$(k,\varepsilon,u)=(5,1,1)\text{ and
}(t,\varepsilon,u)\in\left\\{(2,1,1),(4,-1,3),(8,1,3)\right\\}.$

Also by Theorem 1.1, we get $L_{k+\varepsilon}=2v^{2}d$ and
$5F_{t+\varepsilon}=2v^{2}d$, which further imply:

$(k,\varepsilon,u,v,d)=(5,1,1,3,1)\text{ and
}(t,\varepsilon,u,v,d)\in\left\\{(2,1,1,1,5),(4,-1,3,1,5),(8,1,3,1,85)\right\\}.$

Utilizing these values in

$x=\begin{cases}u\left((2F_{k}-F_{k-2\varepsilon})^{2}+F_{k}L_{k+\varepsilon}\right),\\\
u\left((2L_{t}-L_{t-2\varepsilon})^{2}+5L_{t}F_{t+\varepsilon}\right),\end{cases}$

we see that $x$ is even in each case. This is not possible, and thus (3.14)
has no solution when $p=5$.

We now assume that $p=3$. In this case, we get by Theorem 1.1:

$(x,Y,d,m)=(4u^{3}-3u,2u^{2}-1,(3u^{2}-2)/v^{2},m),$

where $v=d^{m}$ and $u\geq$ is an odd integer. Thus the corresponding
solutions of (1.6) are given by
$(x,y,d,m)=\left(4u^{3}-3u,(2u^{2}-1)^{3/n},(3u^{2}-2)/v^{2},m\right)$, where
$n=3t$ with $t$ a positive integer having only $1,2$ or $3$ its possible
divisors. ∎

## acknowledgement

The author would like to thank Professor Srinivas Kotyada and The Institute of
Mathematical Sciences, Chennai for hospitality and support during the period
when this work was started. The author thanks Professor Kalyan Chakraborty and
Professor Alain Togbé for comments and suggestions which improved the quality
of the paper. This work is supported by the grants SERB MATRICS Project No.
MTR/2017/00100 and SERB-NPDF (PDF/2017/001958), Govt. of India.

## References

  * [1] F. S. Abu Muriefah, F. Luca, S. Siksek and Sz. Tengely, On the Diophantine equation $x^{2}+C=2y^{n}$, Int. J. Number Theory 5 (2009), no. 6, 1117–1128.
  * [2] Y. Bilu, G. Hanrot and F. M. Voutier, Existence of primitive divisors of Lucas and Lehmer numbers (with an appendix by M. Mignotte), J. Reine Angew. Math. 539 (2001), 75–122.
  * [3] J. H. E. Cohn, Square Fibonacci numbers, etc., Fibonacci Quart. 2 1964, no. 2, 109–113.
  * [4] J. H. E. Cohn, Perfect Pell powers, Glasgow Math. J. 38 (1996),19–20.
  * [5] W. Ljunggren, On the Diophantine equation $Cx^{2}+D=2y^{n}$, Math. Scand. 18 (1966), 69–86.
  * [6] I. Pink and Sz. Tengely, Full powers in arithmetic progressions, Publ. Math. Debrecen 57 (2000, 535–545.
  * [7] C. Störmer, Solution complète en nombres entiers $m,n,x,y,k$ de l’équation $m\arctan(x^{-1})+n\arctan(y^{-1})=k\pi/4$, Norske Videnskabsselsk. Skr. Kristiania 1 (1895), no. 11, 21 pp.
  * [8] Sz. Tengely, On the Diophantine equation $x^{2}+a^{2}=2y^{p}$, Indag. Math. (N. S.) 15 (2004), no. 2, 291–304.
  * [9] Sz. Tengely, On the Diophantine equation $x^{2}+q^{2m}=2y^{p}$, Acta Arith. 127 (2007), no. 1, 71–86.
  * [10] P. M. Voutier, Primitive divisors of Lucas and Lehmer sequences, Math. Comp. 64 (1995), 869–888.
  * [11] P. -Z. Yuan, On the Diophantine equation $ax^{2}+by^{2}=ck^{n}$, Indag. Math. (N. S.) 16 (2005), no. 2, 301–320.
  * [12] H. Zhu, M. -H. Le and A. Togbé, On the exponential Diophantine equation $x^{2}+p^{2m}=2y^{n}$, Bull. Aust. Math. Soc. 86 (2012), 303–314.