equation

$\BODY$ (1)

# Objective Bayesian analysis for the differential entropy of the Gamma
distribution

Eduardo Ramos1, Osafu A. Egbon1, Pedro L. Ramos1,
Francisco A. Rodrigues1 and Francisco Louzada1
1Institute of Mathematical Science and Computing, University of São Paulo, São
Carlos, Brazil

###### Summary

The use of entropy related concepts goes from physics, such as in statistical
mechanics, to evolutionary biology. The Shannon entropy is a measure used to
quantify the amount of information in a system, and its estimation is usually
made under the frequentist approach. In the present paper, we introduce a
fully objective Bayesian analysis to obtain this measure’s posterior
distribution. Notably, we consider the gamma distribution, which describes
many natural phenomena in physics, engineering, and biology. We reparametrize
the model in terms of entropy, and different objective priors are derived,
such as Jeffreys prior, reference prior, and matching priors. Since the
obtained priors are improper, we prove that the obtained posterior
distributions are proper and their respective posterior means are finite. An
intensive simulation study is conducted to select the prior that returns
better results in terms of bias, mean square error, and coverage
probabilities. The proposed approach is illustrated in two datasets, where the
first one is related to the Achaemenid dynasty reign period, and the second
data describes the time to failure of an electronic component in the sugarcane
harvest machine.

##  1 Introduction

In recent years, there has been a growing interest in estimating different
metrics of information theory related to parametric distributions. The Shannon
entropy, also known as differential entropy, introduced by Claude Shannon
[30], is an essential quantity that measures the amount of available
information or uncertainty outcome of a random process. Given a density
function $f(x|\alpha,\beta)$, the differential entropy is given by

$H(\alpha,\beta)=\mathbb{E}\left(-\log f(x|\alpha,\beta)\right).$ (2)

The differential entropy depends on the distribution parameters, and, given a
sample, it is necessary to be estimated. The commonly used method to estimate
the parameters is the maximum likelihood approach due to its one-to-one
invariance property. Hence, we need only to estimate the parameters of the
original model and plug-in the entropy function. Under this approach many
authors have derived the estimators for different distributions such as,
Weibull [9], Inverse Weibull [33], Log-logistic [13] and for the exponential
distribution with different shift origin [20], to list a few.

A major drawback of the maximum likelihood inference is that the obtained
estimates are usually biased for small samples [11]. Another concern happens
under small samples when constructing the confidence intervals for the
parameters since such intervals are not precise and may not return good
coverage probabilities. In this case, the maximum likelihood estimation (MLE)
skewness study is essential to assess the quality of the interval [12]. To
overcome these limitations, we can use objective Bayesian methods. In this
context, the inference for the parameters of the gamma distribution have been
discussed earlier under this approach by Miller [23], Sun and Ye [31], Berger
et al. [6], and Louzada and Ramos [21]. Moreover, Ramos et al. [27] revised
the most common objective priors and provided sufficient and necessary
conditions for the obtained posteriors and their higher moments to be proper.

Although the authors have obtained different joint posterior distributions for
the parameters of interest, the obtained posterior means can not be directly
plunged in the Shannon entropy. Under the Bayesian approach, it is necessary
to obtain the posterior distribution of the entropy measure. In this context,
Shakhatreh [29] recently derived different posterior distributions using
objective priors for the entropy assuming a Weibull distribution. On the other
hand, the cited distribution’s entropy expression is not as complicated as the
gamma distribution’s entropy expression. With this in mind, in this paper,
focusing on the gamma distribution, we derive the posterior distributions
using objective priors, such as Jeffreys prior [18], reference priors [7, 5,
6], and matching priors [32], and prove that the obtained posteriors are
proper and can be used to construct the posterior distributions of the Shannon
entropy. Moreover, even if the posterior distribution is proper, the posterior
mean can be infinite, which is undesirable, and thus we shall also prove that
the obtained posterior means for the entropy measure are finite. Finally, the
credibility intervals are obtained to construct accurate interval estimates.

The gamma distribution considered here is a two-parameter family of
distribution among the most well-known distribution used to model different
stochastic processes and to make statistical inferences, and has received
attention from different fields. It surfaces in many areas of applications,
including financial analysis [10], climate analysis [17], reliability analysis
[16], machine learning [19], and physics [14]. Particularly, the gamma
distribution includes the exponential distribution, Erlang distribution, and
chi-square distribution as special cases.

A random variable $X$ follows a gamma distribution, if its probability density
function, parametrized by a shape parameter $\alpha>0$ and scale parameter
$\beta>0$, is given by,

$f(x\,|\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta
x},~{}~{}x>0,$ (3)

where $\Gamma(\phi)=\int_{0}^{\infty}{e^{-x}x^{\phi-1}dx}$ is the gamma
function.

The paper is organized as follows. Section 2 presents the maximum likelihood
estimators for the gamma distribution parameters and the Shannon Entropy
computation. Section 3 presents the objective Bayesian analysis using
objective priors for the Shannon entropy parameter’s reparametrized posterior
distribution. Section 4 provides a simulation study to select the best
objective prior. In Section 5, the methodology is illustrated on a real
dataset. Some final comments are given in Section 6.

##  2 Frequentist approach

The classical inference (frequentist) is a commonly used approach to conduct
parameter estimation of a particular distribution. In this case, the parameter
is treated as fixed, and the MLE is commonly used to obtain the estimates. The
MLE has good asymptotic properties, such as invariance, consistency, and
efficiency. This procedure search the parameter space of $\boldsymbol{\theta}$
where the maximum likelihood
$\hat{\boldsymbol{\theta}}=\sup_{\boldsymbol{\theta}}L(\boldsymbol{\theta}|x)$
is obtained. Here our main aim is to obtain the estimate of a function of the
parameters. Hence, firstly we need to obtain the entropy measure,
mathematically defined as
$H(\boldsymbol{\theta})=E(-logf(x|\boldsymbol{\theta}))$, which quantifies the
amount of uncertainty in the data $x$. Besides, it should be noted that a
higher realization of $H$ indicates more uncertainty.

The entropy $H$ of the gamma density is given by

$\displaystyle H(\alpha,\beta)$
$\displaystyle=-\int^{\infty}_{0}\log\left(\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}exp\\{-\beta
x\\}\right)f(x|\alpha,\beta)dx$ (4)
$\displaystyle=\alpha-\log(\beta)+\log(\Gamma(\alpha))+(1-\alpha)\psi(\alpha),$

where $\psi(k)=\frac{\partial}{\partial k}\log\Gamma(x)$ is the digamma
function.

Now, consider a change of variable by setting $W=\alpha$, which implies
$H=W-\log(\beta)+\log\Gamma(W)+(1-W)\psi(W)$. The aim of the transformation is
to obtain a likelihood of $H$ and $W$ instead of $\alpha$ and $\beta$.
Therefore, if $X_{1}$, $\ldots$, $X_{n}\,$ are a complete sample from (3) then
the likelihood function of $H$ and $W$ is given as

$L(W,H\,|\boldsymbol{x})=\frac{\delta(W,H)^{nW}}{\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}\exp\left\\{-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\},$
(5)

where $\delta(W,H)=\exp\left(W+\log\Gamma(W)+(1-W)\psi(W)-H\right)$.

The log-likelihood function is given by

$\displaystyle
l(W,H\,|\boldsymbol{x})=W\log(\delta(W,H))-\log\Gamma(W)+W\sum_{i=1}^{n}\log(x_{i})-\delta(W,H)\sum_{i=1}^{n}x_{i}.$
(6)

The MLEs for the parameters are obtained by directly maximizing the log-
likelihood function $\ell(\lambda,\phi;\boldsymbol{t})$. Hence, after some
algebraic manipulations the MLEs $\hat{W}$ and $\hat{H}$ are obtained from the
solution of

$\frac{\partial l(W,H\,|\boldsymbol{x})}{\partial
W}=\log(\delta(W,H))-\psi(W)-\sum_{i=1}^{n}\log(x_{i})+\sigma\left(W-\delta(W,H)\sum_{i=1}^{n}x_{i}\right)$
$\frac{\partial l(W,H\,|\boldsymbol{x})}{\partial
H}=-W+\delta(W,H)\sum_{i=1}^{n}x_{i}$

where $\sigma=1+(1-W)\psi^{\prime}(W)$. The solutions for these equations
provide the maximum likelihood estimators for the entropy of the gamma
distributions, $\widehat{H}$ and $\widehat{W}$. Since equation (3) cannot be
solved easily using a closed-form solution, numerical techniques must estimate
the true parameters.

Following [22], the MLEs are asymptotically normally distributed with a joint
bivariate normal distribution given by

$\left(\hat{W}_{\rm MLE},\hat{H}_{\rm MLE}\right)\sim
N_{2}\left[\left(W,H\right),I^{-1}\left(W,H\right)\right]\quad\mbox{ as }\quad
n\to\infty,$

where $I(W,H)$ is the Fisher information matrix for the reparametrized model
given by

$\displaystyle\begin{aligned}
I(W,H)=&\begin{bmatrix}\psi^{\prime}(W)-2\sigma+W\sigma^{2}&1-\sigma W\\\
1-\sigma W&W\end{bmatrix},\end{aligned}$ (7)

and $\psi\,^{\prime}$($W$) is the derivative of $\psi$($W$), called the
trigamma function.

In the present paper, we are only interested in $H$, and thus, given $0<a<1$
and using the element $(I(W,H)^{-1})_{22}$, we can conclude that the
confidence interval for the estimate of the entropy measure with a confidence
level of $100(1-a)\%$ for $H$ is given by

$\hat{H}-{Z_{\frac{a}{2}}}\sqrt{(1-W)^{2}\psi^{\prime}(W)+2-W}<H<\hat{H}+{Z_{\frac{a}{2}}}\sqrt{(1-W)^{2}\psi^{\prime}(W)+2-W},$
(8)

where $a$ is the significance level and ${Z_{\frac{a}{2}}}$ is the
$\frac{a}{2}$-th percentile of the standard normal distribution.

##  3 Bayesian Inference

Here, the parameter $\boldsymbol{\theta}$ is considered as a random variable
and the distribution that represents knowledge about $\boldsymbol{\theta}$ is
refereed as a prior distribution and defined by $\pi(\boldsymbol{\theta})$.
The distribution $\pi(\boldsymbol{\theta})$ provides the knowledge or
uncertainty about $\boldsymbol{\theta}$ before obtaining the sample data
$\boldsymbol{x}$. After the data $x$ is observed, a natural way of combining
the resulting information from the a priori the distribution and the
likelihood function is done by the Bayes’ theorem, resulting in the posterior
distribution of $\boldsymbol{\theta}$ given $\boldsymbol{x}$. In a Bayesian
framework, Ramos $et\;al.$ [27] analyzed the properties of the posterior
distribution of the gamma distribution parameters and stated the conditions
for this distribution to have proper posterior and finite moments.

To obtain the posterior distributions for the $H$ parameter, we can consider
the one-to-one invariance property of the Jeffreys prior, reference prior, and
matching prior, and thus we only need to obtain the Jacobian matrix related to
the reparametrization from $\alpha$ and $\beta$ to $H$ and $W$. After some
algebraic manipulations, we can conclude that the parameters $\beta$ and
$\alpha$ can be written as

$\beta=\exp\left(W+\log(\Gamma(W))+(1-W)\psi(W)-H\right)\quad\mbox{ and
}\quad\alpha=W,$

and thus, from the relations

$\displaystyle\frac{\partial\alpha}{\partial
H}=0,\,\frac{\partial\alpha}{\partial W}=1,\,\frac{\partial\beta}{\partial
H}=-\beta\ \ \mbox{ and }\ \ \frac{\partial\beta}{\partial
W}=\left(1+(1-W)\psi^{(1)}(W)\right)\beta,$

it follows that the Jacobian matriz (J) relative to the change of variable
will be given by

$\displaystyle J=\begin{bmatrix}\frac{\partial\alpha}{\partial
H}&\frac{\partial\alpha}{\partial W}\\\ \frac{\partial\beta}{\partial
H}&\frac{\partial\beta}{\partial W}\end{bmatrix}=\begin{bmatrix}0&1\\\
-\beta&\sigma\beta\end{bmatrix},$ (9)

where $\sigma=1+(1-W)\psi^{\prime}(W)$.

The use of objective priors plays an essential role in Bayesian analysis where
the data provide the dominant information, and the posterior distribution is
not overshadowed by prior information. Such priors allow us to conduct
objective Bayesian inference. On the other hand, in most situations, they are
not proper prior distributions and may lead to improper posterior,
invalidating the analysis since we cannot compute the normalizing constant.
Therefore, we need to check if the obtained posterior (and posterior mean) is
proper (or finite). The priors for the entropy and its related posterior
distributions will be discussed in the next subsections.

Before we derive the priors and posterior distributions, hereafter, we shall
always assume that there are at least two distinct data $t_{i}$, that is,
there exists $1\leq i<j\leq n$ such that $t_{i}\neq t_{j}$. Additionally,
before we proceed, we present below a definition and proposition that will be
used to prove that the obtained posteriors are proper. In the following let
$\overline{\mathbb{R}}=\mathbb{R}\cup\\{-\infty,\infty\\}$ denote the extended
real number line and let $\mathbb{R}^{+}$ denote the strictly positive real
numbers. The following definition is a special case from the one presented in
[25] and will play an important role in proving that the analyzed posterior
distributions and posterior means are proper.

######  Definition 3.1.

Let $a\in\mathbb{\overline{R}}$,
$\operatorname{g}:\mathcal{U}\to\mathbb{R^{+}}$ and
$\operatorname{h}:\mathcal{U}\to\mathbb{R^{+}}$, where
$\mathcal{U}\subset\mathbb{R}$ and suppose that $\lim_{x\to
a}\dfrac{\operatorname{g}(x)}{\operatorname{h}(x)}=c\in\mathbb{R}$. Then, if
$c>0$, we say that $g(x)\underset{x\to a}{\propto}h(x)$.

Regarding the above definition, we have the following proposition from [25].

######  Proposition 3.2.

Let $\operatorname{g}:(a,b)\to\mathbb{R^{+}}$ and
$\operatorname{h}:(a,b)\to\mathbb{R^{+}}$ be continuous functions in
$(a,b)\subset\mathbb{R}$, where $a\in\overline{\mathbb{R}}$ and
$b\in\overline{\mathbb{R}}$, and let $c\in(a,b)$. Then
$\operatorname{g}(x)\underset{x\to a}{\propto}\operatorname{h}(x)$ implies in
$\int_{a}^{c}g(t)\;dt\propto\int_{a}^{c}h(t)\;dt$ and
$\operatorname{g}(x)\underset{x\to b}{\propto}\operatorname{h}(x)$ implies in
$\int_{c}^{b}g(t)\;dt\propto\int_{c}^{b}h(t)\;dt$.

###  3.1 Jeffreys prior

Jeffreys [18] described a procedure to achieve an objective prior, which is
invariant under one-to-one monotone transformations. The invariant property of
the Jeffreys prior has been widely exploited to make statistical inferences
from its posterior distribution numerical analysis. The prior construction is
based on the square root of the determinant of the Fisher information matrix
$I(\alpha,\beta)$. Thus, the Jeffreys prior to the gamma distribution is given
by

$\pi_{1}\left(\alpha,\beta\right)\propto\frac{\sqrt{\alpha\psi^{\prime}(\alpha)-1}}{\beta}.$
(10)

Additionally, from the determinant of the Fisher information, or using the
change of variables over the Jeffreys prior we have

$\pi_{1}\left(H,W\right)\propto\sqrt{W\psi^{\prime}(W)-1}.$ (11)

Finally, the joint posterior distribution for $H$ and $W$ produced by the
Jeffreys prior is

$\displaystyle\pi_{1}(H,W|\boldsymbol{x})\propto\frac{\delta(W,H)^{nW}\sqrt{W\psi^{\prime}(W)-1}}{\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}\exp\left\\{-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\}.$
(12)

######  Theorem 3.3.

The posterior density (12) is proper for all $n\geq 2$.

###### Proof.

Using the change of variables $\exp(-H)=u\Leftrightarrow du=-\exp(-H)dH$ and
denoting $\delta_{1}(W)=\exp(W+\log(\Gamma(W))+(1-W)\psi(W))$ it follows that

$\displaystyle d_{1}(x)$
$\displaystyle\propto\int_{0}^{\infty}\int_{-\infty}^{\infty}\pi_{1}(H,W|\boldsymbol{x})\;dHdW$
$\displaystyle\propto\int_{0}^{\infty}\int_{0}^{\infty}\frac{\delta_{1}(W)^{nW}u^{nW-1}\sqrt{W\psi^{\prime}(W)-1}}{\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}\exp\left\\{-\delta_{1}(W)u\sum_{i=1}^{n}x_{i}\right\\}dudW$
$\displaystyle=\int_{0}^{\infty}\frac{\delta_{1}(W)^{nW}\sqrt{W\psi^{\prime}(W)-1}}{\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}\int_{0}^{\infty}u^{nW-1}\exp\left\\{-\delta_{1}(W)\left(\sum_{i=1}^{n}x_{i}\right)u\right\\}dudW$
$\displaystyle=\int_{0}^{\infty}\sqrt{W\psi^{\prime}(W)-1}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}dW=\int_{0}^{1}g_{1}(W)dw+\int_{1}^{\infty}g_{1}(W)dw,$

where
$g_{1}(W)=\sqrt{W\psi^{\prime}(W)-1}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}>0$
for all $W\in(0,\infty)$. Now, according to [25, 28], we have
$\frac{\Gamma(nW)}{\Gamma(W)^{n}}\underset{W\to 0^{+}}{\propto}W^{n-1}$ and
$\sqrt{W\psi^{\prime}(W)-1}\underset{W\to 0^{+}}{\propto}W^{-1/2}$ and since

$\lim_{W\to
0^{+}}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}=1\Rightarrow\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\underset{W\to
0^{+}}{\propto}1,$

it follows by Proposition 3.2 that

$\int_{0}^{1}g_{1}(W)dW\propto\int_{0}^{1}W^{-1/2}\times 1\times
W^{n-1}\;dW<\infty.$

Moreover, due to [25, 28] we have
$\frac{\Gamma(nW)}{\Gamma(W)^{n}}\underset{W\to\infty}{\propto}n^{nW}W^{(n-1)/2}$
and $\sqrt{W\psi^{\prime}(W)-1}\underset{W\to\infty}{\propto}W^{-1/2}$, and
since $x_{i}$ are not all equal, due to the inequality of the arithmetic and
geometric means we have
$q=\log\left(\frac{\frac{1}{n}\sum_{i=1}^{n}x_{i}}{\sqrt[n]{\prod_{i=1}^{n}{x_{i}}}}\right)>0$
and thus it follows that

$\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}=\left(\frac{\frac{1}{n}\sum_{i=1}^{n}x_{i}}{\sqrt[n]{\prod_{i=1}^{n}{x_{i}}}}\right)^{-nW}n^{-nW}=\exp\left(-nqW\right)n^{-nW}.$

Therefore, from Proposition 3.2 it follows that

$\displaystyle\int_{1}^{\infty}g_{1}(W)dW$
$\displaystyle\propto\int_{1}^{\infty}W^{-1/2}\times\exp\left(-nqW\right)n^{-nW}\times
n^{nW}W^{(n-1)/2}dW$
$\displaystyle=\int_{1}^{\infty}W^{n/2-1}\exp\left(-nqW\right)\;dW=\frac{\Gamma(n/2)}{(nq)^{n/2}}<\infty,$

which concludes the proof. ∎

######  Theorem 3.4.

The posterior mean of $H$ relative to (12) is finite for any $n\geq 2$.

###### Proof.

Doing the change of variables $\exp(-H)=u\Leftrightarrow du=-\exp(-H)dH$ and
denoting $\delta_{1}(W)=\exp(W+\log(\Gamma(W))+(1-W)\psi(W))$, it follows that

$\displaystyle E_{1}[H|x]$
$\displaystyle\propto\int_{0}^{\infty}\int_{-\infty}^{\infty}H\pi_{1}(\alpha,\beta|\boldsymbol{x})\;dHdW$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{\infty}-\log(u)\frac{\delta_{1}(W)^{nW}u^{nW-1}\sqrt{W\psi^{\prime}(W)-1}}{\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}\exp\left\\{-\delta_{1}(W)u\sum_{i=1}^{n}x_{i}\right\\}dudW$
$\displaystyle=\int_{0}^{\infty}\frac{\delta_{1}(W)^{nW}\sqrt{W\psi^{\prime}(W)-1}}{\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}\int_{0}^{\infty}\left(-\log(u)\right)u^{nW-1}\exp\left\\{-\delta_{1}(W)\left(\sum_{i=1}^{n}x_{i}\right)u\right\\}dudW.$

Moreover, from the identity
$\psi(z)\Gamma(z)=\Gamma^{\prime}(z)=\int_{0}^{\infty}\log(t)t^{z-1}e^{-t}dz$
one obtains that

$\int_{0}^{\infty}\log(s)s^{z-1}e^{-as}ds=1/a^{z}\int_{0}^{\infty}\log(t/a)t^{z-1}e^{-t}dt=1/a^{z}\left(\psi(z)\Gamma(z)-\log(a)\Gamma(z)\right)$

and thus, letting $\left|\cdot\right|$ denote the absolute value operator and
letting
$\delta_{2}(W)=\left|\psi(nW)\right|+\left|\log(\Gamma(W))\right|+(1+W)|\psi(W)|+W+\left|\log\left(\sum_{i=1}^{n}x_{i}\right)\right|$
for all $W>0$, and using the triangle inequality we have

$\displaystyle\left|E_{1}[H|x]\right|$
$\displaystyle\propto\left|\int_{0}^{\infty}\left(\psi(nW)-\log\left(\delta_{1}(W)\sum_{i=1}^{n}x_{i}\right)\right)\sqrt{W\psi^{\prime}(W)-1}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}dW\right|$
$\displaystyle\leq\int_{0}^{\infty}\left|\psi(nW)-\log\left(\delta_{1}(W)\sum_{i=1}^{n}x_{i}\right)\right|\sqrt{W\psi^{\prime}(W)-1}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}dW$
$\displaystyle\leq\int_{0}^{\infty}\delta_{2}(W)\sqrt{W\psi^{\prime}(W)-1}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}dW=\int_{0}^{1}h_{1}(W)dW+\int_{1}^{\infty}h_{1}(W)dW,$

where
$h_{1}(W)=\delta_{2}(W)\sqrt{W\psi^{\prime}(W)-1}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}$
for all $W>0$.

We shall now prove that $\delta_{2}(W)\underset{W\to 0^{+}}{\propto}W^{-1}$
and $\delta_{2}(W)\underset{W\to\infty}{\propto}W\log(W)$. Indeed, notice that
$\delta_{2}(W)\geq W>0$ for $W>0$. Moreover, since due to Abramowitz [1] we
have $\lim_{W\to 0^{+}}W\Gamma(W)=1$ and $\lim_{W\to 0^{+}}-W\psi(W)=1$ it
follows that

$\displaystyle\lim_{W\to 0^{+}}\frac{\left|\psi(nW)\right|}{W^{-1}}=\lim_{W\to
0^{+}}\frac{1}{n}\left|(nW)\psi(nW)\right|=\frac{1}{n}$
$\displaystyle\lim_{W\to
0^{+}}\frac{\left|\log\left(\Gamma(W)\right)\right|}{W^{-1}}=\lim_{W\to
0^{+}}\left|W\log(\Gamma(W))-W\log(W)\right|=\left|0\cdot\log(1)-0\right|=0$
$\displaystyle\lim_{W\to
0^{+}}\frac{(1+W)\left|\psi(W)\right|}{W^{-1}}=\lim_{W\to
0^{+}}(1+W)\left|W\psi(W)\right|=1\mbox{ and }$ $\displaystyle\lim_{W\to
0^{+}}\frac{W+\left|\log\left(\sum_{i=1}^{n}x_{i}\right)\right|}{W^{-1}}=\lim_{W\to
0^{+}}\left(W^{2}+W\left|\log\left(\sum_{i=1}^{n}x_{i}\right)\right|\right)=0$

and thus

$\lim_{W\to
0^{+}}\frac{\delta_{2}(W)}{W^{-1}}=\frac{1}{n}+1\Rightarrow\delta_{2}(W)\underset{W\to
0^{+}}{\propto}\frac{1}{W}.$

On the other hand, since due to Abramowitz [1] we have
$\lim_{W\to\infty}\frac{\psi(W)}{\log(W)}=1$, it follows from the L’hopital
rule that

$\lim_{W\to\infty}\frac{\log(\Gamma(W))}{W(\log(W)+1)}=\lim_{W\to\infty}\frac{(\log(\Gamma(W))^{\prime}}{(W(\log(W)+1))^{\prime}}=\lim_{W\to\infty}\frac{\psi(W)}{\log(W)}=1,$

and therefore, considering $W\geq 1$ we have

$\displaystyle\lim_{W\to\infty}\frac{\left|\psi(W)\right|}{W(\log(W)+1)}=\lim_{W\to\infty}\frac{1}{W}\frac{1}{(1+\log(W)^{-1})}\left|\frac{\psi(W)}{\log(W)}\right|=0,$
$\displaystyle\lim_{W\to\infty}\frac{\left|\log(\Gamma(W))\right|}{W(\log(W)+1)}=\lim_{W\to\infty}\left|\frac{\log(\Gamma(W))}{W(\log(W)+1)}\right|=1,$
$\displaystyle\lim_{W\to\infty}\frac{(1+W)\left|\psi(W)\right|}{W(\log(W)+1)}=\lim_{W\to\infty}\left(1+W^{-1}\right)\frac{1}{\left(1+\log(W)^{-1}\right)}\left|\frac{\psi(W)}{\log(W)}\right|=1,\mbox{
and }$
$\displaystyle\lim_{W\to\infty}\frac{W+\left|\log\left(\sum_{i=1}^{n}x_{i}\right)\right|}{W(\log(W)+1)}=\lim_{W\to\infty}\left(\frac{1}{\log(W)+1}+\frac{\left|\log\left(\sum_{i=1}^{n}x_{i}\right)\right|}{W(\log(W)+1)}\right)=0,$

and thus

$\lim_{W\to\infty}\frac{\delta_{2}(W)}{W(\log(W)+1)}=2\Rightarrow\delta_{2}(W)\underset{W\to\infty}{\propto}W\log(W).$

Therefore, combining the obtained proportionality $\delta_{2}(W)\underset{W\to
0^{+}}{\propto}W^{-1}$ with the proportionalities proved in Theorem 3.3 and
using Proposition 3.2 we have

$\int_{0}^{1}h_{1}(W)dW\propto\int_{0}^{1}W^{-1}\times W^{-1/2}\times 1\times
W^{n-1}\;dW<\infty.$

Finally, using the proportionality
$\delta_{2}(W)\underset{W\to\infty}{\propto}W\log(W)$, letting $q>0$ be as in
the proof of Theorem 3.3 and using that $\log(W)+1\leq\exp(\log(W))=W$ for
$W\geq 1$, it follows from the proportionalities proved during Theorem 3.3 and
from Proposition 3.2 that

$\displaystyle\int_{1}^{\infty}h_{1}(W)dW$
$\displaystyle\propto\int_{1}^{\infty}W(\log(W)+1)\times
W^{-1/2}\times\exp\left(-nqW\right)n^{-nW}\times n^{nW}W^{(n-1)/2}dW$
$\displaystyle\leq\int_{1}^{\infty}W^{(n/2+2)-1}\exp\left(-nqW\right)\;dW=\frac{\Gamma(n/2+2)}{(nq)^{n/2+2}}<\infty,$

which concludes the proof. ∎

In order to sample for the posterior distribution we obtain that the marginal
posterior distributions of $W$ is given by

$\pi_{1}(W|\boldsymbol{x})\propto\sqrt{W\psi^{\prime}(W)-1}\frac{\Gamma(nW)}{\Gamma(W)^{n}}\left(\frac{\sqrt[n]{\prod_{i=1}^{n}{x_{i}}}}{\sum_{i=1}^{n}x_{i}}\right)^{nW},$

and the conditional posterior distribution of $H$ is given by

$\displaystyle\pi_{1}(H|W,\boldsymbol{x})\propto\exp\left\\{-nWH-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\}.$

###  3.2 Reference prior

Bernardo [7] discussed a different approach to obtain a new class of objective
priors, named as reference priors. Further, many studies were presented to
develop formal and rigorous definitions to derive such class of prior
distributions under different contexts [2, 3, 4, 5, 6]. The reference prior is
obtained by maximizing the Kullback-Leibler (KL) divergence assuming some
regularity conditions. The idea of the expected posterior information to the
prior allows the data to have the maximum influence on the posterior
distributions. The reference priors have essential properties such as
consistent sampling, consistent marginalization, and one-to-one transformation
invariance [8]. The reference priors may depend on the order of the parameters
of interest. Hence, for the gamma distribution, we have two distinct priors
that are presented below.

####  3.2.1 Reference prior when $\beta$ is the parameter of interest

The reference prior when $\beta$ is the parameter of interest and $\alpha$ is
the nuisance parameter is given by

$\pi_{2}\left(\alpha,\beta\right)\propto\frac{\sqrt{\psi^{\prime}(\alpha)}}{\beta}.$
(13)

Thus, using the Jacobian transformation it follows that the related reference
prior is given by

$\pi_{2}\left(W,H\right)\propto\sqrt{\psi^{\prime}(W)}.$ (14)

Finally, the joint posterior distribution for $H$ and $W$, produced by the
reference prior (21), is given by

$\pi_{2}(W,H|\boldsymbol{x})\propto\delta(W,H)^{nH}\frac{\sqrt{\psi^{\prime}(W)}}{\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{H}}\right\\}\exp\left\\{-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\}.$
(15)

######  Theorem 3.5.

The posterior density (15) is proper for all $n\geq 2$.

###### Proof.

Doing the change of variables $\exp(-H)=u\Leftrightarrow du=-\exp(-H)dH$,
denoting $\delta_{1}(W)=\exp(W+\log(\Gamma(W))+(1-W)\psi(W))$ and proceeding
analogously as in the proof of Theorem 3.3 we have

$\displaystyle d_{2}(w)\propto$
$\displaystyle\int_{0}^{\infty}\int_{-\infty}^{\infty}\pi_{2}(W,H|\boldsymbol{x})\;dHdW\propto\int_{0}^{1}g_{2}(W)dw+\int_{1}^{\infty}g_{2}(W)dw,$

where
$g_{2}(W)=\sqrt{\psi^{\prime}(W)}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}>0$
for all $W\in(0,\infty)$. Now, according to [25, 28], we have
$\frac{\Gamma(nW)}{\Gamma(W)^{n}}\underset{W\to 0^{+}}{\propto}W^{n-1}$ and
$\sqrt{\psi^{\prime}(W)}\underset{W\to 0^{+}}{\propto}W^{-1}$, and since we
proved in Theorem 3.3 that
$\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\underset{W\to
0^{+}}{\propto}1$, it follows from Proposition 3.2 that

$\int_{0}^{1}g_{2}(W)dW\propto\int_{0}^{1}W^{-1}\times 1\times
W^{n-1}\;dW<\infty.$

Moreover, from Abramowitz [1] we have
$\sqrt{\psi^{\prime}(W)}\underset{W\to\infty}{\propto}W^{-1/2}$, which
combined with
$\sqrt{W\psi^{\prime}(W)-1}\underset{W\to\infty}{\propto}W^{-1/2}$ implies in
$\sqrt{\psi^{\prime}(W)}\underset{W\to\infty}{\propto}\sqrt{W\psi^{\prime}(W)-1}$.
Therefore it follows that $g_{2}(W)\underset{W\to\infty}{\propto}g_{1}(W)$.
and by Proposition 3.2 it follows that

$\displaystyle\int_{1}^{\infty}g_{2}(W)dW\propto\int_{1}^{\infty}g_{1}(W)\;dW<\infty,$

which concludes the proof. ∎

######  Theorem 3.6.

The posterior mean of $H$ relative to (15) is finite for all $n\geq 2$.

###### Proof.

Proceeding analogously as in the proof of Theorem 3.4 it follows that

$\displaystyle\left|E_{2}[H|x]\right|$
$\displaystyle\propto\int_{0}^{\infty}\left|\int_{-\infty}^{\infty}H\pi_{2}(H,W|\boldsymbol{x})\;dHdW\right|$
$\displaystyle\leq\int_{0}^{\infty}\delta_{2}(W)\sqrt{\psi^{\prime}(W)}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}dW=\int_{0}^{1}h_{2}(W)dW+\int_{1}^{\infty}h_{2}(W)dW,$

where $\delta_{2}(W)$ is the same as defined in the proof of Theorem 3.4 and

$\displaystyle
h_{2}(W)=\delta_{2}(W)\sqrt{\psi^{\prime}(W)}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}\cdot$

Since in the proof of Theorem 3.4 we showed that $\delta_{2}(W)\underset{W\to
0^{+}}{\propto}W^{-1}$, together with the proportionalities proved in Theorem
3.3 and Proposition 3.2 we have

$\int_{0}^{1}h_{2}(W)dW\propto\int_{0}^{1}W^{-1}\times W^{-1}\times 1\times
W^{n-1}\;dW<\infty.$

Finally, from the proof of Theorem 3.5 we know that
$\sqrt{\psi^{\prime}(W)}\underset{W\to\infty}{\propto}\sqrt{W\psi^{\prime}(W)-1}$,
which implies directly that $h_{2}(W)\underset{W\to\infty}{\propto}h_{1}(W)$,
and thus from Proposition 3.2 it follows that

$\displaystyle\int_{1}^{\infty}h_{2}(W)dW\propto\int_{1}^{\infty}h_{1}(W)\;dW<\infty,$

which concludes the proof. ∎

The marginal posterior distributions of $W$ is given by

$\pi_{2}(W|\boldsymbol{x})\propto\sqrt{\psi^{\prime}(W)}\frac{\Gamma(nW)}{\Gamma(W)^{n}}\left(\frac{\sqrt[n]{\prod_{i=1}^{n}{x_{i}}}}{\sum_{i=1}^{n}x_{i}}\right)^{nW}.$

Moreover, the conditional posterior distribution of $H$ is given by

$\displaystyle\pi_{2}(H|W,\boldsymbol{x})\propto\exp\left\\{-nWH-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\}.$

####  3.2.2 Reference prior when $\alpha$ is the parameter of interest

The reference prior when $\alpha$ is the parameter of interest and $\beta$ is
the nuisance parameter is given by

$\pi_{3}\left(\alpha,\beta\right)\propto\frac{1}{\beta}\sqrt{\frac{\alpha\psi^{\prime}(\alpha)-1}{\alpha}}.$
(16)

Therefore, in terms of the reparametrized model, the reference prior when $W$
is the parameter of interest and $H$ is the nuisance parameter is given by

$\pi_{3}\left(W,H\right)\propto\sqrt{\frac{W\psi^{\prime}(W)-1}{W}}.$ (17)

Finally, the joint posterior distribution for $\alpha$ and $\beta$, produced
by the reference prior (17) is given by

$\pi_{3}(W,H|\boldsymbol{x})\propto\sqrt{\frac{W\psi^{\prime}(W)-1}{W}}\frac{\delta(W,H)^{nW}}{\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}\exp\left\\{-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\}.$
(18)

######  Theorem 3.7.

The posterior density (18) is proper for all $n\geq 2$.

###### Proof.

Since
$\sqrt{\frac{W\psi^{\prime}(W)-1}{W}}=\sqrt{\psi^{\prime}(W)-\frac{1}{W}}\leq\sqrt{\psi^{\prime}(W)}$
it follows that $\pi_{3}(W,H)\leq\pi_{2}(W,H)$ for all $W\in(0,\infty)$ and
$H\in(-\infty,\infty)$ and thus Theorem 3.7 follows directly from Theorem 3.5.
∎

######  Theorem 3.8.

The posterior mean of $H$ relative to (18) is finite for all $n\geq 2$.

###### Proof.

Since
$\sqrt{\frac{W\psi^{\prime}(W)-1}{W}}=\sqrt{\psi^{\prime}(W)-\frac{1}{W}}\leq\sqrt{\psi^{\prime}(W)}$,
it follows that
$|H\pi_{3}(W,H)|=|H|\pi_{3}(W,H)\leq|H|\pi_{2}(W,H)=|H\pi_{2}(W,H)|$ for all
$W\in(0,\infty)$ and $H\in(-\infty,\infty)$ and thus Theorem 3.8 follows
directly from Theorem 3.6 ∎

The marginal posterior distributions of $W$ is given by

$\pi_{3}(W|\boldsymbol{x})\propto\sqrt{\frac{W\psi^{\prime}(W)-1}{W}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}\left(\frac{\sqrt[n]{\prod_{i=1}^{n}{x_{i}}}}{\sum_{i=1}^{n}x_{i}}\right)^{nW}.$

Moreover, the conditional posterior distribution of $H$ is given by

$\displaystyle\pi_{3}(H|W,\boldsymbol{x})\propto\exp\left\\{-nWH-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\}.$

###  3.3 Matching priors

Tibshirani [32] considered a different method to obtain a class of one
parameter non-informative prior distribution with nuisance parameters. Letting
$\pi(\theta_{1},\theta_{2})$ be a prior distribution with the parameter of
interest $\theta_{1}$ and a nuisance $\theta_{2}$, the proposed approach
requires that the resulting credible interval of the posterior distribution
for $\theta_{1}$ have a frequentist coverage accurate to $O$($n^{-1}$), that
is, it requires that

$P\left[\theta_{1}\leq\theta_{1}^{1-\alpha}(\pi;X)|(\theta_{1},\theta_{2})\right]=1-\alpha-O(n^{-1}),$
(19)

where $\theta_{1}^{1-\alpha}(\pi;X)|(\theta_{1},\theta_{2})$ denotes the
$(1-\alpha)$th quantile of the posterior distribution of $\theta_{1}$. The
priors that satisfy (19) up to $O(n^{-1})$ are know as matching priors. Under
parametric orthogonality, Mukerjee & Dey [24] discussed sufficient and
necessary conditions for a class of Tibshirani priors to be a matching prior
up to $o(n^{-1})$. Sun and Ye [31] derived a Berger and Bernardo’s [2] forward
and backward reference prior for a two-parameter exponential family, and
further showed that the reference prior are special cases of the matching
priors. For a gamma distribution, they showed that the reference prior (17) is
a matching prior when $\beta$ is set as a nuisance parameter and $\alpha$ is
the interest parameter, and proved that there exist no matching prior up to
order $O(n^{-1})$. Again, they showed that the reference prior (17) is a
matching prior when $\beta$ is the interest parameter and $\alpha$ is the
nuisance parameter with order $O(n^{-1})$ and proved that there exists a
matching prior up to order $o(n^{-1})$. The cited matching prior is defined as

$\pi_{4}\left(\alpha,\beta\right)\propto\frac{\alpha\psi^{\prime}(\alpha)-1}{\beta\sqrt{\alpha}}.$
(20)

Thus, the reparametrized version of the proposed matching prior is given by

$\pi_{4}\left(H,\delta(W,H)\right)\propto\frac{W\psi^{\prime}(W)-1}{\sqrt{W}}.$
(21)

Finally, the joint posterior distribution for $H$ and $W$, produced by the
matching prior (21) is given by

$\pi_{4}(W,H|\boldsymbol{x})\propto\delta(W,H)^{nH}\frac{(W\psi^{\prime}(W)-1)}{\sqrt{W}\,\Gamma(W)^{n}}\left\\{\prod_{i=1}^{n}{x_{i}^{H}}\right\\}\exp\left\\{-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\}.$
(22)

######  Theorem 3.9.

The posterior density (22) is proper for all $n\geq 2$.

###### Proof.

Doing the change of variables $\exp(-H)=u\Leftrightarrow du=-\exp(-H)dH$,
denoting $\delta_{1}(W)=\exp(W+\log(\Gamma(W))+(1-W)\psi(W))$ and proceeding
analogously as in the proof of Theorem 3.3 we have

$\displaystyle d_{4}(w)\propto$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\pi_{4}(W,H|\boldsymbol{x})\;dHdW\propto\int_{0}^{1}g_{4}(W)dw+\int_{1}^{\infty}g_{4}(W)dw,$

where
$g_{4}(W)=\frac{(W\psi^{\prime}(W)-1)}{\sqrt{W}}\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}>0$
for all $W\in(0,\infty)$. Now, according to [25, 28], we have
$\frac{\Gamma(nW)}{\Gamma(W)^{n}}\underset{W\to 0^{+}}{\propto}W^{n-1}$ and
$\sqrt{W\psi^{\prime}(W)-1}\underset{W\to 0^{+}}{\propto}W^{-1/2}$, which
implies in particular that
$\frac{(W\psi^{\prime}(W)-1)}{\sqrt{W}}\underset{W\to
0^{+}}{\propto}W^{-3/2}$, and since we already proved in Theorem 3.3 that
$\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\underset{W\to
0^{+}}{\propto}1$, it follows by Proposition 3.2 that

$\int_{0}^{1}g_{4}(W)dW\propto\int_{0}^{1}W^{-3/2}\times 1\times
W^{n-1}\;dW<\infty.$

Moreover, due to [25, 28] we have
$\frac{\Gamma(nW)}{\Gamma(W)^{n}}\underset{W\to\infty}{\propto}n^{nW}W^{(n-1)/2}$
and $\sqrt{W\psi^{\prime}(W)-1}\underset{W\to\infty}{\propto}W^{-1/2}$, which
implies in particular that
$\frac{W\psi^{\prime}(W)-1}{\sqrt{W}}\underset{W\to\infty}{\propto}W^{-3/2}$,
and since we already proved in Theorem 3.3 that
$\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}=\exp\left(-nqW\right)n^{-nW}$,
where
$q=\log\left(\frac{\frac{1}{n}\sum_{i=1}^{n}x_{i}}{\sqrt[n]{\prod_{i=1}^{n}{x_{i}}}}\right)>0$,
by Proposition 3.2 it follows that

$\displaystyle\int_{1}^{\infty}g_{4}(W)dW$
$\displaystyle\propto\int_{1}^{\infty}W^{-3/2}\times\exp\left(-nqW\right)n^{-nW}\times
n^{nW}W^{(n-1)/2}dW$
$\displaystyle=\int_{1}^{\infty}W^{(n/2-1)-1}\exp\left(-nqW\right)\;dW=\frac{\Gamma(n/2-1)}{(nq)^{n/2-1}}<\infty,$

which concludes the proof. ∎

######  Theorem 3.10.

The posterior mean of $H$ relative to (22) is finite for all $n\geq 2$.

###### Proof.

Proceeding analogously as in the proof of Theorem 3.4 it follows that

$\displaystyle\left|E_{4}[H|x]\right|$
$\displaystyle\propto\int_{0}^{\infty}\left|\int_{-\infty}^{\infty}H\pi_{4}(H,W|\boldsymbol{x})\;dHdW\right|$
$\displaystyle\leq\int_{0}^{\infty}\delta_{2}(W)(W\psi^{\prime}(W)-1)\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}dW=\int_{0}^{1}g_{4}(W)dW+\int_{1}^{\infty}g_{4}(W)dW,$

where $\delta_{2}(W)$ is given as in the proof of Theorem 3.4 and

$\displaystyle
h_{4}(W)=\delta_{2}(W)(W\psi^{\prime}(W)-1)\frac{\left\\{\prod_{i=1}^{n}{x_{i}^{W}}\right\\}}{\left(\sum_{i=1}^{n}x_{i}\right)^{nW}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}.$

Since in the proof of Theorem 3.4 we showed that $\delta_{2}(W)\underset{W\to
0^{+}}{\propto}W^{-1}$, together with the proportionalities proved in Theorem
3.3 and Proposition 3.2 we have

$\int_{0}^{1}h_{4}(W)dW\propto\int_{0}^{1}W^{-1}\times W^{-3/2}\times 1\times
W^{n-1}\;dW<\infty.$

Moreover, letting $q>0$ as in the proof of Theorem 3.3, since we proved during
the proof of Theorem 3.4 that
$\delta_{2}(W)\underset{W\to\infty}{\propto}W(\log(W)+1)$ and since
$\log(W)+1\leq\exp(\log(W))=W$ for $W\geq 1$ it follows from Proposition 3.2
that

$\displaystyle\int_{1}^{\infty}h_{4}(W)dW$
$\displaystyle\propto\int_{1}^{\infty}W(\log(W)+1)\times
W^{-3/2}\times\exp\left(-nqW\right)n^{-nW}\times n^{nW}W^{(n-1)/2}dW$
$\displaystyle\leq\int_{1}^{\infty}W^{(n/2+1)-1}\exp\left(-nqW\right)\;dW=\frac{\Gamma(n/2+1)}{(nq)^{n/2+1}}<\infty,$

which concludes the proof. ∎

The marginal posterior distributions of $W$ is given by

$\pi_{4}(W|\boldsymbol{x})\propto\frac{W\psi^{\prime}(W)-1}{W^{\frac{1}{2}}}\frac{\Gamma(nW)}{\Gamma(W)^{n}}\left(\frac{\sqrt[n]{\prod_{i=1}^{n}{x_{i}}}}{\sum_{i=1}^{n}x_{i}}\right)^{nW}.$

Moreover, the conditional posterior distribution of $H$ is given by

$\displaystyle\pi_{4}(H|W,\boldsymbol{x})\propto\exp\left\\{-nWH-\delta(W,H)\sum_{i=1}^{n}x_{i}\right\\}.$

##  4 Simulation Study

A Monte Carlo simulation study is conducted to quantify and compare the
different non-informative priors’ impact on the entropy measure’s posterior
distribution. The Bias and Mean Square Error (MSE) were used to identify the
prior that provides the posterior distribution with posterior estimates closer
to the true value. These metrics are given by

$\operatorname{Bias}_{H}=\frac{1}{N}\sum_{i=1}^{N}(\hat{H}_{i}-H)\ \ \mbox{
and }\ \ \operatorname{MSE}_{H}\sum_{i=1}^{N}\frac{(\hat{H}_{i}-H)^{2}}{N},$
(23)

where $N=10,000$ is the number of samples used to estimate the MLE and
posterior quantities of interest. Here, we used the posterior mean as the
Bayes estimate due to its good properties. The estimates of $W$ are not
presented since we only considered $W$ as an auxiliary parameter to conduct
the Jacobian transformation, and therefore we are not interested in its
respective estimates.

In addition to the Bias and MSE, the coverage probabilities $CP$ were also
presented. Such metrics were obtained from the Bayesian credibility intervals
(CI) and the asymptotic confidence intervals of $H$. The nominal level assumed
was 0.95, i.e., we expect an adequate procedure to compute the
confidence/credibility intervals should return coverage probabilities closer
to 0.95. Regarding the Bias and MSE, the best approach among the selected ones
should return the Bias and MSE closest to zero.

The Newton-Raphson iterative method was used to maximize the likelihood in
order to obtain the MLE. For a fair comparison, the initial values used to
start the iterative procedures were the same values as the used to generate
the samples. In real applications, there is a need to set initial values. To
this end, we can use the closed-form maximum a posteriori estimator derived by
Louzada and Ramos [21] given by

$\tilde{\alpha}=\left(\frac{n-2.9}{n}\right)\cfrac{n\sum_{i=1}^{n}t_{i}}{\left(n\sum_{i=1}^{n}t_{i}\log\left(t_{i}\right)-\sum_{i=1}^{n}t_{i}\sum_{i=1}^{n}\log\left(t_{i}\right)\right)}$
(24)

and

$\tilde{\beta}=\frac{1}{n^{2}}\left(n\sum_{i=1}^{n}t_{i}\log\left(t_{i}\right)-\sum_{i=1}^{n}t_{i}\sum_{i=1}^{n}\log\left(t_{i}\right)\right).$
(25)

Therefore, the initial values for $H$ and $W$ are computed from
$\tilde{H}=\tilde{\alpha}-\log(\tilde{\beta})+\log\Gamma(\tilde{\alpha})+(1-\tilde{\alpha})\psi(\tilde{\alpha})$
and $\tilde{W}=\tilde{\alpha}$.

In the Bayesian framework, the posterior distribution’s marginal densities
involve double integrals to obtain the normalizing constants. Therefore, the
MCMC approach was adopted to obtain the posterior estimates. Moreover, the
Metropolis-Hastings algorithm was adopted to simulate quantities of interest
from the posterior densities. The first 500 samples were discarded in the
burn-in stage for each data simulated, and 5000 iterations were further
conducted. It was considered a thin parameter of 5 to avoid significant
autocorrelation among the sample, returning at the end,- 1000 simulated values
for each marginal distribution. The Geweke diagnostics [15] was considered to
confirm the convergence of chains under a confidence level of 95%. The
generated samples were used to estimate the posterior mean and the credibility
intervals, resulting in 10,000 estimates for $H$ and $W$.

The R software (R Core Development Team) was used for the simulation, where
the codes can be obtained on request from the corresponding author. For
$n=(20,\ldots,120)$, only the results sets at
$(\alpha,\beta)=\left(4,2)\right)$ and $(\alpha,\beta)=\left(2,0.5\right)$
were presented, which leaded respectively to $H=1.33$ and $H=2.27$. However,
the results were similar for different $\alpha$ and $\beta$ and therefore were
not presented here. For each sample from the posterior distribution, the
posterior mode and the credible intervals were evaluated for $\alpha$,
$\beta$, and $H$.

Table 1: The Bias(MSE) from the estimates of $\mu$ considering different values of $n$ with $N=100,000$ simulated samples, using the estimation methods: 1 - MLE, 2 - Jeffreys’s rule, 3 - Reference 1 prior, 4 - Reference 2 prior, and 5 - Tibshirani Prior. $\boldsymbol{\theta}$ | n | MLE | Jeffreys’s | Reference 1 | Reference 2 | Tibshirani
---|---|---|---|---|---|---
| 20 | 0.0509(0.0374) | 0.0374(0.0355) | 0.0359(0.0354) | 0.0251(0.0339) | 0.0114(0.0326)
| 30 | 0.0329(0.0234) | 0.0233(0.0225) | 0.0221(0.0224) | 0.0148(0.0218) | 0.0056(0.0213)
| 40 | 0.0260(0.0170) | 0.0185(0.0165) | 0.0176(0.0164) | 0.0121(0.0161) | 0.0052(0.0158)
| 50 | 0.0215(0.0139) | 0.0153(0.0136) | 0.0147(0.0135) | 0.0102(0.0133) | 0.0047(0.0131)
| 60 | 0.0154(0.0113) | 0.0101(0.0111) | 0.0097(0.0111) | 0.0059(0.0109) | 0.0014(0.0108)
$H=1.33$ | 70 | 0.0171(0.0099) | 0.0125(0.0097) | 0.0121(0.0097) | 0.0090(0.0095) | 0.0049(0.0094)
| 80 | 0.0118(0.0083) | 0.0078(0.0081) | 0.0074(0.0081) | 0.0047(0.0081) | 0.0012(0.0080)
| 90 | 0.0107(0.0073) | 0.0073(0.0072) | 0.0068(0.0072) | 0.0044(0.0071) | 0.0012(0.0071)
| 100 | 0.0090(0.0067) | 0.0058(0.0067) | 0.0054(0.0066) | 0.0032(0.0066) | 0.0005(0.0066)
| 110 | 0.0085(0.0061) | 0.0056(0.0061) | 0.0053(0.0061) | 0.0033(0.0060) | 0.0007(0.0060)
| 120 | 0.0084(0.0055) | 0.0056(0.0055) | 0.0054(0.0055) | 0.0036(0.0054) | 0.0012(0.0054)
| 20 | 0.0431(0.0290) | 0.0221(0.0272) | 0.0208(0.0271) | 0.0026(0.0261) | 0.0191(0.0261)
| 30 | 0.0353(0.0204) | 0.0205(0.0193) | 0.0199(0.0194) | 0.0077(0.0187) | 0.0065(0.0185)
| 40 | 0.0254(0.0152) | 0.0140(0.0147) | 0.0135(0.0146) | 0.0044(0.0143) | 0.0061(0.0142)
| 50 | 0.0220(0.0121) | 0.0127(0.0117) | 0.0122(0.0117) | 0.0050(0.0115) | 0.0033(0.0114)
| 60 | 0.0157(0.0097) | 0.0078(0.0095) | 0.0075(0.0095) | 0.0015(0.0093) | 0.0054(0.0093)
$H=2.27$ | 70 | 0.0141(0.0083) | 0.0073(0.0081) | 0.0070(0.0081) | 0.0019(0.0080) | 0.0040(0.0079)
| 80 | 0.0143(0.0072) | 0.0083(0.0070) | 0.0081(0.0070) | 0.0036(0.0069) | 0.0016(0.0069)
| 90 | 0.0117(0.0065) | 0.0064(0.0064) | 0.0061(0.0064) | 0.0022(0.0063) | 0.0024(0.0063)
| 100 | 0.0110(0.0056) | 0.0062(0.0055) | 0.0060(0.0055) | 0.0024(0.0055) | 0.0017(0.0055)
| 110 | 0.0079(0.0053) | 0.0035(0.0052) | 0.0033(0.0052) | 0.0000(0.0052) | 0.0037(0.0052)
| 120 | 0.0072(0.0048) | 0.0031(0.0047) | 0.0030(0.0047) | 0.0000(0.0047) | 0.0035(0.0047)
Table 2: The $CP_{95\%}$ from the estimates of $\mu$ and $\Omega$ considering different values of $n$ with $N=10,000,000$ simulated samples, using the estimation methods: 1 - MLE, 2- Jeffreys’s rule, 3 - Reference 1 prior, 4 - Reference 2 prior, and 5 - Tibshirani Prior. $\boldsymbol{\theta}$ | n | MLE | Jeffreys | Ref $W$ | Ref $H$ | Tibshirani
---|---|---|---|---|---|---
$H=1.33$ | 20 | 0.923 | 0.932 | 0.936 | 0.943 | 0.951
30 | 0.937 | 0.940 | 0.942 | 0.947 | 0.953
40 | 0.941 | 0.946 | 0.946 | 0.952 | 0.955
50 | 0.936 | 0.942 | 0.943 | 0.945 | 0.950
60 | 0.941 | 0.946 | 0.946 | 0.948 | 0.951
70 | 0.937 | 0.945 | 0.942 | 0.946 | 0.946
80 | 0.943 | 0.948 | 0.948 | 0.949 | 0.950
90 | 0.946 | 0.952 | 0.952 | 0.954 | 0.955
100 | 0.945 | 0.950 | 0.948 | 0.950 | 0.951
110 | 0.940 | 0.945 | 0.945 | 0.946 | 0.948
120 | 0.949 | 0.953 | 0.954 | 0.954 | 0.956
$H=2.27$ | 20 | 0.946 | 0.951 | 0.952 | 0.956 | 0.957
30 | 0.938 | 0.945 | 0.946 | 0.951 | 0.953
40 | 0.937 | 0.941 | 0.941 | 0.946 | 0.947
50 | 0.939 | 0.945 | 0.946 | 0.948 | 0.949
60 | 0.941 | 0.951 | 0.951 | 0.951 | 0.954
70 | 0.943 | 0.954 | 0.956 | 0.955 | 0.957
80 | 0.944 | 0.956 | 0.957 | 0.959 | 0.959
90 | 0.945 | 0.954 | 0.956 | 0.956 | 0.956
100 | 0.946 | 0.960 | 0.960 | 0.962 | 0.963
110 | 0.943 | 0.956 | 0.956 | 0.956 | 0.957
120 | 0.948 | 0.960 | 0.960 | 0.960 | 0.961

Tables 1 and 2 present the Bias, MSEs and $CP_{95\%}$ for the MLE and Bayesian
estimators of the entropy measure $H$. In particular, the results revealed
that:

  1. 1.

For all the parameter estimators, the Bias and MSE approach zero for large n,
which implies asymptotic unbiasedness, i.e., the Bias approaches zero, and the
MSE decreases as the number of samples increases.

  2. 2.

The posterior means using the reference prior 1 and 2 were superior to the
posterior means using Jeffreys prior and MLE. However, the posterior mean
using reference prior 2 was consistently superior to the posterior mean using
reference prior 1. The said performance is validated through the coverage
probability informed by the CI. Additionally, the coverage probability was
high for all the estimators, and the credibility of the interval increases
with sample size.

  3. 3.

For all estimators, the highest drop in Bias and MSE was observed when the
sample size increased from $20$ to $30$.

  4. 4.

Overall, the results show that the MLE performed worst, given its high bias
and MSE. On the other hand, the posterior estimates using the matching prior
provided Bayes estimates with smaller Bias and MSE; it was considered the most
adequate prior to estimating $H$.

According to the simulation results, the posterior distribution with
associated matching prior leads to the most precise results with the least
bias and MSE. The cited prior outperforms other objective priors and ML
estimates considered in this study, and therefore should be chosen as the most
appropriate prior for inference. Besides, the posterior estimates obtained
from the matching have superior theoretical properties, such as invariance
under one-to-one parameter transformation, consistent sampling, and
consistency under marginalization. Therefore, we conclude that the posterior
estimates derived from the matching prior distribution are more appropriate
and superior in making inferences of the gamma distribution’s population
parameter. To conduct the Bayesian analysis with the proposed Bayes estimator,
we have presented a function in R that can be used for this purpose, where the
details can be seen in Appendix A.

##  5 Application

###  5.1 Achaemenid dynasty

In this section, the proposed model is applied in the Achaemenid dynasty’s
rule time to quantify the variability in the Persian Empire’s political
institutions. The Achaemenid dynasty of the Achaemenid empire was the royal
house of the ancient Persians who ruled over Persia kingdom. The authority is
passed to the descendant of the same bloodline. The Persian Empire was build
and expanded through military conquest to extend political control to a
broader territory and people. They set out for wars to enlarge territories,
resulting in an imbalance in the political institutions. Despite the
governmental strategic techniques introduced in Cyrus the Great, the Persian
dynasty suffered several reoccurring internal political conflicts,
assassinations, and wars from internal and external entities, which shaped the
political institutions over the years.

Consequentially, the conflicts influenced each emperor’s tenor duration and
the pattern of a new emperor’s ascendancy. A stable government often has a
long time interval between successive emperors, whereas an unstable government
has a short tenor period. The induced patterns can be accounted for using
statistical tools such as entropy. The more frequent new emperors ascend the
throne, the higher the entropy and the more unstable the government. For the
sake of comparison, the proposed model was applied to the Roman Empire
timeline data, which was also analyzed by Ramos $et~{}al.$ [26].

Figure 1: Timeline (BC) containing the data, time series plot of the posterior
distribution of the entropy and autocorrelation plot for the same
distribution.

Figure 2 presents the timeline of the Achaemenid dynasty (top panel), the time
series of the Bayesian estimate of the entropy $H$ (down left panel), and the
autocorrelation plot (downright panel). The time series and the
autocorrelation plot indicate the chain’s convergence, which was also
confirmed by the Geweek test [15]. It is worth mentioning that the Kolmogorov-
Smirnov (KS) test was used (statistic D = 0.21) to confirm if the data can be
assumed to follow a gamma distribution.

Using the posterior distribution obtained from the matching prior, the Bayes
estimate of the Achaemenid dynasty’s entropy is $4.13$ with a $95\%$ credible
interval of $(3.55;4.73)$. Moreover, with the same prior, the posterior
estimate for the Roman Empire is $3.08$ with a $95\%$ credible interval of
$(2.80;3.36)$. Although the estimated entropies obtained for the Roman Empire
were high, the Achaemenid dynasty had a higher entropy, which implies that the
Achaemenid dynasty political institution is more volatile compared to the
Roman Empire. That is, the time between the successive emperor is
significantly different, shorter, and irregular for the Achaemenid dynasty,
which signifies instability in their political institutions. These results
support the historian’s claim that the Achaemenid Empire set out for wars and
consequently were exposed to internal and external insurgence, resulting in
instability and high entropy.

###  5.2 Harvest Sugarcane machine

Sugarcane farming is pertinent to Brazil’s economic growth and has heavily
contributed to its Gross Domestic Product (GDP). The production process
involves an automated harvesting mechanism, and the interest of the sugarcane
farmers is to sustain its harvesting mechanism for an extended period.
Moreover, the production chain must be kept in stable conditions to avoid
fluctuation in production and prevent wastage due to the deterioration of
sugarcane. This application modeled the entropy of the harvest machine failure
times in a Bayesian framework, using the Tibshirani prior, to quantify the
production process’s irregularity. A high entropy indicates a severe
irregularity in the production process. Otherwise, the production process is
steady.

The considered data was collected from January 2015 till August 2017, which
corresponds to 2.5 harvests and describe twenty-one failure times in days of
the suspension of a harvester sugarcane machine: 11, 19, 36, 4, 8, 11, 39, 74,
168, 27, 116, 3, 34, 1, 46, 12, 2, 56, 14, 52, 14.

Figure 2 presents the time series of the Bayesian estimation of the average
entropy (left panel) and the correlation lags (right panel). The convergence
of the estimate was tested using the Geweek test, and the KS test was used
(statistic D = 0.12) to confirm if the data can be assumed to follow a gamma
distribution.

Figure 2: Time series plot of the posterior distribution of the entropy and
autocorrelation plot for the same distribution.

The Bayes estimate of the entropy for the failure time of harvest sugarcane
machine is $4.55$ with a $95\%$ credible interval of $(4.04;5.09)$. The
estimated entropy is significantly high, which implies instability in the
production process. Moreover, the process is disrupted in an irregular time
duration, which causes the harvest output to be unpredictable and unreliable.
The harvesting machine must regularly pass a thorough maintenance check within
the harvester’s life circle to keep a steady production flow.

##  6 Final Remarks

The entropy concept was considered first in statistical thermodynamics, and it
was further modified to be applied in other fields. In information theory, the
Shannon entropy can measure the uncertainty of a random process. In
statistical inference, the Shannon entropy parameters are determined by the
maximum likelihood approach (MLE). However, using this approach, the results
are biased for small samples, and the confidence intervals may not achieve the
desired coverage probabilities as the asymptotic assumptions are not
fulfilled. In the present paper, we introduced a fully objective Bayesian
analysis to obtain the Shannon entropy’s posterior distribution to overcome
this limitation.

We considered objective priors where the obtained posterior distributions are
not overshadowed by prior information. The posterior distributions were
derived assuming the Jeffreys prior, reference priors, and matching priors,
all invariant under one-to-one transformations. Since the obtained priors are
improper, they could lead to improper posteriors, which is undesirable. We
proved that the obtained posteriors are proper distributions to address this
issue and thus can be used to conduct the Bayesian analysis. The posterior
mean was considered a Bayes estimator, and since they may not exist or not be
finite, we also proved that the posterior means are finite for any sample
sizes. Hence, four posterior distributions were proposed to conduct inference.
An intensive simulation study was conducted in order to select one Bayesian
estimator or the MLE. The posterior distribution using the matching prior
returned better results in terms of bias, mean square error, and coverage
probabilities compared with other methods, while the MLE returned the worst
results.

We analyzed the gamma distribution’s particular case, a more flexible and
general model than the exponential distribution, and has been applied to
describe many real phenomena. Although we considered the gamma distribution’s
particular case, our approach is general and can be further extended for any
probability distribution function.

The proposed Bayes estimator was implemented in an R language, available in
the appendix, to estimate the Shannon entropy measure. We apply the
implementation to estimate the entropy related to the Achaemenid dynasty’s
rule time, which returned a higher value compared with the Roman Empire. It
shows that the change in the throne in the Achaemenid dynasty was less
probable than the ones in the Roman Empire and indicates a significant
instability in their political institutions, which probably contributed to its
fall. Further, we analyzed the time until the failure of suspension of the
harvest sugarcane machine, where its entropy was estimated using the Bayesian
approach.

There are a large number of possible extensions of the current work. Other
distributions can be considered assuming the same context, and the Bayes
estimator of the Shannon entropy can be derived. Different entropy measure
types, such as Hartley, Rényi, and Tsallis entropy, could also be estimated
under a Bayesian approach. We plan to explore this line of research in the
future.

## Acknowledgements

Eduardo Ramos acknowledges financial support from São Paulo State Research
Foundation (FAPESP Proc. 2019/27636-9). Pedro L. Ramos acknowledges support
from São Paulo State Research Foundation (FAPESP Proc. 2017/25971-0).
Francisco Rodrigues acknowledges financial support from CNPq (grant number
309266/2019-0). Francisco Louzada is supported by the Brazilian agencies CNPq
(grant number 301976/2017-1) and FAPESP (grant number 2013/07375-0).

## References

  * [1] Milton Abramowitz and Irene A. Stegun.  Handbook of Mathematical Functions.  NBS, Washington, D.C., 10 edition, 1972.
  * [2] James O Berger and José M Bernardo.  Estimating a product of means: Bayesian analysis with reference priors.  Journal of the American Statistical Association, 84(405):200–207, 1989.
  * [3] James O Berger and José M Bernardo.  Ordered group reference priors with application to the multinomial problem.  Biometrika, 79(1):25–37, 1992.
  * [4] James O Berger and José M Bernardo.  Reference priors in a variance components problem.  In Bayesian Analysis in Statistics and Econometrics, pages 177–194. Springer, 1992.
  * [5] James O Berger, José M Bernardo, et al.  On the development of reference priors.  Bayesian statistics, 4(4):35–60, 1992.
  * [6] James O Berger, Jose M Bernardo, Dongchu Sun, et al.  Overall objective priors.  Bayesian Analysis, 10(1):189–221, 2015.
  * [7] Jose M Bernardo.  Reference posterior distributions for bayesian inference.  Journal of the Royal Statistical Society. Series B (Methodological), pages 113–147, 1979.
  * [8] José M Bernardo.  Reference analysis.  Handbook of statistics, 25:17–90, 2005.
  * [9] Youngseuk Cho, Hokeun Sun, and Kyeongjun Lee.  Estimating the entropy of a weibull distribution under generalized progressive hybrid censoring.  Entropy, 17(1):102–122, 2015.
  * [10] Pavel Cizek, Wolfgang Karl Härdle, and Rafał Weron.  Statistical tools for finance and insurance.  Springer Science & Business Media, 2005.
  * [11] Gauss M Cordeiro and Francisco Cribari-Neto.  An introduction to Bartlett correction and bias reduction.  Springer, 2014.
  * [12] Hérbetes de H Cordeiro and Gauss M Cordeiro.  Skewness for parameters in generalized linear models.  Communications in Statistics-Theory and Methods, 30(7):1317–1334, 2001.
  * [13] Yuge Du, Yu Guo, and Wenhao Gui.  Statistical inference for the information entropy of the log-logistic distribution under progressive type-i interval censoring schemes.  Symmetry, 10(10):445, 2018.
  * [14] OE Garcia.  Stochastic modeling of intermittent scrape-off layer plasma fluctuations.  Physical review letters, 108(26):265001, 2012.
  * [15] John Geweke.  Evaluating the accuracy of sampling-based approaches to the calculations of posterior moments.  Bayesian statistics, 4:641–649, 1992.
  * [16] Shanti S Gupta and Shanti S Gupta.  Gamma distribution in acceptance sampling based on life tests.  Journal of the American Statistical Association, 56(296):942–970, 1961.
  * [17] Gregory J Husak, Joel Michaelsen, and Chris Funk.  Use of the gamma distribution to represent monthly rainfall in africa for drought monitoring applications.  International Journal of Climatology: A Journal of the Royal Meteorological Society, 27(7):935–944, 2007.
  * [18] Harold Jeffreys.  An invariant form for the prior probability in estimation problems.  In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 186, pages 453–461. The Royal Society, 1946.
  * [19] Firuz Kamalov and Dmitry Denisov.  Gamma distribution-based sampling for imbalanced data.  Knowledge-Based Systems, 207:106368, 2020.
  * [20] Suchandan Kayal and Somesh Kumar.  Estimation of the shannon’s entropy of several shifted exponential populations.  Statistics & Probability Letters, 83(4):1127–1135, 2013.
  * [21] Francisco Louzada and Pedro Luiz Ramos.  Efficient closed-form maximum a posteriori estimators for the gamma distribution.  Journal of Statistical Computation and Simulation, 88(6):1134–1146, 2018.
  * [22] Helio S Migon, Dani Gamerman, and Francisco Louzada.  Statistical Inference: An Integrated Approach.  CRC, London, 2 edition, 2014.
  * [23] Robert B Miller.  Bayesian analysis of the two-parameter gamma distribution.  Technometrics, 22(1):65–69, 1980.
  * [24] Rahul Mukerjee and Dipak K Dey.  Frequentist validity of posterior quantiles in the presence of a nuisance parameter: higher order asymptotics.  Biometrika, 80(3):499–505, 1993.
  * [25] Pedro L Ramos, Jorge A Achcar, Fernando A Moala, Eduardo Ramos, and Francisco Louzada.  Bayesian analysis of the generalized gamma distribution using non-informative priors.  Statistics, 51(4):824–843, 2017.
  * [26] Pedro L. Ramos, Luciano da F. Costa, Francisco Louzada, and Francisco A. Rodrigues.  Power laws in the roman empire: a survival analysis.  2020\.
  * [27] Pedro Luiz Ramos, Dipak K. Dey, Francisco Louzada, and Eduardo Ramos.  On posterior properties of the two parameter gamma family of distributions.  Anais da Academia Brasileira de Ciências, 2020.
  * [28] Pedro Luiz Ramos, Francisco Louzada, and Eduardo Ramos.  Posterior properties of the nakagami-m distribution using noninformative priors and applications in reliability.  IEEE Transactions on reliability, 67(1):105–117, 2017.
  * [29] Mohammed K Shakhatreh, Sanku Dey, and MT Alodat.  Objective bayesian analysis for the differential entropy of the weibull distribution.  Applied Mathematical Modelling, 2020.
  * [30] Claude E Shannon.  A mathematical theory of communication.  The Bell system technical journal, 27(3):379–423, 1948.
  * [31] Dongchu Sun and Keying Ye.  Frequentist validity of posterior quantiles for a two-parameter exponential family.  Biometrika, 83(1):55–65, 1996.
  * [32] Robert Tibshirani.  Noninformative priors for one parameter of many.  Biometrika, 76(3):604–608, 1989.
  * [33] Jiao Yu, Wenhao Gui, and Yuqi Shan.  Statistical inference on the shannon entropy of inverse weibull distribution under the progressive first-failure censoring.  Entropy, 21(12):1209, 2019.

## Appendix A - Code of the Metropolis-Hasting algorithm within Gibbs

The code in R necessary to obtain the posterior mean of $H$ using the matching
prior is given by:

    library(coda)
    ###########################################################################
    ### Gibbs with Metropolis-Hasting algorithm ###
    ### R: Iteration Number; burn: Burn in; jump: Jump size; b= Control generation values ###
    ### posteriorW and posteriorH: logarithm of posterior densities ###
    ###########################################################################

    MCMC<-function(t,R,burn,jump,cW=1,seH=0.2) {
    posteriorW <- function (W) {
    p<-lgamma(n*W)-n*lgamma(W) +0.5*log(W*trigamma(W)-1)-0.5*log(W)+(W)*sum(log(t))-(n*W)*log(sum(t))
    return(p) }
    posteriorH <- function (H,W) {
    delta<-exp(W+lgamma(W)+(1-W)*digamma(W)-H)
    p<-n*W*log(delta)-delta*sum(t)
    return(p) }
    n<-length(t)
    vH<-length(R+1); vW<-length(R+1)
    beta<-(1/(n*(n-1)))*(n*sum((t)*log(t))-sum(t)*sum(log(t)))
    vW[1]<-(n-2.9)/(((n*sum((t)*log(t))/sum(t))-sum(log(t))))
    vH[1]<-vW[1]-log(beta)+lgamma(vW[1])+(1-vW[1])*digamma(vW[1])
    c1<-rep(0,times=R) ; i<-1
    while (i<=R) {
    prop1<-rgamma(1,shape=cW*vW[i],rate=cW)
    d1<-dgamma(vW[i],shape=cW*prop1,rate=cW,log=TRUE)-dgamma(prop1,shape=cW*vW[i],rate=cW,log=TRUE)
    ratio1<-posteriorW(prop1)-posteriorW(vW[i])+d1
    has<-min(1,exp(ratio1)); u1<-runif(1)
    if (is.finite(has)) {
    if (u1<has) {vW[i+1]<-prop1} else {vW[i+1]<-vW[i]}} else {vW[i+1]<-vW[i]}
    prop2<- rnorm(1,vH[i],seH)
    d2<-dnorm(vH[i], mean=prop2, sd=seH, log = TRUE)-dnorm(prop2, mean=vH[i], sd=seH, log = TRUE)
    ratio2<-posteriorH(prop2,vW[i+1])-posteriorH(vH[i],vW[i+1]) +d2
    has<-min(1,exp(ratio2)); u2<-runif(1)
    if (u2<has & is.double(has)) {vH[i+1]<-prop2 ; c1[i]<-0} else {vH[i+1]<-vH[i] ; c1[i]<-1}
    i<-i+1 }
    vvH<- vH[seq(burn,R,jump)]; H<-mean(vvH)
    ge1<-abs(geweke.diag(vvH)$z[1])
    prli<-quantile(vvH, probs = 0.025, na.rm = FALSE,names = FALSE,type = 7)
    prls<-quantile(vvH, probs = 0.975, na.rm = FALSE,names = FALSE,type = 7)
    return(list(acep=(1-sum(c1)/length(c1)),H=H, LCI_H=prli, UCI_H=prls, Geweke.statistics=ge1))
    }

    ################################################################
    ## Example ###
    palpha<-2; pbeta=0.5; n=50 #Parameters
    t<-rgamma(n, palpha, rate=pbeta) ### data vector

    #Calling the function
    MCMC(t,R=2000,burn=500,jump=5)

    #Output in R
    $acep          ##Aceptance rate
    [1] 0.4795

    $H        ##Posterior mean of H
    [1] 4.880441

    $LCI_H    ##Lower credibility interval of H
    [1] 1.880321

    $UCI_H   ##Upper credibility interval of H
    [1] 2.32128

    $Geweke.statistics ## Geweke Statitics
    1.054945