# Model selection of stochastic simulation algorithm based on generalized
divergence measures

Papa NGOM and B. Don Bosco DIATTA LMA-Laboratoire de Mathématiques Appliquées
Université Cheikh Anta Diop BP 5005 Dakar-Fann Sénégal
e-mail<EMAIL_ADDRESS>
###### Abstract

MCMC methods (Monte Carlo Markov Chain) are a class of methods used to perform
simulations per a probability distribution $P$. These methods are often used
when we have difficulties to directly sample per a given probability
distribution $P$ . This distribution is then considered as a target and
generates a Markov chain $(X_{n})_{n\in\mathbb{N}}$ that, when $n$ is large we
have $X_{n}\sim P$. These MCMC methods consist of several simulation
strategies including the _Independent Sampler (IS)_ , the _Random Walk of
Metropolis Hastings(RWMH)_, the _Gibbs sampler_ , the _Adaptive Metropolis
(AM)_ and _Metropolis Within Gibbs (MWG)_ strategy. Each of these strategies
can generate a Markov chain and is associated with a convergence speed. It is
interesting, with a given target law, to compare several simulation strategies
for determining the best. Chauveau and Vandekerkhove [4] have compared IS and
RWMH strategies using the Kullback-Leibler divergence measure. In our article
we will compare our five simulation methods already mentioned using
generalized divergence measures. These divergence measures are taken in family
of $\alpha$-divergence measures [6] , with a parameter $\alpha$. This is the
Rényi divergence, Tsallis divergence and $D_{\alpha}$ divergence .

###### keywords:

MCMC methods, Metropolis-Hastings algorithm, Gibbs sampler, Adaptive
Metropolis, Metropolis within Gibbs, simulation strategy, target density,
proposal density, $\alpha$-divergence.
AMS Subject Classification : 60J60, 62F03, 62F05, 94A17.

††thanks:  This research was supported, in part, by grants from Université
Cheikh Anta Diop, Dakar (Sénégal). We are grateful to many seminar
participants and to anonymous referees for comments.

##  1 Introduction

In many areas of science, computer is essential. So computers are used for
testing, calculations, simulations … Beyond performance, the computation time
is very important in the choice of calculation methods. In the field of
statistical, the simulation of random variables are required to perform
assessments such integral calculation. Thus we resort to numerical analysis
methods or probabilistic methods. It is the latter that concern us in
particular. The integral form is rewritten like an expectation $\mathbb{E}(X)$
and use the _law of large numbers_ for the Monte Carlo integration and the
_ergodic theorem_ for MCMC methods. We will focus here on MCMC methods with
which one performs sampling according to density $f$ via Markov chains. Among
MCMC methods we have Metropolis-Hastings algorithms, Gibbs sampler and others
adaptive and hybrid algorithms.

Thus the major challenge of these algorithms, like many computational
techniques, is the saving time i.e. their speedy convergence to the stationary
distribution with density $f$. This is due to the fact that the convergence
time is a brand of efficiency. This paper investigates the problem of
selecting a good strategy simulation. We have several simulation strategies
that are _Independence Sampler_ , _Random Walk of Metropolis-Hastings_ [10],
_Gibbs sampler_ , _Metropolis Within Gibbs_ [8], _Adaptive Metropolis_ [7]…The
problem for all these strategies is the time they take to the generated a
Markov chain which converge to the stationary distribution with probability
density function $f$.

The objective of this paper is to compare different strategies for simulation
in order to select the best. Chauveau and Vandekerkhove [4] have given in
their paper a methodology for selecting the best of $k$ candidate strategies
with the Kullback-Leibler divergence. We use in our article generalized
divergence measures that are $\alpha$-divergence [11] to compare simulation
strategies. These divergence measures that we use are _$\alpha$ -divergence_
$D_{\alpha}$, the _Rényi $\alpha$-divergence_ $R_{\alpha}$, the _Tsallis
$\alpha$-divergence_ $T_{\alpha}$. We know that these divergences are all
generalizations of the Kullback-Leibler divergence. We have the Kullback-
Leibler divergence when the parameter $\alpha$ of these divergence measures
tends to 1. If we have two simulation strategies $s_{1}$ and $s_{2}$ of
respective probability densities $p^{n}_{1}$111$p_{1}^{n}$ is density function
of $X_{n}$ of Markov chain generated by strategy $s_{1}$ and $p_{2}^{n}$ at
time $n$ and a target density $f$ which we want to have samples we can use
these divergence measures to compare $s_{1}$ and $s_{2}$. Indeed, if we
consider for example the Rényi divergence, we can compare
$R_{\alpha}(p_{1}^{n},f)$ and $R_{\alpha}(p_{2}^{n},f)$, for each iteration
$n$, to see which of $s_{1}$ and $s_{2}$ is the best.

We will show in Section 2 these three divergence measures. We know that they
belong to the family of Csiszár $\phi$-divergences. In Section 3 we will show
two important results that are convergence to 0 of these divergences and
building estimator for each divergence. We will show in Section 4 the
asymptotic distributions of our estimators. Then we have some application
examples in order to illustrate our methodology (Section 5). Finally, we do
discussion in Section 6 to explain our various examples and see the lesson we
can draw threreof.

##  2 Description of divergence measures

Let $(\cal X,\cal A,\lambda)$ be an arbitrary measure space with $\lambda$
being a finite or $\sigma-$finite measure. Let also $\mu_{1},\mu_{2}$
probability measures on $\cal X$ such that $\mu_{1},\mu_{2}\ll\lambda$
(absolutely continuous).
Denote the Randon-Nikodym derivatives (densities) of $\mu_{i}$ with respect to
$\lambda$ by $p_{i}(x)$:

$p_{i}(x)=\frac{\mu_{i}(dx)}{\lambda(dx)},\quad i=1,2.$

###### Definition 2.1

Kullback-Leibler’s relative divergence (also called relative entropy) between
two probability measures $\mu_{1}$ , $\mu_{2}$ is defined by

$K(\mu_{1},\mu_{2})=\int_{\cal
X}p_{1}(x)\log\big{(}\frac{p_{1}(x)}{p_{2}(x)}\big{)}\lambda(dx)=\mathbb{E}_{\mu_{1}}\left[\log\frac{p_{1}(x)}{p_{2}(x)}\right]$
(2.1)

We can also write $K(p_{1},p_{2})$.

We will focus on in this article to a particular family of _Csiszár
$\phi$-divergence_ that is the family of _$\alpha$ -divergence_ measures. Thus
we have the $\alpha$-divergence [1], the Rényi $\alpha$-divergence and the
Tsallis $\alpha$-divergence of the discrepancy which we will work.

######  Definition 2.2 ($\alpha$-divergence)

The $\alpha$-divergence is defined by

$D_{\alpha}(\mu_{1},\mu_{2})=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\int_{\mathcal{X}}p_{1}^{\alpha}(x)p_{2}^{1-\alpha}(x)\lambda(\mathrm{d}x)\Big{)},\,\alpha>0\
\hbox{and}\ \alpha\neq 1.$ (2.2)

######  Definition 2.3 (Rényi $\alpha$-divergence)

Rényi (1961) for the first time gave one generalization of the relative
entropy given in 2.1. It is defined by

$R_{\alpha}(\mu_{1},\mu_{2})=\frac{1}{\alpha-1}\log\int_{\cal
X}p_{1}^{\alpha}(x)p_{2}^{1-\alpha}(x)\lambda(dx),\,\alpha>0\ \hbox{and}\
\alpha\neq 1.$ (2.3)

######  Definition 2.4 (Tsallis $\alpha$-divergence)

Tsallis $\alpha$-divergence is defined by

$T_{\alpha}(\mu_{1},\mu_{2})=\frac{1}{\alpha-1}\Big{(}\int_{\mathcal{X}}p_{1}^{\alpha}(x)p_{2}^{1-\alpha}(x)\lambda(\mathrm{d}x)-1\Big{)},\,\alpha>0\
\hbox{and}\ \alpha\neq 1$ (2.4)

Chauveau and Vandekerkhove [4] studied the strategies of simulation with
Kullback-Leibler divergence. Their study was the comparison of simulation
strategies by the Kullback-Leibler divergence which allowed to choose the best
among candidate strategies. They used the Metropolis-Hastings (MH) algorithm
to explain their method. Thus, with a target density f, proposal density q and
an initial distribution $p_{i}^{0}$ corresponding to the strategy i of MH
algorithm they assumed several conditions on these densities to provide
certain regularity properties, as the Lipschitz property, on the successive
densities $p_{i}^{n}$ of strategy i. This allowed to obtain an consistent
estimate of the entropy

$\mathcal{H}(p_{i}^{n})=\int_{\mathcal{X}}p_{i}^{n}(x)\log(p_{i}^{n}(x))\lambda(\mathrm{d}x)$

of $p_{i}^{n}$ at time $n$.

In this present paper we propose to study how to find the optimal stochastic
simulation algorithm may come from the MH methods (IS and RWMH), Gibbs
sampler, or recent adaptive (AM) and hybrid (MWG) methods. For this, we will
use our three divergence measures. We see that we can find the Kullback-
Leibler divergence if for each divergence, the parameter $\alpha$ tends to 1.
The interest of these divergences is that its subtracts the target density
$f$, the proposal density $q$ and initial density $p_{i}^{0}$ many assumptions
like in Chauveau and Vandekerkhove [5].

##  3 Convergence and estimators of divergence measures

###  3.1 Convergence of $\alpha$-divergence measures

We will use a result of Holden (1998) to show the geometric convergence of the
MH algorithm under a minoration condition: if there exists $a\in(0,1)$ such
that $q(y|x)\geq af(y)$ for all $x,y\in\mathcal{X}$, then

$\forall
y\in\mathcal{X},\Big{|}\frac{p^{n}(y)}{f(y)}-1\Big{|}\leq(1-a)^{n}\Big{\|}\frac{p^{0}}{f}-1\Big{\|}_{\infty}$
(3.5)

###### Theorem 3.1

If the proposal density of the Metropolis-Hastings algorithm satisfies
$q(y|x)\geq\delta f(y)$, for all $x,y\in\mathcal{X}$, and $\delta\in(0,1)$,
then with $\alpha\in(1,+\infty)$

$D_{\alpha}(p^{n},f)\leq\frac{1}{\alpha(1-\alpha)}(1-(r\nu^{n}+1)^{\alpha}),$
(3.6) $R_{\alpha}(p^{n},f)\leq\frac{\alpha}{\alpha-1}r\nu^{n}$ (3.7)
$T_{\alpha}(p^{n},f)\leq\frac{1}{\alpha-1}\big{(}(r\nu^{n}+1)^{\alpha}-1\big{)}$
(3.8)

where $r=\|\frac{p^{0}}{f}-1\|_{\infty}>0,$ and $\nu=(1-\delta).$

Proof.

a) $\alpha$-divergence $D_{\alpha}$

$D_{\alpha}(p^{n},f)=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\int_{\mathcal{X}}\big{(}\frac{p^{n}(y)}{f(y)}-1+1\big{)}^{\alpha}f(y)\lambda(\mathrm{d}y)\Big{)}$
(3.9)

with $\alpha\in]1,+\infty[$, we have

$1-\int_{\mathcal{X}}\big{(}\frac{p^{n}(y)}{f(y)}-1+1\big{)}^{\alpha}f(y)\lambda(\mathrm{d}y)\geq
1-\int_{\mathcal{X}}\big{(}\big{|}\frac{p^{n}(y)}{f(y)}-1\big{|}+1\big{)}^{\alpha}f(y)\lambda(\mathrm{d}y)$

Therefore using (3.5) we have

$\displaystyle D_{\alpha}(p^{n},f)$ $\displaystyle\leq$
$\displaystyle\frac{1}{\alpha(1-\alpha)}\Big{(}1-\int_{\mathcal{X}}\big{(}\big{|}\frac{p^{n}(y)}{f(y)}-1\big{|}+1\big{)}^{\alpha}f(y)\lambda(\mathrm{d}y)\Big{)}$
(3.10) $\displaystyle\leq$
$\displaystyle\frac{1}{\alpha(1-\alpha)}\Big{(}1-\int_{\mathcal{X}}(r\nu^{n}+1)^{\alpha}f(y)\lambda(\mathrm{d}y)\Big{)}$
$\displaystyle\leq$
$\displaystyle\frac{1}{\alpha(1-\alpha)}\big{(}1-(r\nu^{n}+1)^{\alpha}\big{)}$

b) Rényi $\alpha$-divergence

We will use (3.5) and the same procedure give us

$\displaystyle R_{\alpha}(p^{n},f)$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha-1}\log\int_{\mathcal{X}}\big{(}\frac{p^{n}(y)}{f(y)}\big{)}^{\alpha}f(y)\lambda(\mathrm{d}y)$
(3.11) $\displaystyle\leq$
$\displaystyle\frac{\alpha}{\alpha-1}\log(r\nu^{n}+1)$ $\displaystyle\leq$
$\displaystyle\frac{\alpha}{\alpha-1}r\nu^{n}$

c) Tsallis $\alpha$-divergence

Here too we have

$\displaystyle T_{\alpha}(p^{n},f)$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha-1}\Big{(}\int_{\mathcal{X}}\big{(}\frac{p^{n}(y)}{f(y)}\big{)}^{\alpha}f(y)\lambda(\mathrm{d}y)-1\Big{)}$
(3.12) $\displaystyle\leq$
$\displaystyle\frac{1}{\alpha-1}\big{(}\big{(}r\nu^{n}+1\big{)}^{\alpha}-1)\big{)}$

$\square$

These three divergence measures are all positive, we see that those converge
to zero under the conditions of the Theorem 3.1. This allows us to see that
the densities $p^{n}$ of random variables $X_{n}$ of the Markov chain converge
to the stationary distribution $f$ when $n$ goes to infinity.

###  3.2 Estimators for the three divergence measures

We first make an estimator of $\alpha$-divergence $D_{\alpha}$, in the same
way we can have an estimator of divergences $T_{\alpha}$ and Rényi divergence
$R_{\alpha}$. Suppose $\lambda$ is the Lebesgue measure on
$\mathcal{X}=\mathbb{R}^{d}$. Then we have

$D_{\alpha}(p^{n},f)=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\int\Big{(}\frac{f(x)}{p^{n}(x)}\Big{)}^{1-\alpha}p^{n}(x)\mathrm{d}x\Big{)}$

We can initially think to write the integral like a mathematical expectancy
and apply the method of Monte Carlo integration. We would have then

$D_{\alpha}(p^{n},f)=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\mathbb{E}\Big{(}\big{(}\frac{f(X)}{p^{n}(X)}\big{)}^{1-\alpha}\Big{)}\Big{)}$

Using the strong law of large numbers we will have the following estimator

$\hat{D}_{\alpha}(p^{n},f)=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\frac{1}{N}\sum_{i=1}^{N}\big{(}\frac{f(X_{i})}{p^{n}(X_{i})}\big{)}^{1-\alpha}\Big{)}$

which converges almost surely to $D_{\alpha}(p^{n},f)$.

But very often encountered in practice, cases where the density $f$ and the
densities $p^{n}$ are not analytically known. There are also cases where $f$
is known to a multiplicative constant: the Bayesian context.

It is often difficult to obtain the analytical form of $p^{n}$. If $f$ is
known analytically, it is not useful to estimate it. But let us consider the
more general case where $p^{n}$ and $f$ are not analytically known . Then we
use the methods of nonparametric estimation of probability densities. Póczos
and Schneider [3] then proposed in their paper an estimator of these densities
based on _k-NN method_(k nearest neighbor).

####  3.2.1 k-NN based density estimator

This is a method of density estimation proposed by Loftsgaarden and
Quesenberry [9]. It has a set of data training
$E=\left\\{x_{1},x_{2},\ldots,x_{N}\right\\}$ where $x_{i}$ are independent
observations with common probability density function $g$. For a data $z$, the
problem is the estimation of a probability density $g$ at $z$. The principle
is to find among $x_{i}$ the kth nearest neighbor of $z$ considering euclidean
distance and we are interested in the distance between the two points
($x_{i,k}$ and $z$). If $x_{i,k}$ and $z$ are points of $\mathbb{R}^{d}$ this
distance $d(z,x_{i,k})$ allows us to calculate the volume of the hypersphere
with center $z$ and radius $r=d(z,x_{i,k})$.

The proposed estimator is as follows

$\hat{g}_{k,N}(z)=\frac{k/N}{V(Hs(z,d(z,x_{i,k})))}=\frac{k/N}{\pi^{d/2}r^{d}/\Gamma((d/2)+1)}$

where $V(Hs(z,d(z,x_{i,k})))=\pi^{d/2}r^{d}/\Gamma((d/2)+1)$ is the volume of
d-dimensional ball around $z\in\mathbb{R}^{d}$ with radius $r>0$, $\Gamma(.)$
is the Gamma function. Loftsgaarden and Quesenberry (1965) suggested after
several experiments to take $k$ equals to $\sqrt{N}$. We will round it if
necessary to the nearest whole.

We find in [11] theorems of convergence for density’s estimators (k-NN
estimators), showing the consistency of these estimators.

######  Theorem 3.2 (k-NN density estimators, convergence in probability)

Let $k(n)$ an integer that denotes the $k(n)$th nearest neighbor in the sample
of size n. If $\lim_{n\rightarrow\infty}k(n)=\infty$, and
$\lim_{n\rightarrow\infty}n/k(n)=\infty$, then
$\hat{g}_{k(n),N}(x)\longrightarrow_{p}g(x)$ for almost all $x$.

######  Theorem 3.3 (Almost sure convergence in sup norm)

If
$\lim_{n\rightarrow\infty}k(n)/\log(n)=\infty$ and
$\lim_{n\rightarrow\infty}n/k(n)=\infty$, then

$\lim_{n\rightarrow\infty}\sup_{x}|\hat{g}_{k(n),N}(x)-g(x)|=0,\quad\textrm{almost
surely.}$

a) Application to $p^{n}$

Let $X_{1:N}=(X_{1},...,X_{N})$ be a simulated sample according to the law
with density $p^{n}$ with $X_{i}$ i.i.d. We choose one $X_{i}$ on the sample.
Let $\rho_{k,i}$ the euclidean distance between $X_{i}$ and its kth nearest
neighbor on the sample.

The estimate of $p^{n}$ at $X_{i}$ is

$\hat{p}_{k,n,N}(X_{i})=\frac{k/(N-1)}{V(Hs(X_{i},\rho_{k,i}))}$ (3.13)

$V(Hs(X_{i},\rho_{k,i}))=c\,\rho_{k,i}^{d}$, where $c$ stands for the volume
of a d-dimensional unit ball.

Hence

$\hat{p}_{k,n,N}(X_{i})=\frac{k}{(N-1)c\,\rho_{k,i}^{d}}$ (3.14)

b) Application to $f$

We should also have the i.i.d. simulations according to the density $f$. But
the principle in this article is that we don’t know how to do i.i.d.
simulations directly from $f$. Then we apply the MCMC simulation algorithms.
After several iterations (large $n$) we can begin to recover samples. However,
we take the precaution that two successive samples are separated by $n_{0}$
iterations. We get now the samples $Y_{1:M}=(Y_{1},...,Y_{M})$ of random
variables which are ’independent’ with density $f$. Let $\gamma_{k,i}$ the
euclidean distance between $X_{i}$ and its kth nearest neighbor on the sample
$Y_{1:M}$.

The estimator of $f$ at $X_{i}$ is

$\hat{f}_{k,N}(X_{i})=\frac{k}{M\,c\,\gamma_{k,i}^{d}}$ (3.15)

####  3.2.2 Construction of estimators

Póczos and Schneider (2011) proposed an estimator of the integral part that is
common to the three divergence measures that we have presented. However, their
first estimator is asymptotically biased. They showed that multiplying by a
constant, they obtained an asymptotically unbiased estimator. Starting from it
we construct our estimators for the $\alpha$-divergence $D_{\alpha}$, the
Rényi $\alpha$-divergence $\mathcal{R}_{\alpha}$ and the Tsallis
$\alpha$-divergence $T_{\alpha}$.

a) Estimator for $\alpha$-divergence $D_{\alpha}$

Begin by noting
$M_{\alpha}(p^{n},f)=\int(p^{n}(x))^{\alpha}(f(x))^{1-\alpha}\mathrm{d}x$
common integral part of the three divergence measures.

Once the densities $p^{n}$ and $f$ estimated at $X_{i}$ , we have the
following estimator of $\hat{D}_{\alpha,N}(p^{n},f)$.

$\displaystyle\hat{D}_{\alpha,N}(p^{n},f)$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha(1-\alpha)}\Big{(}1-\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{\hat{f}_{k,N}(X_{i})}{\hat{p}_{k,n,N}(X_{i})}\Big{)}^{1-\alpha}\Big{)}$
(3.16) $\displaystyle=$
$\displaystyle\frac{1}{\alpha(1-\alpha)}\Big{(}1-\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{(N-1)\rho_{k,i}^{d}}{M\gamma_{k,i}^{d}}\Big{)}^{1-\alpha}\Big{)}$

We have now

$\hat{M}_{\alpha,N}(p^{n},f)=\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{(N-1)\rho_{k,i}^{d}}{M\gamma_{k,i}^{d}}\Big{)}^{1-\alpha}$

that is an asymptotically biased estimator of $M_{\alpha}(p^{n},f)$ (Póczos
and Schneider (2011)). When we multiply it by a constant

$B_{k,\alpha}=\frac{(\Gamma(k))^{2}}{\Gamma(k-\alpha+1)\Gamma(k+\alpha-1)}$
(3.17)

independent of $p^{n}$ and $f$, with $\Gamma(.)$ the Gamma function. We obtain
an asymptotically unbiased estimator.

We obtain now a new asymptotically unbiased estimator of $D_{\alpha}(p^{n},f)$

$\hat{D}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{(N-1)\rho_{k,i}^{d}}{M\gamma_{k,i}^{d}}\Big{)}^{1-\alpha}\times
B_{k,\alpha}\Big{)}$ (3.18)

b) Estimator for Tsallis $\alpha$-divergence

Applying the same process we obtain an asymptotically unbiaised estimator
$\hat{T}_{\alpha,N,k}(p^{n},f)$ of $T_{\alpha}(p^{n},f)$

$\hat{T}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha-1}\Big{(}\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{(N-1)\rho_{k,i}^{d}}{M\gamma_{k,i}^{d}}\Big{)}^{1-\alpha}\times
B_{k,\alpha}-1\Big{)}$ (3.19)

c) Estimator for Rényi divergence

The Rényi measure divergence has a characteristic due to the presence of the
logarithmic function that takes $M_{\alpha}(p^{n},f)$ as argument. We have not
shown in this paper that its estimator is unbiased or asymptotically unbiased.
However, we propose an estimator

$\hat{R}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha-1}\log\Big{(}\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{(N-1)\rho_{k,i}^{d}}{M\gamma_{k,i}^{d}}\Big{)}^{1-\alpha}\times
B_{k,\alpha}\Big{)}$ (3.20)

Now give some theorems which show that
$\hat{M}_{\alpha,N}(p^{n},f)B_{k,\alpha}$ is asymptotically unbiased (Póczos
and Schneider [11]). Define first the following function

$F(x,p,\delta,\omega)=\sum_{j=1}^{k-1}\Big{(}\frac{1}{j!}\Big{)}^{\omega}\Gamma(\gamma+j\omega)\Big{(}\frac{p(x)+\delta}{p(x)-\delta}\Big{)}^{j\omega}\big{(}p(x)-\delta\big{)}^{-\gamma}\big{(}(1-\delta)\omega\big{)}^{-\gamma-j\omega}$
(3.21)

######  Theorem 3.4 (Asymptotic unbiasedness)

Assume that
$0<\gamma=1-\alpha<k$, $p$ is bounded away from 0, $p$ is uniformly Lebesgue
approximable, $\exists\delta_{0}$ such that
$\forall\delta\in(0,\delta_{0}),\>\int F(x,p,\delta,1)p(x)\mathrm{d}x<\infty$,
$\int\|x-y\|^{\gamma}p(y)\mathrm{d}y<\infty$ for almost all
$x\in\mathbb{R}^{d}$,
$\int\int\|x-y\|^{\gamma}p(y)p(x)\mathrm{d}y\mathrm{d}x<\infty$ and that $q$
is bounded from above. Then

$\lim_{n,m\rightarrow\infty}\mathbb{E}\Big{(}\hat{M}_{\alpha,N}(p^{n},f)B_{k,\alpha}\Big{)}=M_{\alpha}(p^{n},f)$
(3.22)

i.e., the estimator is asymptotically unbiased.

The following theorem provide conditions under which estimator
$\hat{M}_{\alpha,N}(p^{n},f)$ is L2 consistent (Póczos and Schneider (2011)).

######  Theorem 3.5 ($L_{2}$ consistency)

We have the following assumptions: $k\geq 2$, $0<\gamma=1-\alpha<(k-1)/2$, $p$
is bounded away from 0, $p$ is uniformly Lebesgue approximable,
$\exists\,\delta_{0}$ such that $\forall\,\delta\,\in\,(0,\delta_{0})\>\int
F(x,p,\delta,1/2)p(x)\mathrm{d}x<\infty$,
$\int\|x-y\|^{\gamma}p(y)\mathrm{d}y<\infty$ for almost all
$x\in\mathbb{R}^{d}$,
$\int\int\|x-y\|^{\gamma}p(y)p(x)\mathrm{d}y\mathrm{d}x<\infty$, and that $q$
is bounded above. Then

$\lim_{n,m\rightarrow\infty}\mathbb{E}\Big{(}\big{(}\hat{M}_{\alpha,N}(p^{n},f)B_{k,\alpha}-M_{\alpha}(p^{n},f)\big{)}^{2}\Big{)}=0$
(3.23)

that is, the estimator is $L_{2}$ consistent.

The last theorem show the consistency $\hat{M}_{\alpha,N}(p^{n},f)$ which is
the estimator of the common integral part to the three differences. Based on
this, estimators $\hat{\mathcal{R}}_{\alpha,N,k}(p^{n},f)$,
$\hat{D}_{\alpha,N,k}(p^{n},f)$ and $\hat{T}_{\alpha,N,k}(p^{n},f)$ are also
consistent.

##  4 Asymptotic distribution of divergences estimators

We seek to know the asymptotic distribution of our estimators. The asymptotic
distribution of these estimators is studied under the assumption of continuity
of densities $p^{n}$ and $f$. If densities $p^{n}$ and $f$ are continuous,
their estimators will also be continuous.

###  4.1 Asymptotic distribution of $\alpha$-divergence estimator

Recall that the asymptotically unbiased estimator of

$D_{\alpha}(p^{n},f)=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\int\Big{(}\frac{f(x)}{p^{n}(x)}\Big{)}^{1-\alpha}p^{n}(x)\mathrm{d}x\Big{)}$

is

$\hat{D}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{(N-1)V(Hs(X_{i},\rho_{k,i}))}{mV(Hs(X_{i},\gamma_{k,i}))}\Big{)}^{1-\alpha}\times
B_{k,\alpha}\Big{)}.$

Let

$h(X_{i})=\Big{(}\frac{(N-1)V(hs(X_{i},\rho_{k,i}))}{mV(Hs(X_{i},\gamma_{k,i}))}\Big{)}^{1-\alpha}$

and

$\bar{Y}_{N}=\frac{1}{N}\sum_{i=1}^{N}h(X_{i})$

We know that $X_{i}$ are i.i.d. $\sim p^{n}$ . Function $h$ which is the ratio
of two continuous functions is also continuous. We have therefore the
independence of $h(X_{i}),\,i=1,\ldots,N$. We can then apply the central limit
theorem. Then we will have

$\displaystyle\frac{\bar{Y}_{N}-\mathbb{E}(\bar{Y}_{N})}{\sqrt{Var(\bar{Y}_{N})}}$
$\displaystyle\xrightarrow{d}$ $\displaystyle\mathcal{N}\big{(}0,1\big{)}$
$\displaystyle\bar{Y}_{N}\times
B_{k,\alpha}-B_{k,\alpha}\mathbb{E}(\bar{Y}_{N})$
$\displaystyle\xrightarrow{d}$
$\displaystyle\mathcal{N}\Big{(}0,\sigma^{2}B^{2}_{k,\alpha}\Big{)}$

with $\sigma^{2}=\lim_{N\to\infty}Var(\bar{Y}_{N})$ we have a normal
distribution. But $\bar{Y}_{N}$ is asymptotically biased (Póczos and Schneider
[11]). Multiplyind it by a constant $B_{k,\alpha}$ as stated by these authors
gives us asymptotically unbiased estimator of
$M_{\alpha}(p^{n},f)=\int(p^{n}(x))^{\alpha}(f(x))^{1-\alpha}\mathrm{d}x$. We
will have

$\bar{Y}_{N}B_{k,\alpha}\xrightarrow{d}\mathcal{N}\big{(}M_{\alpha}(p^{n},f)\,,\sigma^{2}B_{k,\alpha}^{2}\big{)}$

Hence

$\hat{D}_{\alpha,N,k}(p^{n},f)\xrightarrow{d}\mathcal{N}\Big{(}D_{\alpha}(p^{n},f)\,,\frac{\sigma^{2}B_{k,\alpha}^{2}}{\alpha^{2}(1-\alpha)^{2}}\Big{)}$
(4.25)

with $N,\,M\to\infty$ and $\sigma^{2}<\infty$.

###  4.2 Asymptotic distribution of Tsallis $\alpha$-divergence estimator

For Tsallis divergence

$T_{\alpha}(p^{n},f)=\frac{1}{\alpha-1}\Big{(}\int\Big{(}\frac{f(x)}{p^{n}(x)}\Big{)}^{1-\alpha}p^{n}(x)\mathrm{d}x-1\Big{)}$

its estimator

$\hat{T}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha-1}\Big{(}\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{(N-1)V(Hs(X_{i},\rho_{k,i}))}{M\>V(Hs(X_{i},\gamma_{k,i}))}\Big{)}^{1-\alpha}\times
B_{k,\alpha}-1\Big{)}$

follows a normal distribution. As the previous estimator, we have the same
procedure. We have now,

$\hat{T}_{\alpha,N,k}(p^{n},f)\xrightarrow{d}\mathcal{N}\Big{(}T_{\alpha}(p^{n},f)\,,\frac{\sigma^{2}B^{2}_{k,\alpha}}{(\alpha-1)^{2}}\Big{)}$
(4.26)

with $N,M\to\infty$ and $\sigma^{2}<\infty$.

###  4.3 Asymptotic distribution for Rényi $\alpha$-divergence estimator

Now back to our estimator of Rényi divergence.

$\hat{R}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha-1}\log\Big{(}\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{\hat{f}_{k,N}(X_{i})}{\hat{p}_{k,n,N}(X_{i})}\Big{)}^{1-\alpha}\times
B_{k,\alpha}\Big{)}$

We had

$h(X_{i})=\Big{(}\frac{\hat{f}_{k,N}(X_{i})}{\hat{p}_{k,n,N}(X_{i})}\Big{)}^{1-\alpha}\,\textrm{and}\quad\bar{Y}_{N}=\frac{1}{N}\sum_{i=1}^{N}h(X_{i})$

$X_{i}$ are i.i.d. and if $h$ is continuous, then $h(X_{i})$ are also i.i.d.

According to the central limit theorem, we will have

$\frac{\bar{Y}_{N}-\mathbb{E}(\bar{Y}_{N})}{\sqrt{Var(\bar{Y}_{N})}}\xrightarrow{d}\mathcal{N}(0,1)$

then,

$\sqrt{N}\big{(}\bar{Y}_{N}-\mathbb{E}(\bar{Y}_{N})\big{)}\xrightarrow{d}\mathcal{N}\big{(}0,N\sigma^{2}\big{)},$

if $\mathbb{E}(\bar{Y}_{N})<\infty$, $\sigma^{2}<\infty$ and $N$ large. We can
apply _delta method_ and we have

$\sqrt{N}\big{(}\log(\bar{Y}_{N})-\log(\mathbb{E}(\bar{Y}_{N}))\big{)}\xrightarrow{d}\mathcal{N}\Big{(}0,\frac{N\sigma^{2}}{[\mathbb{E}(\bar{Y}_{N})]^{2}}\Big{)}.$

Hence

$\hat{R}_{\alpha,N,k}(p^{n},f)\xrightarrow{d}\mathcal{N}\Big{(}\frac{1}{\alpha-1}\log[B_{k,\alpha}\mathbb{E}(\bar{Y}_{N})]\,,\frac{\sigma^{2}}{(\alpha-1)^{2}[\mathbb{E}(\bar{Y}_{N})]^{2}}\Big{)}$
(4.27)

From Póczos at Schneider [11], $B_{k,\alpha}\mathbb{E}(\bar{Y}_{N})\rightarrow
M_{\alpha}(p^{n},f)$, when $N,M\rightarrow\infty$. From which

$\hat{R}_{\alpha,N,k}(p^{n},f)\xrightarrow{d}\mathcal{N}\Big{(}R_{\alpha}(p^{n},f)\,,\frac{\sigma^{2}}{(\alpha-1)^{2}[\mathbb{E}(\bar{Y}_{N})]^{2}}\Big{)}$
(4.28)

##  5 Examples

We will now illustrate our methodology with simple examples. That is why we
will limite ourselves to one-dimensional and two-dimensional cases. In the
following examples , we will use divergence measures to compare proposal
densities corresponding to a given simulation strategy. This difference of
densities may appear at their parameters. The comparison will be for different
parameters values. The proposal density considered as optimum is that which
the divergence measure (function of $n$) between successive densities $p^{n}$
arising and the target density $f$ tends to $0$ faster. We will mainly compare
as we have already told simulation strategies.

###  5.1 One-dimensional case

####  5.1.1 Target density $f$ fully known

In the case where the target density $f$ is analytically known the estimators
of our divergence measures are slightly modified. This change applies at the
constant $B_{k,\alpha}$. This constant is replaced by another in order to get
there also asymptotically unbiased estimators. Then consider

$\hat{M}_{\alpha,N}(p^{n},f)=\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{\hat{f}_{k,N}(X_{i})}{\hat{p}_{k,n,N}(X_{i})}\Big{)}^{1-\alpha}$

If we replace the density’s estimator $f$ by $f$ itself we will have

$\hat{M}^{\prime}_{\alpha,N}(p^{n},f)=\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{f(X_{i})(N-1)c\rho_{k,i}^{d}}{k}\Big{)}^{1-\alpha}$

which is an another estimator of the common integrale part
$M_{\alpha}(p^{n},f)$ but asymptotically biased. As in the proof of Theorem
3.4 (Póczos at Schneider (2011)) we can simply check that the new estimator
multiplied by the constant
$Q_{k,\alpha}=\frac{\Gamma(k-\alpha+1)}{k^{1-\alpha}\Gamma(k)}$ provides an
asymptotically unbiased estimator of $M_{\alpha}(p^{n},f)$.

In one-dimension we are in the set $\mathbb{R}$, therefore $d=1$. Consequently
the d-dimensional hypersphere around $X_{i}$ boils down to the line segment
with $X_{i}$ middle. The real $c$ represents the length measuring of a line
segment with $X_{i}$ middle. This segment has 2 as length measuring since the
distance from $X_{i}$ at each end is 1. As suggested by Loftsgaarden and
Quesenberry [9] we will take $k$ equals to the nearest integer of
$\sqrt{N-1}$. We have now

$\hat{M}^{\prime}_{\alpha,N}(p^{n},f)=\frac{1}{N}\sum_{i=1}^{N}\Big{(}\frac{2\rho_{k,i}(N-1)f(X_{i})}{k}\Big{)}^{1-\alpha}$
(5.29)

We can now have the asymptotically unbiased estimators

$\hat{D}^{\prime}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha(1-\alpha)}\Big{(}1-\hat{M}^{\prime}_{\alpha,N}(p^{n},f)\times
Q_{k,\alpha}\Big{)}$
$\hat{T}^{\prime}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha-1}\Big{(}\hat{M}^{\prime}_{\alpha,N}(p^{n},f)\times
Q_{k,\alpha}-1\Big{)}$

respectively of $D_{\alpha}(p^{n},f)$ and $T_{\alpha}(p^{n},f)$. For the
estimator of the Rényi divergence

$\hat{R}^{\prime}_{\alpha,N,k}(p^{n},f)=\frac{1}{\alpha-1}\log\Big{(}\hat{M}^{\prime}_{\alpha,N}(p^{n},f)\times
Q_{k,\alpha}\Big{)}$

we don’t have results regarding the presence or absence of bias. But the most
important is the fact that three divergence estimators
$\hat{D}^{\prime}_{\alpha,N,k}(p^{n},f)$,
$\hat{T}^{\prime}_{\alpha,N,k}(p^{n},f)$ and
$\hat{R}^{\prime}_{\alpha,N,k}(p^{n},f)$ are both consistent.

First we show an example that compares two proposal densities222Compare
densities is equivalent to compare the respective probability distributions,
so we will use interchangeably densities or probability distributions for a
given strategy. We choose the Independence Sampler (IS), which is one of the
strategies of the Metropolis Hastings algorithm, to compare these densities.

a) Independence Sampler (IS): comparison of proposal densities

For a given simulation strategy the choice of good proposal density is
important. Indeed, for achieving satisfactory simulation results it is
important to choice a good proposal density. For simplicity, we consider
densities that differ by the value of their parameters. Thus we take
respective distribution densities $\mathcal{N}(-3,2)$ and $\mathcal{N}(0,3)$;
the target density is from normale standard distribution $\mathcal{N}(0,1)$.
It is found by looking at Figure 1 that the dashed curve converges very
rapidly to 0 while the solid curve is slow to converge.

Figure 1: Comparison of two proposal densities,using the Independence Sampler
and the Chi-squared divergence ($D_{\alpha}$ with $\alpha=2$)

b) IS - RWMH333Random Walk of Metropolis-Hastings : comparison of simulation
strategies

After finding a good proposal distribution for each strategy, we compare here
the two main strategies of MH algorithm that are IS and RWMH. It may happen in
an experiment that the IS strategy trumps RWMH strategy and in another
experiment the opposite occurs. Everything depends on the instrumental
distribution but also the target distribution to some extent even if the
initial law is the same for both strategies.

Here the target distribution is a gaussian mixture
$0.4\,\mathcal{N}(-8,2)+0.6\,\mathcal{N}(0,6)$. For the _Independence Sampler_
we use the proposal distribution $\mathcal{N}(-2.5,15)$ whereas for RWMH
method we propose to take $\mathcal{N}(x,15)$. This distribution has as mean
equal to the current element $X_{n}=x$. We show the comparison of IS and RWMH
strategies (Fig. 2). Note that the curve associated with IS strategy is below
curve associated with RWMH. However, the two curves converge very quickly to
0.

Figure 2: Comparison of two simulation strategies: Independence Sampler vs
Random Walk of Metropolis-Hastings, using the Chi-squared divergence
($D_{\alpha}$ with $\alpha=2$)

####  5.1.2 Target density $f$ is not known completely

In most real situations, the density $f$ is not known analytically. This is
the case, for example Bayesian context where $f$ is the density of the
posterior. Then $f$ is written as $f=c\varphi$ where $c$ is unknown constant.
[11]

– Adaptive Metropolis444This simulation strategy is proposed by Haario et al.
(2001) \- RWMH

Here the target density is known up to a constant

$f(x)\varpropto\exp(-x^{2})(2+\sin(5x)+\sin(2x))$

Present some _Adaptive Metropolis_ (AM) strategy proposed by Haario et al
(2001). First recall that the stochastic process generated by this simulation
method is not a Markov chain. However it has well ergodicity properties. The
assumptions for this are that the target density is bounded from above and has
a bounded support. Describe the algorithm now.

The target density has a support $E\subset\mathbb{R}^{d}$. Suppose, that at
time $t$ we have sampled the states $X_{0},X_{1},\ldots,X_{t-1}$, where
$X_{0}$ is the initial state. Then a candidate point $Y$ is sampled from the
(asymptotically symmetric) proposal distribution $q_{t}(.|X_{0},...,X_{t-1})$,
which now may depend on the whole history $(X_{0},X_{1},...,X_{t-1})$. The
candidate point $Y$ is accepted with probability

$\alpha(X_{t-1},Y)=\min\Big{(}1,\frac{\pi(Y)}{\pi(X_{t-1})}\Big{)}$

in which case we set $X_{t}=Y$, and otherwise $X_{t}=X_{t-1}$. Observe that
the chosen probability for the acceptance resembles the familiar acceptance
probability of the Metropolis algorithm. The proposal distribution
$q_{t}(.|X_{0},\ldots,X_{t-1})$ employed in the AM algorithm is a Gaussian
distribution with mean at the current point $X_{t-1}$ and covariance

$C_{t}=\left\\{\begin{array}[]{ll}C_{0}\,,&t\leq t_{0}\\\
S_{d}Cov(X_{0},\ldots,X_{t-1})+S_{d}\epsilon
I_{d}\,,&t>t_{0}\,.\end{array}\right.$

where $S_{d}$ is a parameter that depends only on dimension $d$ and
$\epsilon>0$ is a constant that we may choose very small compared to the size
of $E$. Here $I_{d}$ denotes the $d$-dimensional identity matrix. The
covariance $C_{t}$ may be viewed as a function of $t$ variables from
$\mathbb{R}^{d}$ having values in uniformly positive definite matrices.

For these two simulation methods (AM and RWMH) we see that the respective
divergence measures stand very close and very quickly all tend to 0 (Fig. 3).
Indeed, near 0 divergence measures which are in this case the _Hellinger
divergence_ measures $D_{1/2}$ are slight oscillations and intersect.

Figure 3: Comparison of two simulation strategies: Adaptive Metropolis vs
Random Walk of Metropolis-Hastings, using the Hellinger divergence
($D_{\alpha}$ with $\alpha=1/2$)

###  5.2 Two-dimensional case

We consider here the only case of a known function up to a constant. We choose
now a sample $X_{1},\ldots,X_{n}$ which are i.i.d.555independent identically
distributed such that $X_{i}\sim\mathcal{N}(m,\sigma^{2})$. So we have the
following likelihood

$L(x|m,\sigma^{2})\propto(\sigma^{2})^{(-n/2)}\exp\Big{(}-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(x_{i}-m)^{2}\Big{)},$

the prior distributions are

$m\sim\mathcal{N}(m_{0},\sigma_{0}^{2})$
$\sigma^{2}\sim\mathcal{IG}(\alpha,\beta),$

the full posterior density is known up to a constant

$\Pi(m,\sigma^{2}|x)\propto(\sigma^{2})^{-\frac{n}{2}-(\alpha+1)}\exp\Big{(}-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(x_{i}-m)^{2}-\frac{(m-m_{0})^{2}}{2\sigma_{0}^{2}}-\frac{\beta}{\sigma^{2}}\Big{)},$

the conditional distributions of parameters are

$m|\sigma^{2},x\sim\mathcal{N}(M,\Sigma^{2})$

where

$M=\frac{\sigma_{0}^{2}\sum_{i=1}^{n}x_{i}+\sigma^{2}m_{0}}{\sigma^{2}+n\sigma_{0}^{2}}\quad\textrm{and}\quad\Sigma^{2}=\frac{\sigma^{2}\sigma_{0}^{2}}{\sigma^{2}+n\sigma_{0}^{2}}$
$\sigma^{2}|m,x\sim\mathcal{IG}\Big{(}\frac{n}{2}+\alpha,\frac{1}{2}\sum_{i=1}^{n}(x_{i}-m)^{2}+\beta\Big{)}$

So let’s compare firstly RWMH and Gibbs sampler and secondly the RWMH and
Metropolis Within Gibbs. $\Pi(m,\sigma^{2}|x)$ will be our target density.

a) RWMH - Gibbs sampler

If RWMH applied both in dimension $d\geq 1$ , the Gibbs sampler on it only
applied in dimension $d>1$. However, we limit ourselves here in dimension 2.
Note that this method (Gibbs sampler) has been used by Geman (1984) to
generate observations from a Gibbs distribution (Boltzmann distribution). It
is a particular form of the MCMC method, because of its effectiveness, is
widely used in many fields of Bayesian analysis. Thus to simulate according to
a probability density $f(\theta)$ with $\theta=(\theta_{1},\ldots,\theta_{p})$
one can use the following idea

Initialisation: generating a vector $\theta=(\theta_{1},\ldots,\theta_{p})$
according to a initial proposal law $\Pi_{0}$.

Simulate following the conditional distributions

$\theta_{i}|\theta_{1},\ldots,\theta_{i-1},\theta_{i+1},\ldots,\theta_{p}\sim
f_{i}(\theta_{i}|\theta_{1},\ldots,\theta_{i-1},\theta_{i+1},\ldots,\theta_{p})$

$i=1,2,\ldots,p$.

We see that the curve corresponding to the RWMH strategy is well above that
representing the Gibbs sampler (Fig. 4). The two curves have only one common
point that is their origin. As mentioned in the legend, the dashed curve is
associated with RWMH while the solid curve is associated with the Gibbs
sampler. However, the curves denote divergence measures, we will talk about it
in Section 6 (Discussion).

Figure 4: Comparison of two simulation strategies: Gibbs sampler vs Random
Walk of Metropolis-Hastings, using the Tsallis divergence with $\alpha=0.99$

b) RWMH - Metropolis Within Gibbs

_Metropolis Within Gibbs_ is a hybrid simulation method that combines stages
of the Gibbs sampler and Metropolis Hastings method. It is used in some cases
where we have conditional distributions for which we can’t have samples
directly. There are several versions of this sampler, so we present the
following.

Assume that $\pi(.)$ is the target distribution666$\pi(.)$ and $q(.)$ are the
probability densities functions , $\pi(.|x_{-i})$ denote now the conditional
distribution of $Z|Z_{-i}=z_{-i}$ where $Z\sim\pi$.
$X_{n}:=(X_{n,1},\ldots,X_{n,d})$;
$X_{n,-i}:=(X_{n,1},\ldots,X_{n,i-1},X_{n,i+1},X_{n,d})$;
$\alpha:=(\alpha_{1},\ldots,\alpha_{d})$.

Now we have the following algorithm.

Algorithm
1\. Choose coordinate $i\in\\{1,\ldots,d\\}$ according to selection
probabilities $\alpha$, that is, with $\mathbb{P}(i=j)=\alpha_{j}$ .
2\. Draw $Y\sim q(X_{n-1,i},.)$.
3\. Accepte the candidate $Y$ with probability

$\min\Big{(}1,\frac{\pi(Y|X_{n-1,-i})q(Y,X_{n-1,i})}{\pi(X_{n-1,i}|X_{n-1,-i})q(X_{n-1,i},Y)}\Big{)}$

and set $X_{n}:=(X_{n-1,1},\ldots,X_{n-1,i-1},Y,X_{n-1,i+1},\ldots,X_{n-1,d})$
otherwise reject $Y$ and set $X_{n}=X_{n-1}$.

Starting from a common point, the curves stay away from the value 0 (Fig. 5).
The solid curve is associated with _Metropolis Within Gibbs_ strategy and the
dashed curve associated with _RWMH_ strategy. Here we use the _Rényi
divergence_ measure in order to obtain our curves.

Figure 5: Comparison of two simulation strategies:_Metropolis within Gibbs_ vs
_Random Walk of Metropolis-Hastings_ , using the Rényi $\alpha$-divergence
with $\alpha=0.3$

##  6 Discussion

###  6.1 Indepence Sampler : comparison of proposal densities

The two curves represent measurements of _Chi-square divergence_ $D_{2}$.
Indeed, the dashed curve indicates the difference between the densities
$p^{n}_{1}$ from the IS strategy with proposal law $\mathcal{N}(0,3)$ and the
target density $f$ of $\mathcal{N}(0,1)$, noted by $D_{2}(p^{n}_{1},f)$. The
solid curve denotes the difference between the densities $p^{n}_{2}$ produced
by the IS strategy with proposal distribution $\mathcal{N}(-3,2)$ and the
target distribution $\mathcal{N}(0,1)$, noted by $D_{2}(p^{n}_{2},f)$.

It is clear that $D_{2}(p^{n}_{1},f)$ and $D_{2}(p^{n}_{2},f)$ are functions
of the number of iterations $n$ and have the same origin. The fact that they
have the same starting point can be explained by relevance to begin with a
same state $x_{0}$ or same initial law $\pi_{0}$ to compare two simulation
strategies. The two curves overlap up to the first iteration. Indeed, the IS
strategy with $\mathcal{N}(0,3)$ is more efficient because this distribution
has a mean $m=0$ like the mean of target distribution. Its variance
$\sigma^{2}=3$ is greater than the variance of the target distribution which
is $\sigma^{2}=1$. It follows that its support covers the support of the
target density. The other strategy has a proposal distribution
$\mathcal{N}(-3,2)$. This law has a mean $m=-3$ and a variance $\sigma^{2}=2$
that make its density function is shifted to the left (see figure) and is
therefore not adequate to cover the support of the target density. So that the
chain generated by this simulation strategy will soon converge to the target
distribution.

###  6.2 Independence Sampler - Random Walk of Metropolis Hastings

We said earlier that the two curves in Figure 2 converge rapidly to 0.We use
also, here _Chi-square divergence_.Dashed curve represent also the divergence
measure between the densities $p^{n}$ of RWMH strategy and target density
$f(x)=0.4f_{1}(x,-8,2)+0.6f_{2}(x,0,6)$ which is Gaussian mixture density.
Similarly, the solid curve represents divergence measure between the densities
$p^{n}$ from the IS strategy and the target density $f(x)$.

Although the two divergence measures converge rapidly to $0$, we see that the
divergence which represents the IS strategy is below the divergence associated
with RWMH strategy. Therefore IS strategy is more efficient, even if RWMH is
also good method.

Note that the two curves have the same starting point and overlap untill the
first iteration. From there, the curve of IS strategy approximates much to
$0$. This is due to the fact this method (IS) has, here, a good proposal
distribution $\mathcal{N}(-2.5,15)$. With a mean $m=-2.5$ and a variance
$\sigma^{2}=15$, its density function is centered relatively to the target
density $f(x)$. Thus the support of the target density is well covered by the
proposal density.

Regarding the RWMH strategy, its has a normal proposal distribution to each
iteration $n$ with a mean equal to the current state $X_{n}=x_{n}$ and
variance also equal to $\sigma^{2}=15$. That’s why even if its has not a
support that covers all the time the target distribution’s support, its does
not so far away. Hence, this is also a good strategy.

###  6.3 Adaptive Metropolis - RWMH

For this comparison, we use the _Hellinger divergence_ measure $(D_{1/2})$.
The solid curve then describes the divergence measure between densities
$p^{n}$ (AM) and the target density $f(x)$. The dashed curve described the
divergence between the densities $p_{n}$ (RWMH)and also target density. The
two curves merge between the initial state and the first iteration. From
there, they don’t move away from each other, but coexist around value $0$.
Thus, the two corresponding simulation strategies are all very efficient.

The effectiveness of these two strategies is explained in first concerns the
AM strategy, by the fact that it adapts its proposal density to the target
density. The adjustment mechanism is performed at the variance of instrumental
density. Indeed, if $X_{1},X_{2},\ldots,X_{n-1}$ have already been simulated
and one desire to obtain the point (or vector) $X_{n}$ at time $n$, we
generate with a proposal distribution with mean equal the current state of
$X_{n-1}$ and covariance martice described above.

Recall that we use the AM strategy here in one dimension. Thus, the variance
$\sigma^{2}=5$ chosen for our instrumental Gaussian between the initial time
and time $t=15$ ( this is the period prior to the iterative update of the
variance ) yields samples which describe almost support the target density.
Added to this, from the 16th iteration adjustments allow the variance to have
good samples. Which explains why it has a very fast convergence. Concerning
RWMH method, with same variance $\sigma^{2}=5$ for the proposal density, we
also obtain good samples that converge very quickly to the stationary density.

###  6.4 Gibbs Sampler - RWMH

As we said earlier, we have in Figure 4 two curves designating measures of
Tsallis divergence with $\alpha=0.99$. The dashed curve (RWMH) is above the
other and far from the value $0$ after the first 50 iterations. It shows,
here, that RWMH strategy is ineffective . This inefficiency is due to the
covariance matrix $M$ of the proposal distribution. This matrix which is
implementation parameter is not optimum. The solid curve (Gibbs sampler) tends
rapidly to 0 (after the 7th iteration). The strength of this sample is mainly
due to the fact that the components of the vector $X_{n}$ are generated
directly from the simulated conditional distributions of the target
distribution.

###  6.5 Metropolis Within Gibbs - RWMH

The Rényi divergence measures between the respective densities $p^{n}$ both
strategies (MWG and RWMH) and the stationary distribution are held away from
the value 0 even after a large number of simulations. What makes us to say
that two strategies are not effective.

As in the previous comparison (Gibbs Sampler vs RWMH) RWMH strategy is always
slow in computation time. In fact it is the same strategy with the
instrumental Gaussian distribution with parameter matrix $M$, but evaluated
using the Rényi divergence measure with $\alpha=0.3$. We see thereby that we
have a wide range of divergence measures allowing us to evaluate a given
strategy simulation. Respecting the MwG method in our case, even after 30,000
iterations, the created process (Markov chain) does not converge (Fig 5); even
if we are limited here in 1,000 iterations.

This method (MwG) is often used in cases where one wants to use the Gibbs
sampler and for some conditional distributions he can not simulate directly.
Then one introduce in the algorithm a few steps of RWMH sampler, offering
instrumental distributions having for target laws the conditional
distributions. In our case (Fig. 5) we used a version that systematically
applies the steps of the Metropolis sampler to all conditional distributions,
so that the algorithm will delay to converge.

###  6.6 Conclusion

We have, therefore, shown that with various divergence measures we can compare
two different simulation strategies. To achieve this, we used
$\alpha$-divergence measures that we have described in Section 2. We showed
the convergence of these divergence measures and proposed estimators for these
ones (Section 3). In Section 4 we gave the asymptotic distribution of each
estimator. Then we gave some examples for the implementation of simulation
strategies that have been described (Section 5). Finally we discussed in
Section 6 of the causes that make a simulation strategy is more powerful than
the other.

## References

  * [1] Y. D. Kim A. Cichocki, H. Lee and S. Choi.  Non-negative matrix factorization with $\alpha$-divergence.  Pattern Recognition Letters, 2008.
  * [2] N. K. Verma B. Mehta and P. Sircar.  Performance analysis of alpha divergence in nonnegative matrix factorization of monaural musical sounds.  International Journal of Engineering, Science and Technology, 3(6):273–282, 2011.
  * [3] J. Xiong B. Póczos and J. Schneider.  Nonparametric Divergence Estimation with Applications to Machine Learning on Distributions.  arXiv:1202.3758v1 [cs.LG], pages 599–608, 2012.
  * [4] D. Chauveau and P. Vandekerkhove.  How to compare mcmc simulation strategies?  hal-00019174 (version 3), 2007.
  * [5] D. Chauveau and P. Vandekerkhove.  Smoothness of metropolis-hastings algorithm and application to entropy estimation.  ESAIM: Probability and Statistics, 2012.
  * [6] A. Cichocki and S. Amari.  Families of Alpha- Beta- and Gamma- Divergences: Flexible and Robust Measures of Similarities.  Entropy, pages 1–41, 2010.
  * [7] E. Saksman H. Haario and J. Tamminen.  An adaptive Metropolis algorithm.  Bernoulli, 7(2):223–242, 2001.
  * [8] G. O. Roberts K. Latuszyński and J. S. Rosenthal.  Adaptive Gibbs sampler and related MCMC methods.  The Annals of Applied Probability, 23(1):66–98, 2013.
  * [9] D. O. Loftsgaarden and C. P. Quesenberry.  A Nonparametric Estimate of a Multivariate Density Function.  Annals of Mathematical Statistic, 36(3):1049–1051, 1965.
  * [10] Annie MILLET.  Méthodes de monte-carlo.  Universités Paris 7 et Paris 1, Cours de Master 2: spécialité Modélisation Aléatoire.
  * [11] B. Póczos and J. Schneider.  On the Estimation of Divergences.  In Proceedings of the 14th International Conference on Artificial Intelligence and Statistics (AISTATS), volume 15, pages 609–617, 2011\.
  * [12] T. van Erven and P. Harremoës.  Rényi Divergence and Majorization.  arXiv: 1001.4448v3 [cs.IT], 2010.