Mireuta Matei 2011 Memoire

Université de Montréal
Étude de la performance d’un algorithme Metropolis-Hastings avec ajustement

directionnel
par
Matei Mireuta
Département de mathématiques et de statistique

Faculté des arts et des sciences
Mémoire présenté à la Faculté des études supérieures

en vue de l’obtention du grade de Maître ès sciences (M.Sc.)
en mathématiques
août, 2011
c Matei Mireuta, 2011.

Université de Montréal
Faculté des études supérieures
Ce mémoire intitulé:
Étude de la performance d’un algorithme Metropolis-Hastings avec ajustement

directionnel
présenté par:
Matei Mireuta
a été évalué par un jury composé des personnes suivantes:
Pierre Lafaye de Micheaux, président-rapporteur

Mylène Bédard, directrice de recherche
Jean-François Angers, membre du jury
Mémoire accepté le: . . . . . . . . . . . . . . . . . . . . . . . . . .

RÉSUMÉ
Les méthodes de Monte Carlo par chaîne de Markov (MCMC) sont des outils très popu-
laires pour l’échantillonnage de lois de probabilité complexes et/ou en grandes dimen-
sions. Étant donné leur facilité d’application, ces méthodes sont largement répandues
dans plusieurs communautés scientifiques et bien certainement en statistique, particuliè-
rement en analyse bayésienne. Depuis l’apparition de la première méthode MCMC en
1953, le nombre de ces algorithmes a considérablement augmenté et ce sujet continue
d’être une aire de recherche active.
Un nouvel algorithme MCMC avec ajustement directionnel a été récemment déve-
loppé par Bédard et al. (IJSS, 9 :2008) et certaines de ses propriétés restent partiellement
méconnues. L’objectif de ce mémoire est de tenter d’établir l’impact d’un paramètre clé
de cette méthode sur la performance globale de l’approche. Un second objectif est de
comparer cet algorithme à d’autres méthodes MCMC plus versatiles afin de juger de sa
performance de façon relative.
Mots clés: échantillonneur indépendant, algorithme Metropolis-Hastings de type
marche aléatoire, taux de convergence, algorithme Metropolis adaptatif, candidats
multiples, diagnostics de convergence.
ABSTRACT
Markov Chain Monte Carlo algorithms (MCMC) have become popular tools for sam-
pling from complex and/or high dimensional probability distributions. Given their rel-
ative ease of implementation, these methods are frequently used in various scientific
areas, particularly in Statistics and Bayesian analysis. The volume of such methods has
risen considerably since the first MCMC algorithm described in 1953 and this area of
research remains extremely active.
A new MCMC algorithm using a directional adjustment has recently been described
by Bédard et al. (IJSS, 9:2008) and some of its properties remain unknown. The objec-
tive of this thesis is to attempt determining the impact of a key parameter on the global
performance of the algorithm. Moreover, another aim is to compare this new method to
existing MCMC algorithms in order to evaluate its performance in a relative fashion.
Keywords: independent sampler, random walk Metropolis-Hastings algorithms,
convergence rate, adaptive Metropolis algorithm, multiple proposals, convergence
diagnostics.
TABLE DES MATIÈRES
RÉSUMÉ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
TABLE DES MATIÈRES . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
LISTE DES TABLEAUX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
LISTE DES FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
LISTE DES ANNEXES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x
DÉDICACE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
REMERCIEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii
INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii
CHAPITRE 1 : INTRODUCTION AUX ALGORITHMES MCMC . . . . 1

1.1 Construction des algorithmes MCMC . . . . . . . . . . . . . . . . . . 1
1.2 Algorithmes MCMC élémentaires . . . . . . . . . . . . . . . . . . . . 4
1.3 Convergence des algorithmes MCMC . . . . . . . . . . . . . . . . . . 7
1.4 Taux de convergence des algorithmes MCMC . . . . . . . . . . . . . . 9
1.5 Théorème central limite des algorithmes MCMC . . . . . . . . . . . . 11
1.6 Diagnostics de convergence . . . . . . . . . . . . . . . . . . . . . . . . 15
1.6.1 Graphique de la trace d’un paramètre d’intérêt . . . . . . . . . 15
1.6.2 Densité (histogramme) de valeurs générées et chaînes parallèles 16
1.6.3 Taux d’acceptation . . . . . . . . . . . . . . . . . . . . . . . . 16
1.6.4 Autocorrélation et distance de saut carrée moyenne . . . . . . . 17
vi
1.6.5 Statistiques de convergence . . . . . . . . . . . . . . . . . . . 20

1.7 Variations d’algorithmes MCMC . . . . . . . . . . . . . . . . . . . . . 20
CHAPITRE 2 : ALGORITHME DE TYPE METROPOLIS-HASTINGS AVEC

AJUSTEMENT DIRECTIONNEL . . . . . . . . . . . . . 22
2.1 L’algorithme DA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2 Étapes de l’algorithme DA . . . . . . . . . . . . . . . . . . . . . . . . 28
2.3 Impact de la valeur de λ sur la performance de l’algorithme DA . . . . 30
2.3.1 Exemple 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3.2 Exemple 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.4 Convergence de l’algorithme DA . . . . . . . . . . . . . . . . . . . . . 39
CHAPITRE 3 : COMPARAISON ENTRE L’ALGORITHME DA ET DES

MÉTHODES LOCALES . . . . . . . . . . . . . . . . . . 50
3.1 Algorithmes RWMH . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.2 Algorithme Metropolis adaptatif . . . . . . . . . . . . . . . . . . . . . 58
3.3 Algorithmes Metropolis avec essais multiples . . . . . . . . . . . . . . 60
3.3.1 Algorithme MTM . . . . . . . . . . . . . . . . . . . . . . . . 61
3.3.2 Algorithme MTM hit-and-run . . . . . . . . . . . . . . . . . . 62
3.3.3 Algorithme Metropolis avec rejet différé . . . . . . . . . . . . . 63
3.4 Conclusion et analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
BIBLIOGRAPHIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
LISTE DES TABLEAUX
3.I Résumé des résultats de méthodes locales (exemple 2) (Période de

chauffe 10 000, Itérations 4 000 000). . . . . . . . . . . . . . . . 67
3.II Résumé des résultats de méthodes locales (exemple 2) avec ajus-
tement de la matrice de covariance (Période de chauffe 10 000,
Itérations 4 000 000). . . . . . . . . . . . . . . . . . . . . . . . . 70
3.III Résumé des résultats de méthodes locales (exemple 1) (Période de
chauffe 10 000, Itérations 4 000 000). . . . . . . . . . . . . . . . 72
LISTE DES FIGURES
1.1 Graphique de la trace de paramètres en fonction des itérations. . . 16

1.2 Graphique en boîte de valeurs générées par deux chaînes parallèles
avec valeur initiale différente. . . . . . . . . . . . . . . . . . . . 17
1.3 Graphique de l’autocorrélation de valeurs générées. . . . . . . . . 18
2.1 Graphique de la densité cible (ligne pleine) et instrumentale (ligne

pointillée) lorsque la direction choisie est vers la droite (+1). Le
point d’intersection est u prop ∗ ∗
j+1 · s = s . . . . . . . . . . . . . . . . 29
2.2 Graphique de l’écart-type σ de s(β = 1) en fonction de λ (algo-
rithme DA, exemple 1). . . . . . . . . . . . . . . . . . . . . . . . 32
2.3 Graphique du taux d’acceptation en fonction de λ (algorithme DA,
exemple 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.4 Graphique de la DSCM en fonction de λ (algorithme DA, exemple
1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.5 Graphique des degrés de liberté moyens (avec un écart-type) en
fonction de λ (algorithme DA, exemple 1). . . . . . . . . . . . . 35
2.6 Graphique de la DSCM des algorithmes IS et DA (exemple 1). . . 36
2.7 Graphique du taux d’acceptation des algorithmes IS et DA (exemple
1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.8 Graphique de l’écart-type σ de s(β = 1) des algorithmes IS et DA
(exemple 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.9 Graphique de la DSCM en fonction de λ (algorithme DA, exemple
2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.10 Graphique du taux d’acceptation en fonction de λ (algorithme DA,
exemple 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
ix
2.11 Graphique de l’écart-type σ de la valeur-p pour H0 : β6 = −0, 01

(algorithme DA, exemple 2). . . . . . . . . . . . . . . . . . . . . 42
2.12 Graphique des degrés de liberté moyens en fonction de λ (algo-
rithme DA, exemple 2). . . . . . . . . . . . . . . . . . . . . . . . 43
2.13 Graphique de la fonction de densité (2.12). . . . . . . . . . . . . 45
2.14 Graphique des degrés de liberté proposés, de la valeur s(x = 1)
obtenue (avec un écart-type) et de l’autocorrélation en fonction de
λ (algorithme DA, exemple 2.12). . . . . . . . . . . . . . . . . . 47
3.1 Graphique des valeurs-p obtenues par la méthode DA (ligne pleine),

la méthode RWMH (lignes pointillées) ainsi que les approxima-
tions de troisième ordre (astérisques) en fonction des hypothèses
H0 sur β6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.2 Graphique de l’autocorrélation des valeurs générées par la mé-
thode RWMH. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.3 Graphique de la DSCM en fonction du facteur multiplicatif de la
matrice de covariance (algorithme RWMH, exemple 2). . . . . . 54
3.4 Graphique du taux d’acceptation en fonction du facteur multipli-
catif de la matrice de covariance (algorithme RWMH, exemple 2).
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
LISTE DES ANNEXES
Annexe I : Programmes R . . . . . . . . . . . . . . . . . . . . . . . . xvi

À Marina, Ioana, Neculai, Hunter et Askim
REMERCIEMENTS
Je tiens avant tout à remercier ma directrice de recherche, Mme Mylène Bédard, qui
a été d’un immense soutien pendant toute la durée de ce projet. Grâce à sa perspicacité
et à son expertise remarquables, elle a été une grande source d’inspiration et de conseils.
Elle a toujours su se montrer disponible, à l’écoute et a su partager ses connaissances
de manière exceptionnelle. J’aimerais aussi remercier les membres du jury pour la cor-
rection de ce mémoire et pour leurs commentaires qui sauront améliorer la qualité de ce
travail.
INTRODUCTION
En statistique, un défi récurrent est d’obtenir des descriptions adéquates d’une va-
riable aléatoire étant donné sa fonction de densité. En effet, il arrive souvent qu’il y ait
un intérêt pour l’espérance d’une variable aléatoire, pour sa variance, ses quantiles ou
toute autre mesure significative. Parfois, il est nécessaire de calculer une valeur-p, une
valeur-s ou tout simplement d’acquérir un échantillon d’une distribution donnée. Bien
certainement, il existe une multitude d’approches à ces problèmes et chaque contexte
particulier déterminera la meilleure solution.
Afin de répondre à certaines de ces questions, il est possible d’utiliser des mé-
thodes d’intégration numérique, aussi appelées quadratures. Dans le cas unidimension-
nel, il existe des méthodes classiques telles que l’algorithme de Newton-Cotes ou des
approches plus puissantes telles que la quadrature de Gauss-Legendre. Ces méthodes
peuvent être modifiées afin d’accommoder les cas multidimensionnels qui sont beau-
coup plus fréquents, par exemple en traitant une intégrale multiple comme une série
d’intégrales unidimensionnelles. D’un autre côté, il est possible d’utiliser certains al-
gorithmes développés récemment tel que la méthode introduite par Genz (1972). Ces
approches se comportent de façon excellente lorsque la dimension de l’intégrale reste
relativement petite, mais sont peu utiles en dimension modérée ou grande. En effet, la
quantité de ressources nécessaires augmente de façon exponentielle avec la dimension
du problème, menant à une durée d’exécution déraisonnable au-delà d’une dimension de
4 ou 5.
Il est également possible d’employer des méthodes de simulation, aussi appelées
méthodes Monte Carlo. L’idée de base est d’obtenir un échantillon aléatoire de la dis-
tribution d’intérêt et ensuite d’estimer les quantités voulues de façon empirique en se
servant de l’échantillon généré. Cette approche repose sur deux conditions qui sont la
capacité de facilement générer des valeurs de la distribution en question et la capacité
de produire un gros échantillon afin d’obtenir des résultats fiables. Il existe un nombre
xiv
considérable d’algorithmes de simulation, comme entre autres, la méthode du rejet ou la

méthode d’échantillonnage d’importance.
Les méthodes de Monte Carlo par chaîne de Markov (MCMC) constituent une des
approches les plus utilisées dans la communauté statistique. Ces méthodes emploient une
chaîne de Markov auxiliaire dont la distribution stationnaire est la distribution d’intérêt.
L’algorithme Metropolis-Hastings, introduit par Metropolis et al. (1953) et généralisé
par Hastings (1970), est considéré comme l’une des premières méthodes MCMC. De-
puis, le nombre de tels algorithmes ainsi que le nombre de publications portant sur leur
convergence et leur application a augmenté de façon remarquable. Le principal attrait
des approches MCMC est leur facilité d’application à des distributions d’intérêt com-
plexes et/ou en grandes dimensions. En plus, un autre grand avantage est le fait que la
constante de normalisation de la densité d’intérêt ne doit généralement pas être spéci-
fiée. Cet aspect représente un attrait important pour la communauté bayésienne, car cette
constante est souvent difficilement obtenable pour une fonction de densité à postériori.
Les algorithmes de simulation autre que les méthodes MCMC requièrent, en général,
au moins une approximation ou des bornes pour cette valeur, rendant leur application
plus ardue. Pour ces raisons, les algorithmes MCMC ont largement gagné en popularité
dans les cinquante dernières années, laissant quelque peu en arrière-plan les méthodes
d’échantillonnage classiques ou bien les méthodes par quadrature.
En première partie, ce mémoire fera une description générale de l’approche MCMC.
Il sera question de la construction de tels algorithmes, de leur convergence ainsi que
du théorème central limite s’appliquant aux méthodes MCMC. Certains algorithmes
de base seront aussi décrits en détails et finalement, plusieurs mesures empiriques de
convergence seront présentées.
La deuxième partie traitera d’un algorithme Metropolis avec ajustement direction-
nel développé récemment par Bédard et Fraser (2008). Cette méthode sera d’abord dé-
crite d’un point de vue théorique ainsi que d’un point de vue pratique. Étant donné la
nouveauté de cette approche, plusieurs de ses propriétés restent partiellement mécon-
xv
nues. Notre travail tentera donc d’élucider quelques aspects de la performance et de la

convergence de cette méthode. Pour ce faire, nous utiliserons, en un premier lieu, deux
exemples réalistes afin de démontrer l’importance du choix d’un certain paramètre λ
intrinsèque à cette approche. Ensuite, nous construirons un exemple unidimensionnel
afin de montrer que cette méthode peut parfois présenter des problèmes de convergence
importants.
Enfin, le dernier chapitre portera sur une comparaison de la performance de l’al-
gorithme Metropolis avec ajustement directionnel et d’autres approches MCMC. Notre
comparaison sera faite surtout avec des méthodes locales, plus précisément les algo-
rithmes Metropolis adaptatifs, avec essais multiples et avec rejet différé. En un même
temps, nous complèterons également l’analyse d’un exemple décrit dans Bédard et Fra-
ser (2008).
CHAPITRE 1
INTRODUCTION AUX ALGORITHMES MCMC
Les algorithmes MCMC représentent une méthode alternative pour l’échantillon-

nage de variables aléatoires ayant des densités complexes et/ou en grandes dimensions.
Ce chapitre servira d’introduction aux approches MCMC et détaillera certaines des prin-
cipales propriétés théoriques. En un premier lieu, il sera question de la construction de
ces algorithmes et de quelques applications de base. Ensuite, les sections suivantes por-
teront sur leur convergence ainsi que sur certaines propriétés asymptotiques. Finalement,
la dernière section sera consacrée à des mesures de convergence empiriques.
1.1 Construction des algorithmes MCMC
Un espace X et une fonction de densité πu (x), possiblement non normalisée, mais

R
satisfaisant 0 < X πu (x) < ∞ sont donnés. Bien que d’autres situations soit possibles, le
présent mémoire considérera seulement le cas le plus courant, c’est-à-dire où X est un
sous-ensemble ouvert de Rn et la mesure sous-jacente est Lebesgue. Dans ce cas, pour
tout A ⊆ X mesurable, la mesure de probabilité π(·) sur l’espace X est définie par
R
πu (x)dx
π(A) = RA .
X πu (x)dx
Afin de produire un échantillon de la distribution π(·), l’approche MCMC consiste
à générer une chaîne de Markov qui aura comme loi stationnaire π(·). En d’autres mots,
il s’agit de trouver des probabilités de transition P(x, dy) = P(X j+1 ∈ dy|X j = x) (dy est
un voisinage infinitésimal de y) pour la chaîne telles que
Z
π(dx)P(x, dy) = π(dy) ∀x, y ∈ X. (1.1)
x∈X
Conséquemment, si la chaîne est simulée assez longtemps, la distribution du j-ième

état X j est approximativement stationnaire L (X j ) ≈ π(·) (sujet à certaines conditions
2
mentionnées à la section 1.3) et une première observation aléatoire Z1 = X j est alors

obtenue. Il est ensuite possible d’utiliser une seconde valeur de départ (possiblement la
même) et de simuler de nouveau la chaîne de manière identique afin d’obtenir Z2 . Le
processus peut être répété en fonction de la taille échantillonnale souhaitée.
Bien que cette approche paraisse attrayante, construire une chaîne de Markov avec
une telle propriété peut sembler une tâche ardue à priori. En fait, en utilisant certains
concepts clé, il est étonnamment simple de bâtir de telles chaînes. Une idée importante
est celle de la réversibilité.
Définition 1. Une chaîne de Markov avec probabilités de transition P(x, dy) sur un
espace X est dite réversible par rapport à une distribution π(·) si
π(dx)P(x, dy) = π(dy)P(y, dx) ∀x, y ∈ X. (1.2)
Ensuite, il est facile de voir que lorsqu’une chaîne de Markov est réversible par rap-
port à π(·), la distribution π(·) est stationnaire pour cette chaîne :
Z Z
π(dx)P(x, dy) = π(dy)P(y, dx)
x∈X x∈X Z
= π(dy) P(y, dx) = π(dy) . (1.3)

x∈X
Donc, en théorie, il suffit de construire une chaîne de Markov réversible par rapport
à π(·) et de simuler suffisamment d’états afin d’obtenir une observation aléatoire Z1 ∼
π(·). Une manière simple d’assurer la réversibilité de la chaîne est d’utiliser l’algorithme
de Metropolis-Hastings (Metropolis et al. (1953); Hastings (1970)).
Cette méthode se sert d’une chaîne de Markov auxiliaire, simple à implémenter, di-
sons Q(x, ·), dont la densité de transition n’est pas nécessairement normalisée, Q(x, dy) ∝
q(x, y)dy. En premier, X0 = x0 est choisi comme état initial. Ensuite, étant donné un état
X j = x, un nouvel état y est proposé selon la loi Q(x, ·) et une variable indépendante B
est générée, où B ∼ Bernoulli(α) avec probabilité de succès définie par :

πu (y)q(y, x)
α(x, y) = min 1, . (1.4)
πu (x)q(x, y)
3
Si B = 1, l’état y est accepté et devient le nouvel état de la chaîne (X j+1 = y). Si B = 0,

l’état y est rejeté et x est posé comme nouvel état (X j+1 = x).
De façon pratique, ce processus crée une nouvelle chaîne de Markov {X j } qui a
la propriété additionelle d’être réversible par rapport à π(·). En fait, la réversibilité est
assurée par la variable aléatoire B à travers la forme de la probabilité de succès (1.4).
Pour vérifier la propriété (1.2), supposons que x 6= y (le cas x = y est trivial) et déno-
R
tons c = X πu (x)dx. Il s’en suit alors que :
π(dx)P(x, dy) = c−1 πu (x)dx q(x, y)α(x, y)dy

−1 πu (y)q(y, x)
= c πu (x)q(x, y) min 1, dxdy
πu (x)q(x, y)
= c−1 min (πu (x)q(x, y), πu (y)q(y, x)) dxdy,
qui est symétrique en x et y (voir Metropolis et al. (1953)).

Donc, la construction d’une chaîne de Markov avec distribution stationnaire π(·)
est bel et bien faisable. Par conséquent, il suffit de choisir un état initial et ensuite de
simuler une chaîne de Markov selon l’algorithme Metropolis-Hastings plusieurs fois afin
de générer un échantillon aléatoire provenant de π(·). En fait, bien qu’il soit possible
de générer un échantillon aléatoire en réinitialisant la chaîne pour chaque observation
souhaitée, il est souvent préférable de garder la queue d’une seule chaîne. En pratique,
une première tranche d’observations, appelée période de chauffe (Burn-in), est laissée
de côté parce que l’on considère que l’algorithme n’a pas encore convergé. Au-delà de
la période de chauffe, les valeurs sont considérées comme un échantillon provenant de
la distribution π(·). En vérité, les valeurs ainsi générées sont dépendantes, mais cela
n’affecte généralement pas de façon considérable les quantités à calculer (souvent des
espérances) et le processus de simulation est beaucoup plus efficace, exigeant moins de
ressources. La distribution π(·) et la densité cπu (x) sont appelées distribution et densité
cible tandis que Q(x, ·) et cq(x, ·) sont appelées distribution et densité instrumentales
(c étant la constante de normalisation). Afin de ne pas introduire une notation superflue,
π(·) se référera à la densité ainsi qu’à la distribution cible et q(x, ·) se référera à la densité
4
ainsi qu’à la distribution instrumentale au-delà de ce chapitre introductoire.

Une autre caractéristique importante de l’approche Metropolis-Hastings est le fait
que la constante de normalisation de la densité cible ne doit pas nécessairement être
connue puisque cette densité n’apparaît que sous forme de ratio dans l’expression (1.4).
Cela est très important dans le contexte bayésien, où cette constante est souvent incalcu-
lable. En effet, un des objectifs de l’approche bayésienne est d’inférer sur la distribution
à postériori d’un paramètre θ . En utilisant une densité à priori p(θ ), la vraisemblance
des réalisations observées L(x|θ ) et en appliquant le théorème de Bayes, la forme de la
densité à postériori de θ est donnée par :
πu (θ ) ∝ L(x|θ )p(θ ) .
R
La constante de normalisation est θ L(x|θ )p(θ )dθ . Dans la plupart des algorithmes
de simulation Monte Carlo, le calcul de cette constante est requis pour générer des obser-
vations aléatoires de la distribution à postériori. Puisque cette étape n’est pas nécessaire
dans l’approche MCMC, l’utilisation de ce type d’algorithme est extrêmement répandue
dans la communauté bayésienne.
1.2 Algorithmes MCMC élémentaires
Il existe une multitude de façons de construire la chaîne de Markov instrumentale,

mais cette section ne discutera que des plus usitées. La plupart des méthodes suivantes
sont présentées pour le cas unidimensionnel, mais la généralisation en n dimensions est
immédiate.
(a) Algorithme de type Metropolis-Hastings avec marche aléatoire (random walk Metro-
polis-Hastings (RWMH))
La méthode RWMH comporte une chaîne de Markov instrumentale dont la densité
de transition dépend de l’état présent, plus précisément elle satisfait q(x, y) = q(y − x).
Par exemple, on peut choisir une densité instrumentale N(x, σ 2 ). Donc, étant donné un
5
état présent x, l’algorithme génère un état potentiel y provenant d’une distribution nor-
male centrée à x et avec variance σ 2 . Si le nouvel état y est accepté, le prochain état
potentiel sera généré selon une distribution N(y, σ 2 ). Sinon, x sera posé comme état
actuel et un autre état y éventuel sera proposé selon une loi N(x, σ 2 ). Ce type d’algo-
rithme est très versatile puisque son application nécessite peu d’information au sujet
de la distribution cible. L’exemple ci-haut peut théoriquement s’appliquer à n’importe
quelle distribution cible donnée, bien qu’il se comportera mieux dans certaines condi-
tions que dans d’autres, et son efficacité peut généralement être améliorée en ajustant la
variance σ 2 . Une autre forme usitée pour la distribution instrumentale est la distribution
uniforme, par exemple Uni f orme(x − 1, x + 1), mais d’autres choix sont aussi possibles.
(b) Algorithme symétrique de type Metropolis-Hastings
Certains algorithmes, comme la méthode RWMH avec loi instrumentale N(x, σ 2 )
et bien d’autres, ont une propriété additionnelle conférée par la symétrie de la densité
instrumentale, notamment le fait que q(x, y) = q(y, x). Cette caractéristique permet de
simplifier le calcul de α :

πu (y)
α(x, y) = min 1, .
πu (x)
(c) Échantillonneur Metropolis-Hastings indépendant (Independent Sampler (IS))

La méthode IS est une approche où la distribution instrumentale est indépendante
de l’état présent, c’est-à-dire que q(x, y) = q(y). On peut choisir, par exemple, une dis-
tribution instrumentale Student f , où f représente les degrés de liberté. Étant donné un
état actuel x, un nouvel état y est proposé selon une loi Student f (y) indépendamment
de x. Cette nouvelle valeur est ensuite acceptée avec probabilité (1.4). Cette méthode
requiert en général une connaissance plus approfondie de la distribution cible pour une
convergence optimale. Par exemple, il est préférable que les régions de forte densité de
la distribution instrumentale coïncident avec les régions de forte densité de la distribu-
tion cible. D’un autre côté, il est souhaitable que l’épaisseur des queues des deux densités
soit relativement similaire afin de ne pas sous-estimer certaines régions, particulièrement
6
dans le calcul de valeurs-p.

(d) Échantillonneur de Gibbs (Gibbs Sampler (GS))
Dans certaines situations, il peut être plus simple de générer des observations à partir
de densités conditionnelles plutôt que de densités conjointes. L’approche GS se base
sur cette idée et peut être un algorithme très puissant pourvu que l’on puisse aisément
exprimer les densités conditionnelles.
On suppose que πu (·) est une densité n-dimensionnelle sur un sous-espace
X de Rn et on exprime les points de la façon x = (x1 , . . . , xn ). Pour un
(1)
état présent x(0) , la méthode GS génère une observation x1 pour la première
(0) (0) (0) (0) (1)
composante selon π(x1 |x2 , x3 , x4 , . . . , xn ), ensuite une observation x2 selon
(1) (0) (0) (0) (1)
π(x2 |x1 , x3 , x4 , . . . , xn ) et ainsi de suite jusqu’à l’observation xn . Cette méthode
génère une chaîne de Markov (voir par exemple Gelman et al. (1995)) qui aura comme
loi stationnaire π(·). En effet, les probabilités de transition satisfont (1.1) :
Z
π(dx)P(x, dy)
x∈Rn Z Z Z
= ... πu (x1 , x2 , x3 , . . . , xn ) ∗ πu (y1 |x2 , x3 , . . . , xn ) ∗ πu (y2 |y1 , x3 , . . . , xn ) ∗ . . . ∗
xn xn−1 x1
πu (yi |y1 , y2 , . . . , yi−1 , xi+1 , . . . , xn ) ∗ . . . ∗ πu (yn |y1 , y2 , . . . , yn−1 )dx1 dx2 . . . dxn dy
Z Z Z
= ... πu (x2 , x3 , . . . , xn ) ∗ πu (y1 |x2 , x3 , . . . , xn ) ∗ πu (y2 |y1 , x3 , . . . , xn ) ∗ . . . ∗
xn xn−1 x2
Z Z Z
= ... πu (y1 , x2 , x3 , . . . , xn ) ∗ πu (y2 |y1 , x3 , . . . , xn ) ∗ . . . ∗
xn xn−1 x2
7
Z Z Z
= ... πu (y1 , x3 , . . . , xn ) ∗ πu (y2 |y1 , x3 , . . . , xn ) ∗ . . . ∗
xn xn−1 x3
Z Z Z
= ... πu (y1 , y2 , x3 , . . . , xn ) ∗ πu (y3 |y1 , y2 , x4 . . . , xn ) ∗ . . . ∗
xn xn−1 x3
= ...
Z
= πu (yn |y1 , y2 , . . . , yn−1 ) πu (y1 , y2 , . . . , yn−1 , xn )dxn dy
xn
= π(y1 , . . . , yn )dy = π(y)dy.
L’algorithme GS peut s’implémenter de façon déterministe comme l’exemple ci-haut

ou bien de façon aléatoire, c’est-à-dire que la composante mise à jour à chaque itération
est choisie au hasard parmi les n composantes. En général, l’échantillonneur de Gibbs
aléatoire est réversible, mais celui déterministe ne l’est pas (Liu et al. (1995)).
1.3 Convergence des algorithmes MCMC
Il est maintenant clair qu’il existe plusieurs façons simples de construire des chaînes
de Markov avec distribution stationnaire π(·). Toutefois, même si la distribution sta-
tionnaire est connue, elle n’est pas nécessairement unique et la chaîne n’y converge pas
dans tous les cas. À titre illustratif, deux exemples tirés de Roberts et Rosenthal (2004)
démontrent la nécessité de l’irréductibilité et de l’apériodicité de la chaîne.
Exemple : Supposons que X = {1, 2, 3} et que π{1} = π{2} = π{3} = 1/3. Aussi,
P(1, {1}) = P(1, {2}) = P(2, {1}) = P(2, {2}) = 1/2 et P(3, {3}) = 1. Ici, P(x, {y}) =
P(X j+1 = y|X j = x). La distribution stationnaire pour cette chaîne est π(·), mais si X0 est
1, alors X j ∈ {1, 2} pour tout j et donc P(X j = 3) ne converge pas vers π({3}).
Cet exemple illustre une chaîne de Markov dite réductible, puisque l’état 3 ne peut
jamais être atteint à partir de l’état 1 ou 2. Dans ce cas, la distribution stationnaire n’est
pas unique et dépend du point de départ de la chaîne. Lorsque X est dénombrable, la
8
caractéristique d’irréductibilité signifie que chaque état a une probabilité non nulle d’être
éventuellement atteint à partir de n’importe quel autre état. Dans un contexte où X est
non dénombrable, une condition similaire, quoique moins robuste, peut être définie.
Définition 2. Une chaîne de Markov est φ -irréductible s’il existe une mesure φ , non-
nulle et σ -finie, sur l’espace X telle que pour tout A ⊆ X avec φ (A) > 0 et pour tout
x ∈ X, il existe j ∈ N tel que P j (x, A) > 0, où P j (x, A) = P(X j ∈ A|X0 = x).
Exemple : Supposons encore que X = {1, 2, 3} et que π{1} = π{2} = π{3} = 1/3.
Aussi, P(1, {2}) = P(2, {3}) = P(3, {1}) = 1. La distribution stationnaire est encore π
et la chaîne est irréductible. Toutefois, si X0 = 1 alors P(X j = 1) oscillera entre 0 et 1 à
tous les multiples de 3, et donc ne convergera pas vers 1/3.
Ce deuxième exemple démontre le concept de périodicité, c’est-à-dire qu’un état
n’est atteignable que dans un nombre précis de pas ou tout multiple de ce nombre.
La condition d’irréductibilité (ou de φ -irréductibilité) n’est donc pas suffisante et une
deuxième condition, celle de l’apériodicité, est requise.
Définition 3. Une chaîne de Markov est apériodique s’il n’existe pas de sous-ensembles
disjoints X1 , X2 , . . . , Xd ⊆ X avec d > 1 et π(Xi ) > 0 pour tout i, tel que P(x, Xi+1 ) = 1
pour tout x ∈ Xi (1 ≤ i ≤ d − 1) et P(x, X1 ) = 1 pour tout x ∈ Xd .
Ensemble, ces deux conditions garantissent la convergence d’une chaîne de Markov

vers la distribution stationnaire (unique) π(·).
Théorème 1. Une chaîne de Markov φ -irréductible et apériodique avec distribution

stationnaire π sur un espace X avec σ -algèbre dénombrablement générée, converge
asymptotiquement vers π :
lim ||P j (x, ·) − π(·)|| = 0,

j→∞
pour π-presque partout x ∈ X. De plus, lim j→∞ P j (x, A) = π(A) pour tout A mesurable.
9
Ici, || · || dénote la distance de variation totale entre deux mesures quelconques

µ1 (·), µ2 (·) :
||µ1 (·) − µ2 (·)|| = supA⊆X |µ1 (A) − µ2 (A)|.
Théorème 2. Sous les mêmes conditions que le théorème précédent, pour toute fonction
h : X → R une forme de loi forte des grands nombres existe
N
1
lim
N→∞ N
∑ h(X j ) = Eπ (h(x)) ≡ π(h)
j=1
presque toujours si π(|h|) < ∞.
La démonstration de ces théorèmes est assez longue et sera omise dans ce mémoire,
mais elle se retrouve dans son entièreté dans Roberts et Rosenthal (2004). Les condi-
tions de ces théorèmes sont souvent satisfaites dans le cas des algorithmes MCMC. Par
construction, la réversibilité des chaînes MCMC garantit l’existence d’une distribution
stationnaire π. Les conditions de φ -irréductibilité et d’apériodicité doivent être vérifiées,
mais généralement elles seront satisfaites dans presque tous les cas pratiques.
1.4 Taux de convergence des algorithmes MCMC
Bien que le théorème 1 garantisse la convergence asymptotique des chaînes MCMC,

il n’y a aucune indication quant à la vitesse de cette convergence. Lorsque cela est pos-
sible, il est souhaitable d’obtenir des bornes analytiques pour la vitesse de convergence
d’un algorithme. Il existe certaines propriétés qui peuvent s’avérer utiles dans un nombre
limité de cas.
Définition 4. Une chaîne de Markov avec distribution stationnaire π(·) est dite unifor-
mément ergodique s’il existent 0 < ρ < 1 et M < ∞ tels que
||P j (x, ·) − π(·)|| ≤ Mρ j , j = 1, 2, 3, . . . ∀x ∈ X.

10
Définition 5. Un sous-ensemble C ⊆ X est dit petit s’il existe un entier positif n0 ∈ N,

un réel ε > 0 et une mesure de probabilité ν(·) sur X tels que
Pn0 (x, ·) ≥ εν(·) x ∈ C
et donc, Pn0 (x, A) ≥ εν(A) pour tout x ∈ C.
De façon intuitive, cette définition implique que les probabilités de transition en n0

pas à partir de C ont une composante ε en commun. Un exemple tiré de Roberts et
Rosenthal (2004) peut permettre de clarifier cette définition.
−(y−x)2
On suppose que la densité de transition q(x, y) ∝ e 2 et que

 |x|− 12 , |x| < 1
πu (x) = .
 0, sinon
Les probabilités de transition satisfont

πu (y)
P(x, dy) = q(x, y) dy min{1, }
πu (x)
Il est facile de voir qu’un voisinage de 0 n’est pas un ensemble petit, puisque πu (x)
n’est pas bornée et la probabilité de transition peut être arbitrairement près de 0. D’un
autre côté, un ensemble compact C où πu (x) < k < ∞ pour tout x ∈ C est petit. En effet,
si D est un autre ensemble compact avec mesure de Lebesgue et π positives et que
infx∈C,y∈D q(x, y) = ε > 0 alors
πu (y) πu (y)
P(x, dy) ≥ ε dy min{1, } ≥ ε dy min{1, }.
πu (x) k
Cette dernière est une mesure positive indépendante de x et donc C est bel et bien
petit.
Théorème 3. Si une chaîne de Markov possède une distribution stationnaire π(·) et

l’ensemble d’états X est un ensemble petit pour n0 ∈ N, ε > 0 et ν(·), alors la chaîne
est uniformément ergodique et ||P j (x, ·) − π(·)|| ≤ (1 − ε)b j/n0 c pour tout x ∈ X (ici, b·c
représente la fonction partie entière).
11
Démonstration. Voir Roberts et Rosenthal (2004).
Ce dernier théorème est important puisqu’il permet d’obtenir des bornes analytiques
pour la distance de variation totale entre P j (x, ·) et π(·) et donc pour le taux de conver-
gence de l’algorithme. En sachant n0 et ε, on peut trouver j tel que, par exemple,
||P j (x, ·) − π(·)|| ≤ 0, 01. Dans ce cas, on sera certain que la différence entre la dis-
tribution stationnaire et la distribution des états à la j-ième itération sera d’au plus 0,01
(ou tout autre nombre choisi par l’utilisateur). Ce théorème est applicable dans certaines
situations particulières, comme par exemple au chapitre 2 du présent mémoire.
Il existe une deuxième propriété relative à la vitesse de convergence qui peut aussi
s’avérer utile.
Définition 6. Une chaîne de Markov avec distribution stationnaire π(·) est dite géomé-
triquement ergodique s’il existent 0 < ρ < 1 et M(x) < ∞ pour π-presque partout x ∈ X
tels que
||P j (x, ·) − π(·)|| ≤ M(x)ρ j , j = 1, 2, 3, . . .
Dans ce mémoire, la propriété de l’ergodicité géométrique restera qualitative, mais

il est possible d’arriver à des bornes quantitatives qui peuvent être utiles dans certains
contextes, quoiqu’assez difficilement applicables (voir Roberts et Rosenthal (2004)).
1.5 Théorème central limite des algorithmes MCMC
Le théorème central limite pour les méthodes MCMC est vastement applicable et
peut permettre de comprendre et de quantifier les erreurs Monte Carlo. En général, l’in-
térêt d’une méthode MCMC est non seulement d’obtenir un échantillon de la distribution
cible mais aussi d’en estimer certains paramètres comme par exemple : la moyenne, la
variance, les quantiles et ainsi de suite. Donc, on construira une fonction des valeurs x j
générées et on sera intéressé à la comparer avec sa valeur théorique. Plus précisément, si
une chaîne de Markov possède une loi stationnaire π(·) et h : X → R est une fonction avec
12
1
∑Nj=1 h(X j ). Il
R
π(h) = x∈X h(x)π(dx) < ∞, on s’intéressera à l’estimateur naturel h = N
faut rappeler que cet estimateur est justifiable par la loi forte des grands nombres pour
les chaînes de Markov (théorème 2).
Définition 7. Étant donnée une fonction h : X → R et une chaîne de Markov avec distri-
bution stationnaire π(·), la fonction h satisfait un théorème central limite si
N
1
lim √
N→∞ N
∑ [h(Xi) − π(h)] ∼ N(0, σ 2) et σ 2 < ∞.
j=1
Il est immédiatement clair qu’une telle propriété est intéressante dans un contexte
MCMC, car elle permet de construire des intervalles de confiance arbitrairement précis
pour un certain paramètre. De plus, elle consent une certaine quantification objective des
erreurs Monte Carlo permettant au lecteur indépendant de tirer des conclusions auto-
nomes.
La seule difficulté apparaît lorsque l’on se rend compte de la complexité de la va-
riance de la distribution limite qui est due à la corrélation inhérente de la chaîne de
Markov. En général, σ 2 6= Varπ (h) et donc la variance empirique de h(X j ), qui est un
estimateur naturel, n’est pas applicable. Toutefois, il existe plusieurs autres façons d’es-
timer σ 2 , mais pour le moment la forme de cette variance est présentée.
Théorème 4. Pour une chaîne de Markov avec distribution stationnaire π(·), si une
fonction h : X → R satisfait un théorème central limite, alors la variance σ 2 de la dis-
tribution limite normale est donnée par :
 !2 
N
1
σ 2 = lim E  ∑ [h(X j ) − π(h)] 
N→∞ N
j=1
ou de façon équivalente par σ 2 = τVarπ (h), où τ = ∑k∈Z Corr(X0 , Xk ).
Une démonstration de ce théorème peut être trouvée dans Chan et Geyer (1994).
La quantité τ, aussi appelée le temps d’autocorrélation intégrée, est une mesure de la
corrélation de la chaîne. De par la deuxième expression pour σ 2 , il est clair que la condi-
tion σ 2 < ∞ de la définition 7 est satisfaite lorsque non-seulement Varπ (h) < ∞ (ou
13
π(h2 ) < ∞), mais aussi τ < ∞. En effet, même si Varπ (h) < ∞, il peut arriver qu’une
chaîne MCMC est tellement sous-optimale qu’elle reste confinée à une seule région et
τ = ∞ ( voir par exemple Roberts (1999)).
L’utilité du théorème central limite pour les chaînes MCMC est maintenant claire,
mais les conditions garantissant cette propriété n’ont pas encore été mentionnées. Les
théorèmes suivants présentent des conditions suffisantes pour assurer un théorème cen-
tral limite pour une fonction h et les démonstrations sont disponibles dans Kipnis et
Varadhan (1986); Cogburn (1972); Roberts et Rosenthal (1997).
Théorème 5. Pour une chaîne de Markov réversible, φ -irréductible et apériodique, un

théorème limite central existe pour une fonction h si σ 2 < ∞.
Les conditions de ce théorème sont similaires à celles qui garantissent la convergence

d’un algorithme MCMC (théorème 1), donc elles seront satisfaites la plupart du temps.
Cependant, il existe des algorithmes MCMC, comme la version déterministe de l’échan-
tillonneur de Gibbs, qui sont convergents même si non réversibles. En plus, l’existence
d’un théorème central limite repose sur le fait que σ 2 < ∞, ce qui nécessite la vérification
de π(h2 ) < ∞ et τ < ∞. Toutefois, si le taux de convergence de l’algorithme est connu,
il est possible de restreindre les conditions sur π(h2 ) seulement.
Théorème 6. Pour une chaîne de Markov avec distribution stationnaire π(·) et unifor-
mément ergodique, un théorème central limite existe pour une fonction h si π(h2 ) < ∞.
Théorème 7. Pour une chaîne de Markov réversible et géométriquement ergodique, un

théorème central limite existe pour une fonction h si π(h2 ) < ∞.
Il semble donc que le théorème central limite existe pour une quantité considérable de
méthodes MCMC. Cependant, même si cette convergence est garantie, il reste à trouver
un estimateur pour la variance de la distribution limite. Comme mentionné, l’estimateur
usuel par la variance empirique est inadéquat à cause de la corrélation entre les états
d’une chaîne de Markov. Il existe plusieurs techniques pour estimer σ 2 , mais dans ce
14
mémoire il sera question d’une des méthodes les plus populaires, soit celle de la moyenne
par séries (Batch Means). Cette technique est à la fois simple et versatile. Supposons que
nous disposons d’un échantillon de taille N généré par un algorithme MCMC. L’idée est
de séparer cet échantillon en a séries de taille b, tel que N = ab. On définit ensuite
1 lb−1
Yl = ∑ h(xi)
b i=(l−1)b
pour l = 1, . . . , a et l’estimateur de la variance σ 2 par

!
2 1 1 a 2 b a
σ̂MS = N ∑ (Y l − hN ) = ∑ (Y l − hN )2
a a − 1 l=1 a − 1 l=1
où hN représente 1a ∑al=1 (Yl ) = 1

N ∑Nj=1 (h(x j )).
Cette méthode ne produit pas un estimateur cohérent en général. Cependant, sous
certaines conditions relativement souples établies par Damerdji (1994), cet estimateur
2 calculé, un intervalle de confiance
peut converger même presque toujours. Une fois σ̂MS
√
peut être généré de la manière habituelle avec demi-largeur ta−1 σ̂MS / N, où ta−1 est
le quantile associé au niveau de confiance désiré d’une distribution Student avec a −
1 degrés de liberté. Il est à noter qu’il existe différentes versions de la méthode par
moyenne de séries, notamment celle par moyenne de séries chevauchées (Overlapping
Batch Means), avec des propriétés similaires.
Ces méthodes supposent que les moyennes de chaque série sont approximativement
indépendantes et identiquement distribuées et donc l’estimateur empirique de la variance
de ces moyennes est justifié. Une autre approche est la méthode de Geyer dont l’idée
est d’estimer τ et d’utiliser cette estimation en conjonction avec la variance empirique
de h(x j ) pour obtenir σ̂ 2 selon le théorème 4. Pour cette technique et d’autres méthodes
spectrales, le lecteur est invité à consulter les références suivantes Flegal et Jones (2010);
Geyer (1992).
15
1.6 Diagnostics de convergence
En pratique, il est difficile de prouver l’ergodicité uniforme ou géométrique d’un

algorithme MCMC. En effet, dans la plupart des contextes réalistes il est impossible
ou très compliqué de connaître la vitesse de convergence de la méthode utilisée afin de
savoir quand terminer la simulation. Pour cette raison, plusieurs techniques d’analyse
de sorties MCMC ont été développées afin de juger de la performance d’un algorithme
donné. Ces méthodes s’appellent collectivement « diagnostics de convergence » et il en
existe un nombre considérable. Il est à noter que leur application nécessite une certaine
expérience de l’utilisateur. En plus, elles ne donnent aucune garantie rigoureuse quant à
la convergence et peuvent parfois induire l’utilisateur en erreur (voir Matthews (1993)).
Néanmoins, elles constituent une partie importante de l’implémentation de l’approche
MCMC et sont souvent la première façon de valider un algorithme donné. En effet, tel
que mentionné à la section 1.2, il existe plusieurs méthodes MCMC et elles auront une
performance différente d’un contexte à l’autre. L’utilisation de diagnostics de conver-
gence peut éliminer les algorithmes inefficaces et orienter l’utilisateur vers les méthodes
appropriées. Il existe différentes mesures ou diagnostics et ce mémoire en présentera
quelques-uns des plus couramment utilisés.
1.6.1 Graphique de la trace d’un paramètre d’intérêt
Une façon simple de vérifier la convergence d’un algorithme MCMC est de produire
le graphique de la trace d’un paramètre d’intérêt, c’est-à-dire sa valeur par rapport à
l’itération. Si la méthode a convergé vers la distribution stationnaire, le paramètre devra
varier aléatoirement dans le temps et le graphique devra ressembler à du bruit de fond. Si
le graphique présente des tendances quelconques ou une stagnation pendant de longues
périodes, cela est une forte indication que la distribution stationnaire n’est pas atteinte.
Par exemple, à la figure 1.1, le graphique de la trace de gauche ne montre à priori aucun
problème de convergence, tandis que celui de droite indique un algorithme sous-optimal.
16
Figure 1.1 – Graphique de la trace de paramètres en fonction des itérations. À gauche un

paramètre qui semble fluctuer aléatoirement. À droite un paramètre qui stagne.
1.6.2 Densité (histogramme) de valeurs générées et chaînes parallèles
L’histogramme des valeurs générées par un algorithme MCMC peut donner une idée
de la convergence de l’algorithme. Généralement, une densité empirique avec plusieurs
modes indique un possible problème de convergence. Aussi, il est possible de comparer
des chaînes parallèles de mêmes dimensions mais avec valeurs initiales différentes. Si
l’algorithme est optimal, la densité empirique (ou bien un diagramme en boîte) ne devrait
pas dépendre de la valeur initiale. La figure 1.2, présente à gauche deux chaînes paral-
lèles qui convergent vers la même distribution et à droite un algorithme sous-optimal.
1.6.3 Taux d’acceptation
La nature de l’approche MCMC requiert une étape d’acceptation ou de refus d’une

valeur générée de la distribution instrumentale afin de garantir la réversibilité du proces-
sus. Le taux d’acceptation est une mesure facilement extraite de l’implémentation infor-
matique d’un algorithme donné. En général, un taux d’acceptation trop élevé (> 85%)
17
Figure 1.2 – Graphique en boîte de valeurs générées par deux chaînes parallèles avec
valeur initiale différente. À gauche les deux chaînes semblent converger vers la même
distribution. À droite elles convergent vers deux distributions distinctes.
indique une acceptation de sauts qui sont souvent petits et donc une mauvaise explo-
ration générale de l’espace. D’un autre côté, un taux d’acceptation trop petit (< 15%)
signifie le rejet de la plupart des sauts et donc encore une exploration de l’espace in-
adéquate. Le chapitre 3 abordera ces concepts et quelques lignes directrices générales
quant au taux d’acceptation optimal en plus de détails. Évidemment, ce taux dépend
souvent du contexte, mais généralement des valeurs extrêmes peuvent être indicatrices
d’un problème de convergence.
1.6.4 Autocorrélation et distance de saut carrée moyenne
Tel que mentionné plus tôt, les valeurs générées par un algorithme MCMC possèdent
une corrélation inhérente puisqu’elles sont dépendantes au moins par paires. Une ma-
nière de quantifier cette corrélation est d’utiliser un estimateur empirique de l’autocor-
rélation de la chaîne.
Définition 8. L’autocorrélation entre deux temps j et l d’un processus stochastique

chronologique est définie par

E (X j − µ j )(Xl − µl )
Acorr( j, l) = ,
σ j σl
où µ j , µl , σ j , σl sont les moyennes et écarts-types aux temps j et l.
18
Dans un contexte MCMC, les moyennes et variances sont invariables dans le temps
lorsque la chaîne a convergé vers la distribution stationnaire. Donc, l’autocorrélation
pour un processus MCMC peut s’exprimer à travers un laps de temps (ou intervalle) k
qui sépare deux temps donnés :

E (X j − µ)(X j+k − µ)
Acorr(k) = .
σ2
Un estimateur empirique naturel de l’autocorrélation d’une chaîne MCMC peut donc
être défini par
∑N−k
j=1 (x j − x)(x j+k − x)
ρ̂k = Corr(X
d j , X j+k ) = .
∑Nj=1 (x j − x)2
Même s’il y a convergence, il est normal que l’autocorrélation soit forte pour un in-
tervalle k petit, mais elle devrait diminuer fortement au fur et à mesure que k augmente.
Une autocorrélation qui ne s’estompe pas signifie généralement une convergence sous-
optimale. À titre d’exemple, à la figure 1.3, le corrélogramme de gauche n’indique pas à
priori un problème de convergence, mais celui de droite indique une situation probléma-
tique.
Figure 1.3 – Graphique de l’autocorrélation de valeurs générées. À gauche l’autocor-

rélation diminue de façon satisfaisante. À droite elle demeure excessive même pour un
grand intervalle.
La distance de saut carrée moyenne (DSCM) est une autre mesure de la performance
d’un algorithme MCMC. Comme le nom l’indique, elle mesure la moyenne de la dis-
tance carrée entre deux états consécutifs. Elle est souvent utilisée comme mesure de
19
performance lorsque plusieurs méthodes s’offrent à l’expérimentateur ou lorsqu’un pa-

ramètre de la distribution instrumentale est ajusté. Intuitivement, une grande DSCM in-
dique un algorithme qui explore bien l’espace et est préférable à une DSCM inférieure,
toutes autres mesures égales par ailleurs.
Définition 9. Étant une chaîne de Markov avec distribution stationnaire π(·). L’espé-
rance de la distance de saut carrée (EDSC) entre deux états consécutifs j et j + 1 est
définie par
EDSC = E (X j − X j+1 )2 .

(1.5)
La distance de saut carrée moyenne (DSCM) est l’estimateur naturel empirique de

(1.5). Pour un algorithme MCMC avec N itérations, la DSCM est définie par
1 N−1
DSCM = ∑ (x j − x j+1)2.
N j=0
Enfin, il est à noter qu’il existe une relation entre l’autocorrélation et l’espérance de
la distance de saut carrée pour un paramètre donné lorsque X j , X j+1 ∼ π(·).
E (X j − X j+1 )2 = E X j2 − 2E X j X j+1 + E X j+1

2
= σ 2 + µ 2 − 2E X j X j+1 + σ 2 + µ 2

= 2σ 2 − 2 E X j X j+1 − µ 2

= 2σ 2 − 2Cov(X j , X j+1 )
= 2σ 2 (1 − Acorr(1)). (1.6)
Donc la DSCM est intimement reliée à l’autocorrélation d’intervalle 1, justifiant son

utilité et son mérite en pratique. Dans le contexte multidimensionnel, il est possible de
définir la DSCM comme la distance euclidienne carrée entre deux états consécutifs et
dans ce cas, elle est équivalente à la somme des DSCM de chaque composante.
20
1.6.5 Statistiques de convergence
Plusieurs tests statistiques de sorties MCMC ont été développés afin d’obtenir un
niveau de confiance plus objectif pour la convergence. Les plus courants sont les tests de
Gelman-Rubin (Gelman et Rubin (1992)), de Geweke (Geweke (1992)) et de Raftery-
Lewis (Raftery et Lewis (1992)). Le test de Gelman-Rubin est analogue à une analyse
de variance, c’est-à-dire qu’il examine l’écart entre la variance dans les chaînes et la
variance entre les chaînes ayant des valeurs de départ différentes. Le test de Geweke
se base sur l’écart entre la moyenne des premières Na observations et la moyenne des
Nb observations subséquentes d’une seule chaîne telle que N = Na + Nb . Finalement,
le test de Raftery-Lewis indique le nombre d’itérations requis pour estimer les quantiles
d’une fonction des valeurs générées ainsi que la dépendance entre ces valeurs. La théorie
motivant ces tests ne fait pas l’objet de ce mémoire, mais le lecteur est invité à consulter
les références mentionnées pour de plus amples détails.
Les diagnostics de convergence présentés s’utilisent en pratique de façon complé-
mentaire. Cela dit, il existe de nombreux exemples où un problème de convergence ne
peut être détecté par ces méthodes Cowles et Carlin (1996). Néanmoins, les diagnostics
de convergence restent la première tentative à postériori afin de juger de la convergence
d’un algorithme. Ces méthodes peuvent être incorporées manuellement dans l’implé-
mentation d’une méthode donnée ou bien utilisées directement sur les sorties MCMC à
l’aide de certains logiciels (ex : le greffon CODA pour le logiciel R).
1.7 Variations d’algorithmes MCMC
Il existe plusieurs autres approches MCMC qui sont une variation des algorithmes
élémentaires décrits plus tôt. L’objectif de ces méthodes est principalement d’améliorer
la convergence des algorithmes existants. Ces techniques sont souvent plus complexes
et nécessitent de plus amples efforts pour leur implémentation. Cependant, elles offrent
une alternative intéressante dans plusieurs contextes où les algorithmes traditionnels se
21
comportent de façon médiocre. Le chapitre 2 de ce mémoire présentera un algorithme

avec ajustement directionnel récemment développé dans Bédard et Fraser (2008) qui est
une variante de l’échantillonneur indépendant. Ensuite, le chapitre 3 présentera un algo-
rithme adaptatif, développé par Haario et al. (2001), qui met à jour certains paramètres
de la distribution instrumentale basé sur les valeurs générées antérieurement. En plus, ce
chapitre étudiera d’autres variations comme l’algorithme à essais multiples (Multiple-
Try MCMC) et l’algorithme avec rejet différé (Delayed Rejection MCMC). Les proprié-
tés théoriques ainsi que les motivations derrière ces méthodes seront abordées dans ces
sections.
CHAPITRE 2
ALGORITHME DE TYPE METROPOLIS-HASTINGS AVEC AJUSTEMENT

DIRECTIONNEL
L’algorithme de type Metropolis-Hastings avec marche aléatoire (RWMH) est une

des méthodes les plus versatiles et les plus couramment utilisées. La fonction de densité
instrumentale est de façon répandue une fonction de densité normale ou uniforme, ce
qui requiert seulement une optimisation de la variance ou du support respectivement.
Cependant, l’application de la méthode RWMH dans presque n’importe quel contexte
avec sensiblement les mêmes stratégies d’optimisation résulte en une convergence qui
peut être extrêmement lente selon la situation.
L’échantillonneur indépendant (IS) est une méthode dont la densité instrumentale est
indépendante de l’état présent. L’application de cette approche nécessite une connais-
sance plus approfondie de la distribution cible afin de bénéficier d’une convergence
raisonnable. En général, plus la densité instrumentale se rapproche de la densité cible
plus l’algorithme converge rapidement, mais moins il est facile de générer des valeurs
proposées. Cependant, avec un choix adéquat de la densité instrumentale, nécessitant
souvent considérablement plus d’efforts, la convergence de l’algorithme IS est souvent
supérieure à celle de l’algorithme RWMH.
L’échantillonneur indépendant avec ajustement directionnel (DA)(Bédard et Fraser
(2008)) est une nouvelle approche qui vise à combiner la versatilité de l’algorithme
RWMH avec la haute performance de l’algorithme IS pour des densités lisses et uni-
modales. Cette méthode utilise une densité instrumentale Student centrée au mode de la
densité cible et dont les queues sont ajustées de façon directionnelle afin d’approximer le
plus possible la densité cible. En un premier lieu, ce chapitre présentera une introduction
détaillée de la méthode DA. Ensuite, la deuxième section sera consacrée à l’ajustement
d’un paramètre clé de l’algorithme et finalement, la dernière partie portera sur la conver-
23
gence de cette approche.
2.1 L’algorithme DA
Étant donnée une densité cible π(·) n-dimensionnelle, unimodale et lisse et en

supposant qu’il est possible d’en trouver le mode, il peut être avantageux de tenir
compte de cette information dans le choix d’une densité instrumentale. Plus précisé-
ment, il serait préférable pour la densité instrumentale en question d’être centrée au
mode x̂ = arg supx π (x) et d’avoir la même courbure que la densité cible au mode. La
matrice hessienne évalue la courbure de π à un point particulier et est donnée par
 
∂ 2π ∂ 2π ∂ 2π
 ∂ x12 ∂ x1 ∂ x2 · · · ∂ x1 ∂ xn 
 2
∂ 2π ∂ 2π 

 ∂ π · · ·
 ∂ x2 ∂ x1 ∂ x22 ∂ x2 ∂ xn 
.. .. .. .
...

. . .
 
 
 2 2 2

∂ π ∂ π ∂ π
∂ xn ∂ x ∂ xn ∂ x · · · ∂ x2
1 2 n
Cette matrice ne sera typiquement pas définie positive au mode. Par conséquent, il
est préférable d’utiliser une autre forme familière, Ĥ = −∂ 2 log π (x) /∂ x∂ x0 |x=x̂ , c’est-
à-dire le négatif de la matrice hessienne du logarithme de la densité évaluée au mode.
À ce point, il serait possible d’implémenter un algorithme IS avec une densité instru-
mentale normale ajustée pour coïncider avec x̂ et Ĥ −1 . Il est à noter que Ĥ = Σ−1 dans
le cas d’une densité normale multivariée. En effet, si y ∼ N(µ, Σ), alors
∂2
(−n/2) 1
( − 12 (y−µ)0 Σ−1 (y−µ))

Ĥ = − log (2π) |Σ| e
2
∂ y∂ y0
∂2

1 0 −1
= (y − µ) Σ (y − µ)
∂ y∂ y0 2
1 ∂
= (Σ−1 (y − µ) + Σ−10 (y − µ))
2 ∂ y0
1 −10
= (Σ + Σ−1 )0
2
1 −1
= 2Σ = Σ−1 . (2.1)
2
24
Cette approche, par la nature des queues courtes de la densité normale, peut sous-
estimer les régions se trouvant dans les queues de la densité cible, si ces dernières
sont plus épaisses. À l’inverse, l’utilisation d’une distribution instrumentale Cauchy, par
exemple, peut résulter inutilement en une trop grande proportion de sauts rejetés si les
queues de la densité cible sont plus minces. Afin de remédier à ces problèmes potentiels,
l’approche DA vise à construire une densité instrumentale qui n’est pas rigide et qui est
adaptée à la densité cible.
Pour ce faire, la densité instrumentale Student f (0, In ) est utilisée, dont la forme ca-
nonique est :

f +n
Γ 2 − f +n
q f (t) = 1 + t12 + . . . + tn2 2
π n/2 Γ 2f

f +n
Γ 2 f +n
= 1 + t0 t − 2 , (2.2)
π n/2 Γ 2f
où t0 = (t1 , . . . ,tn ).
Aussi, il est à noter qu’il est facile de générer une variable aléatoire t0 à partir de va-
riables indépendantes z1 , . . . , zn distribuées selon une loi normale standard et une variable
aléatoire χ 2f distribuée selon une loi chi-carrée avec f degrés de liberté et indépendante
des zi :

0 z1 zn
t = ,..., . (2.3)
χf χf
La matrice hessienne du négatif du log de la densité Student canonique à t = 0 est

égale à ( f + n) In et en général, par un exercice similaire à (2.1), Ĥ = ( f + n)Σ−1 au
mode d’une densité Student f (µ, Σ). Afin de faire correspondre le mode et la matrice
hessienne à ceux de la densité cible, une possibilité est d’utiliser une densité instrumen-
tale Student f x̂, ( f + n) Ĥ −1 avec valeurs proposées (y) générées de la façon suivante :

y = x̂ + ( f + n)1/2 Ĥ −1/2 t,
25
où Ĥ 1/2 est la matrice racine carrée de droite telle que Ĥ = (Ĥ 1/2 )0 (Ĥ 1/2 ) et t est donné
par (2.3).
Dans ce cas, la forme de la densité des valeurs proposées y est alors donnée par

f +n − f +n
Γ 2
(y − x̂)0 Ĥ (y − x̂) 2 Ĥ 1/2
q f (y) = 1+ .
π Γ f
n/2
2
f +n ( f + n)n/2
Cependant, d’une façon équivalente et plus simple, il est aussi possible de reparamé-
triser les densités cible et instrumentale en posant x∗ = Ĥ 1/2 (x − x̂) et T = ( f + n)1/2 t.
Cela implique que x̂∗ = 0, Ĥ ∗ = In et T̂ = 0, Ĥ = In . De plus, la densité cible repara-
métrisée est proportionnelle à l’originale, ce qui signifie que le ratio d’acceptation (1.4)
et l’approximation directionnelle (2.6) qui sera détaillée sous peu seront calculables à
partir de π(·) :

π ∗ (x∗ ) = π x̂ + Ĥ −1/2 x∗ Ĥ −1/2 ∝ π x̂ + Ĥ −1/2 x∗ = π (x) .
En outre, il est possible à tout moment d’accéder à une variable x à partir de x∗ par
la transformation inverse : x = x̂ + Ĥ −1/2 x∗ .
Cette transformation implique que la nouvelle densité instrumentale est donnée par

Γ f +n
2
f +n
T0 T − 2
q f (T) = 1+ ( f + n)−n/2 . (2.4)
n/2
π Γ f f + n
2
Donc, plusieurs choix existent quant à la reparamétrisation des densités cible et ins-
trumentale. Celle qui a été décrite à l’instant sera utilisée dans la description ainsi que
dans l’implémentation de l’algorithme DA tout au long de ce mémoire.
Il serait maintenant possible de spécifier un degré de liberté global pour la densité
instrumentale et d’utiliser un algorithme IS. Cependant, l’approche DA vise justement
à en faire davantage en optimisant la densité instrumentale dans chaque direction. Pour
ce faire, il est intéressant d’exprimer les états de la chaîne de Markov en termes de deux
composantes : la direction et la distance par rapport au mode. Ainsi, la direction d’un
26
q
état présent x j est donnée par u j = x j / x j où |x| = x12 + . . . + xn2 et sa distance radiale
par s j = x j . Donc, en posant x j = u j · s j , il est possible d’exprimer la densité cible en

terme de la densité de chaque composante, c’est-à-dire π̃ u j , s j = π̃ s j u j π̃ u j . De
façon similaire, une valeur proposée peut également être exprimée en termes de direction
et distance, y j+1 = u prop prop
j+1 · s j+1 et la densité instrumentale correspondante devient

prop prop prop prop prop
q̃ u j+1 , s j+1 = q̃ s j+1 u j+1 q̃ u j+1 .
Le but de cette technique est de diviser l’étape de proposition d’un nouvel état
en deux parties afin de permettre un ajustement directionnel. Donc, au lieu de géné-
rer de façon conventionnelle y j+1 selon q(y j+1 ), cette méthode permet de générer en

premier une direction u prop
j+1 selon q̃ u prop
j+1 et ensuite une distance radiale s prop
j+1 selon

q̃ s prop prop
j+1 u j+1 . À ce point, il reste seulement à poser la forme des deux densités ins-
trumentales. Il est bon de rappeler que l’objectif est de choisir ces densités de façon à
ce qu’elles approchent le mieux possible les densités cibles correspondantes. Générale-
ment, une bonne approximation de la distribution directionnelle cible est une distribu-
tion uniforme sur une sphère de rayon 1 dans Rn . En d’autres mots, on considère que
π̃ (u) ≈ q̃ (u) = 1/An , où An = 2π n/2 /Γ (n/2) représente l’aire de la sphère. Ensuite,
dans un contexte général, la densité q̃ (s |u) est choisie comme la meilleure approxima-
tion disponible de la densité conditionnelle π̃ (s |u).

Dans le contexte présent, la densité instrumentale globale q y j+1 est voulue de la
forme Student f . Par conséquent, en utilisant l’approximation par une distribution uni-
forme sur la sphère de rayon 1, une direction proposée est générée par :
(z1 , . . . , zn )
(u prop 0
j+1 ) = , où z = (z1 , . . . , zn ) et zi ∼ N(0, 1). (2.5)
|z|
Ensuite, une distance radiale proposée est générée par s prop

j+1 = ( f + n)
1/2 |z|/χ où
f
|z| est le même vecteur ayant servi à obtenir la direction. En utilisant cette méthode, il
27

est clair, de (2.3), que la distribution instrumentale globale q y j+1 est une loi Student f
avec densité donnée par (2.4).
Enfin, l’ajustement directionnel est effectué par la sélection du degré de liberté f . Le
choix optimal est différent d’une direction à l’autre afin d’approximer le mieux possible

la densité conditionnelle π̃ s prop
j+1 u prop
j+1 . En fait, f est choisi tel que les décroissances
des densités cible et instrumentale à un point s∗ relativement au mode soient les mêmes.
En d’autres mots, le ratio de la densité cible à une distance donnée s∗ et au mode est égal
à celui de la densité instrumentale

π ∗ u prop
j+1 · s ∗ π x̂ + Ĥ −1/2 u prop
j+1 · s ∗ q f prop u prop
j+1 · s ∗
j+1
= =
π ∗ (0) π (x̂) q f prop (0)
j+1
prop
f j+1 +n
 0 − 2
prop ∗
u j+1 · s u prop ∗
j+1 · s 
= 1 + (2.6)

prop
f j+1 +n

où q f prop est donnée par (2.4).

j+1
prop
Une solution analytique de f j+1 est ardue, mais en pratique il est suffisant de choi-
sir la valeur la plus appropriée parmi les entiers entre 1 (distribution de Cauchy) et 50
(distribution presque normale).
n o
En posant r2 = 2 log π x̂ + Ĥ −1/2 u propj+1 · s∗ /π (x̂) , Q2 =
0
u prop
j+1 · s∗ u prop ∗
j+1 · s prop
et f j+1 + n = f , il suffit de trouver la meilleure solu-
prop
tion f j+1 ∈ {1, 2, . . . , 50} de l’équation
Q2

f log 1 + = r2 . (2.7)
f
Donc, la méthode DA utilise une densité instrumentale globale qui peut être vue
comme un amalgame de densités Student f , chacune ayant comme support une direction
donnée à partir du mode et un degré de liberté f spécifié par cette direction. L’approche
28
générale peut être illustrée de façon simpliste pour le cas unidimensionnel par la figure
2.1, reproduite de Bédard et Fraser (2008).
Le seul paramètre libre restant est s∗ , c’est-à-dire la distance où le ratio (2.6) des
densités cible et instrumentale est évalué. En un premier lieu, on considère le cas où π
est une densité normale multivariée avec composantes indépendantes distribuées selon
une loi normale standard. Dans ce cas, les degrés de liberté proposés par la méthode
DA seront fixes à 50 dans chaque direction et la distribution instrumentale sera approxi-
mativement π. Donc, dans ce contexte, u prop
j+1 peut être généré selon (2.5), et la distance
radiale proposée peut être approximée par (s prop prop
j+1 ) = |z|. La distribution de (s j+1 ) sera
√
donc χ avec n degrés de liberté et une espérance approximative de n. Pour cette rai-
son, les auteurs de la méthode DA suggèrent d’exprimer la valeur s∗ en termes de cette
√
espérance et d’un paramètre de rodage λ , c’est-à-dire s∗ = λ n. Dans Bédard et Fraser
(2008), il est aussi suggéré d’utiliser une valeur de λ autour de 2 ou 3.
En fait, le choix de ce paramètre peut avoir un impact considérable sur la perfor-
mance de la méthode DA et une grande partie de ce chapitre sera consacrée à ce sujet.
Cependant, une brève description de l’algorithme est tout d’abord nécessaire afin de
clarifier et de résumer les étapes de cette approche.
2.2 Étapes de l’algorithme DA
En un premier lieu, il faut déterminer les valeurs x̂, Ĥ 1/2 et Ĥ −1/2 . Dans le logiciel R,
ces quantités peuvent être évaluées à l’aide des fonctions nlm, fdHess, pdMat et pdFactor
(greffon nlme) .
Ensuite, l’algorithme procède comme suit :
1 Déterminer la direction de la valeur initiale reparamétrisée x∗0 = Ĥ 1/2 (x0 − x̂),

x∗0
c’est-à-dire |x∗0 | et trouver f0 ∈ {1, 2, · · · , 50}, la meilleure solution de (2.7).
2 Pour j ≥ 0
29
1.0
0.8
0.6
0.4
0.2
●
0.0
−2 0 2 4
Figure 2.1 – Graphique de la densité cible (ligne pleine) et instrumentale (ligne poin-
tillée) lorsque la direction choisie est vers la droite (+1). Le point d’intersection est
u prop ∗ ∗
j+1 · s = s .
a) Générer u prop
j+1 à partir de la relation (2.5) et stocker la valeur |z j+1 |.
prop
b) Trouver f j+1 ∈ {1, 2, · · · , 50}, la meilleure solution de (2.7).
c) Obtenir une distance radiale s prop prop

j+1 = ( f j+1 + n)
1/2 |z
j+1 |/χ f prop , où χ f prop est
j+1 j+1
une variable indépendante distribuée selon une loi Chi et |z j+1 | est la valeur
trouvée en (a).
d) Poser y∗j+1 = u prop prop

j+1 · s j+1 .
e) Calculer la probabilité d’acceptation

x̂ + −1/2 y∗ ∗
∗ ∗
π Ĥ j+1 q f j x j
α x j , y j+1 = 1 ∧ , (2.8)
π x̂ + Ĥ −1/2 x∗j q f prop y∗j+1
j+1
où q f (x) est donnée par (2.4).
f) Générer une valeur r j+1 d’une loi U[0, 1].

30

g) Si r j+1 ≤ α x∗j , y∗j+1 , accepter la valeur proposée et poser x∗j+1 = y∗j+1 ,
prop
f j+1 = f j+1 . Sinon, rejeter la valeur proposée et poser x∗j+1 = x∗j , f j+1 = f j .
h) Obtenir x j+1 = x̂ + Ĥ −1/2 x∗j+1 .
i) Poser j = j + 1 et revenir à l’étape (a), pour N itérations.
2.3 Impact de la valeur de λ sur la performance de l’algorithme DA
Il a été brièvement mentionné que le paramètre λ pouvait avoir un impact sur la

performance de la méthode DA. Le reste de ce chapitre se penchera sur ce sujet. En
premier, deux exemples tirés de l’article Bédard et Fraser (2008) seront analysés en
fonction du paramètre λ et enfin un exemple unidimensionnel sera construit afin de
démontrer les problèmes de convergence potentiels de l’algorithme.
2.3.1 Exemple 1
Dans un premier problème de l’article Bédard et Fraser (2008), on cherche à obtenir

un échantillon d’une distribution de paramètres de régression. À l’origine, cet exemple
a été étudié dans Bédard et al. (2007) et se base sur un modèle de régression yi = α +
β xi + σ zi , où y est la variable réponse, x la variable explicative et z représente l’erreur
distribuée selon une loi Student7 . Les données suivantes ont été en fait générées avec
α = 0, β = 1, σ = 1 :
xi -3 -2 -1 0 1 2 3
.
yi -2,68 -4,02 -2,91 0,22 0,38 -0,28 0,03
On suppose qu’il y a un intérêt particulier pour l’hypothèse β = 1 et que, d’un point

de vue bayésien, la distribution à postériori des trois paramètres (α, β , τ = log σ ) est
donnée par
7 −4
(y0i − α − β xi )2

0 −7τ

π1 α, β , τ y dα dβ dτ = c e ∏ 1+ 2τ
dα dβ dτ. (2.9)
i=1 7e
31
Dans ce contexte, π1 sera la densité cible et l’objectif est une estimation de la valeur-
s pour l’hypothèse β = 1, c’est-à-dire s(β = 1) = P(β ≥ 1). Étant donné un échantillon
de la distribution cible, une approximation de cette quantité est donnée par la valeur-s
empirique :
N N
1 1
s(β ) =
N ∑ I βj ≥ β = N ∑I βj ≥ 1 ,
j=1 j=1
où N est la taille de l’échantillon généré. Pour cet exemple, il a été démontré dans Bédard
et Fraser (2008) que la méthode DA offre une meilleure performance comparativement à
d’autres méthodes RWMH et IS sélectionnées par les auteurs, résultant en une variance
empirique de s(β ) inférieure.
√
Afin de mesurer l’impact du paramètre s∗ = λ n sur la performance de l’algorithme
DA, trois mesures empiriques ont été calculées pour différentes valeurs de λ . En premier,
l’écart-type de la valeur-s pour β = 1 a été obtenue. Dans ce contexte, une légère variante
de la mesure de variance par moyenne de séries présentée au chapitre 1 a été utilisée.
Une période de chauffe initiale de 10 000 itérations a été écartée et ensuite, 4 000 000
valeurs générées ont été divisées en 4000 séries de 1000 observations. Dans chacune
de ces séries, les cinquante premières valeurs ont été rejetées afin de minimiser encore
davatage la corrélation entre les séries. Pour des raisons de cohérence, ces paramètres
seront fixes sauf avis contraire pour tous les algorithmes simulés dans le cadre de ce
mémoire. Donc, en gardant la même notation qu’au chapitre 1, a = 4000, b = 950 et la
variance par moyenne de séries est donnée par :
950 4000
σ̂λ2 = ∑ (sl − sN )2
3999 l=1
1 4000
où sl est la valeur-s empirique de la l-ième série et sN = 4000 ∑l=1 sl est la valeur-s
moyenne globale. La figure 2.2 présente l’écart-type empirique de s(β ) et démontre
qu’une valeur de λ ∈ [2; 4] semble minimiser cette mesure.
32
0,03
0,025
0,02
σ 0,015
0,01
0,005
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.2 – Graphique de l’écart-type σ de s(β = 1) en fonction de λ (algorithme DA,

exemple 1).
Ensuite, le taux d’acceptation de l’algorithme a été utilisé à titre de deuxième mesure

empirique. La figure 2.3 indique que le taux d’acceptation semble être maximal autour
des mêmes valeurs λ .
Un grand taux d’acceptation est généralement insuffisant pour garantir une bonne
performance d’un algorithme MCMC. Cependant, dans le cas présent, si l’algorithme
explore bien l’espace et ne sous-estime pas les queues tel qu’argumenté dans Bédard et
Fraser (2008), on peut présumer qu’un meilleur taux d’acceptation est préférable.
Selon la figure 2.4, qui présente la distance de saut carrée moyenne, une valeur de λ
de 3 ou 4 semble encore être idéale. Comme mentionné à l’équation (1.6), la DSCM est
reliée à l’autocorrélation de la chaîne. On peut donc aussi s’attendre à une autocorréla-
tion totale d’intervalle 1 plus faible lorsque l’algorithme emploie une valeur de λ entre
2 et 4.
D’un autre côté, la densité instrumentale de l’algorithme DA est ajustée selon la
33
0,8
0,7
0,6
0,5
Taux
0,4
d'acceptation
0,3
0,2
0,1
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.3 – Graphique du taux d’acceptation en fonction de λ (algorithme DA, exemple

1).
direction choisie et une mesure potentiellement intéressante est la moyenne des degrés
de liberté de la densité instrumentale en fonction de la valeur de λ . Pour cet exemple,
il a été montré dans Bédard et Fraser (2008) que la moyenne des degrés de liberté de
la densité instrumentale pour une valeur de λ = 2, se situe à 28,57. Cependant, lorsque
l’on suit son évolution, on peut discerner une tendance de cette moyenne à augmenter en
fonction du paramètre λ . La figure 2.5 présente ces résultats. Ce graphique indique que la
vitesse moyenne de décroissance des queues par rapport au mode n’est pas constante. En
d’autres mots, la densité cible décroît de façon relativement lente au début, s’apparentant
à une densité Student27 , mais au fur et à mesure que la distance par rapport au mode
augmente, la vitesse de décroissance des queues augmente et à λ = 20, la densité cible
s’apparente plutôt à une densité Student33 . Cette idée est importante et sera explorée tout
au long de ce chapitre.
Tel que mentionné auparavant, l’algorithme DA a été comparé et jugé supérieur à des
34
0,8
0,7
0,6
DSCM 0,5
0,4
0,3
0,2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.4 – Graphique de la DSCM en fonction de λ (algorithme DA, exemple 1).
algorithmes RWMH et IS spécifiques implémentés dans Bédard et Fraser (2008). Dans

un même ordre d’idées, il peut être intéressant de comparer la méthode DA avec des algo-
rithmes IS se situant aux extrémités de la famille des distributions Student, par exemple
q1 (x) = Student1 (x̂, ( f + n)Ĥ −1 )(Cauchy) et q50 (x) = Student50 (x̂, ( f + n)Ĥ −1 ) (qui
sera, pour des raisons de simplicité, appelée normale).
La supériorité de l’approche DA devient évidente lorsque les trois algorithmes sont
comparés selon les diagnostics de convergence présentés précédemment. La figure 2.6
démontre que la distance de saut carrée moyenne est plus élevée pour n’importe quelle
valeur de λ pour la méthode DA.
Ensuite, la figure 2.7 indique que le taux d’acceptation est plus élevé pour la mé-
thode DA que pour l’échantillonneur indépendant avec distribution instrumentale Cau-
chy. Étant donné que ces deux algorithmes ne sous-estiment pas les queues de la densité
cible, nous déduisons que la méthode avec le taux d’acceptation le plus élevé devrait être
préférée. Notons qu’ici, le taux d’acceptation doit être interprété avec prudence dans le
35
35
33
31
Degrés de liberté
29
moyens
27
25
23
0 5 10 15 20 25
Paramètre λ
Figure 2.5 – Graphique des degrés de liberté moyens (avec un écart-type) en fonction de
λ (algorithme DA, exemple 1).
cas de l’échantillonneur indépendant avec distribution instrumentale normale. En effet,

cet algorithme aura tendance à sous-estimer les queues de l’exemple considéré et le taux
d’acceptation s’en trouvera (artificiellement) gonflé.
Finalement, la figure 2.8 révèle que l’écart-type de s(β = 1) obtenu par la méthode
DA est inférieur à celui obtenu par l’algorithme IS avec distribution Cauchy pour une
valeur λ ≤ 10 et supérieur pour λ > 10. Cependant, l’algorithme IS avec densité instru-
mentale normale résulte en un écart-type 6 fois plus élevé que les deux autres approches.
Ce dernier résultat indique un problème potentiel de convergence dans le cas de l’algo-
rithme IS avec densité instrumentale normale. En fait, pour cet exemple, cette méthode
converge beaucoup plus lentement que les deux autres approches pour des raisons qui
seront abordées plus loin dans ce chapitre.
Donc, ces résultats réitèrent que l’approche DA, malgré un effort computationnel
accru, peut mener à des résultats supérieurs à d’autres approches IS lorsque peu d’in-
36
0,8
0,7
0,6
DSCM 0,5 DA
0,4 Normale
Cauchy
0,3
0,2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.6 – Graphique de la DSCM des algorithmes IS et DA (exemple 1).
formations au sujet de la densité cible sont disponibles. De plus, il semble qu’un choix
optimal existe pour le paramètre λ et puisse avoir un impact sur la performance de la
méthode. L’écart-type de la valeur-s a aussi été calculé pour d’autres valeurs de β avec
des résultats identiques.
2.3.2 Exemple 2
Le second exemple analysé est aussi tiré de l’article Bédard et Fraser (2008), à l’ori-
gine étudié dans Cox et Snell (1981). Dans ce contexte, on cherche encore à obtenir des
mesures relativement à des paramètres de régression. Les données concernent 32 réac-
teurs à eau légère situés aux États-Unis. La variable réponse est le coût d’investissement
d’un nouveau réacteur et il y 10 variables explicatives dont 7 sont considérées signifi-
catives. Le modèle de régression est posé par y = Xβ + σ z, où X est une matrice 32x7
contenant les variables explicatives et z ∼ Student4 représente l’erreur. La distribution
37
0,8
0,7
0,6
0,5
Taux
0,4 DA
d'acceptation
0,3 Normale
0,2 Cauchy
0,1
0
1 3 5 7 9 11 13 15 17 19
Paramètre λ
Figure 2.7 – Graphique du taux d’acceptation des algorithmes IS et DA (exemple 1).
de la variable réponse est

7
−7 yi − α − xi β
f (y; α, β , σ ) dy = σ ∏h dyi ,
i=1 σ
où h(·) est la densité Student4
D’un point de vue classique, après une reparamétrisation dont les détails seront omis
(voir Bédard et Fraser (2008)), la distribution cible est donnée par
n
π2 b, a d0 db da = c ∏ h ea di0 + Xi b ea(n−r) db da .

(2.10)
i=1
La densité π2 est en huit dimensions, b représente le vecteur des sept coefficients

de régression et s = exp(a) représente l’écart-type de l’erreur de régression. Le terme di0
représente le i-ème résidu observé standardisé ; ces résidus satisfont d0 = y0 − Xb0 /s0

où b0 dénote les estimateurs par moindres carrés et (s0 )2 = ∑ni=1 (yi − ŷi )2 /(n − r), avec
n = 32 et r = 7. Le terme Xi dénote la i-ème rangée de la matrice des variables expli-
catives. Le lecteur est invité à consulter Bédard et Fraser (2008) pour de plus amples
38
0,07
0,06
0,05
0,04
σ 0,03
DA
Normale
0,02
Cauchy
0,01
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.8 – Graphique de l’écart-type σ de s(β = 1) des algorithmes IS et DA (exemple

1).
détails ainsi que le tableau complet des valeurs observées (yi ) et explicatives (xi ). Pour
l’objectif de ce chapitre, la forme générale (2.10) de la distribution cible est suffisante.
Dans ce contexte, il y a un intérêt particulier pour le sixième coefficient de régression,
c’est-à-dire la dépendance du coût d’investissement par rapport à la sixième variable
explicative qui dans ce cas correspond au nombre de réacteurs construits par un ingénieur
donné. L’objectif est d’obtenir une estimation de la valeur-p pour l’hypothèse nulle β6 =
−0, 01. De manière semblable à l’exemple précédent, un échantillon sera généré de la
distribution cible et une valeur-p empirique sera utilisée comme estimateur.
Pour cet exemple, il a été démontré dans Bédard et Fraser (2008) que la performance
de l’algorithme avec ajustement directionnel est supérieure à celle des algorithmes IS et
RWMH particuliers considérés dans cet article. En fait, l’approche RWMH peut présen-
ter de graves problèmes de convergence pour cet exemple et les détails en seront abordés
au chapitre 3.
39
Afin de déterminer l’effet du paramètre λ sur la performance de l’algorithme DA,

les mêmes diagnostics de convergence empiriques qu’à la section 2.3.1 ont été utilisés.
Les figures 2.9 et 2.10 démontrent que la distance de saut carrée moyenne et le taux
d’acceptation sont maximisés lorsque λ est petit, autour de 2 ou 3. La figure 2.11 montre
que l’écart-type de la valeur-p calculée dépend de façon moins marquée de λ , mais
atteint un minimum pour λ = 2 ou 3. Les relations de l’écart-type de la valeur-p en
fonction de λ pour d’autres valeurs de β6 ont été omises puisque similaires.
De manière semblable à l’exemple 2.3.1, la figure 2.12 de la moyenne des degrés de
liberté des densités instrumentales donne un aperçu du comportement des queues de la
densité cible. Dans ce cas, les queues de la densité cible sont nettement moins épaisses
qu’à l’exemple 2.3.1. En plus, il semble que les queues ont un taux de décroissance
rapide jusqu’à λ = 13 et ce taux ralentit à partir de λ > 13.
Dans cette situation, il existe un choix optimal pour le paramètre λ , tout comme pour
l’exemple 2.3.1, et il se situe autour des mêmes valeurs λ ∈ [2; 4].
2.4 Convergence de l’algorithme DA
Dans la section 2.3, trois diagnostics de convergence ont été utilisés pour juger de
l’importance du paramètre λ dans deux exemples réalistes. À ce point, il semble que
le choix de ce paramètre peut légèrement influencer la performance de l’algorithme et
qu’un choix sécuritaire semble être situé entre 2 et 4. Dans cette section, un exemple
sera construit pour démontrer que le choix de ce paramètre peut parfois être crucial et se
trouver très loin de ces valeurs « sécuritaires ».
Cependant, un résultat théorique est d’abord nécessaire. Au chapitre 1, des notions
d’uniformité ergodique et géométrique ont été présentées de façon générale. Bien que
leur application soit souvent difficile en pratique, dans le cas d’un algorithme IS et par
conséquent dans le cas de la méthode DA, il existe certaines conditions simples qui
garantissent ces propriétés.
40
1350
1300
1250
1200
DSCM
1150
1100
1050
1000
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.9 – Graphique de la DSCM en fonction de λ (algorithme DA, exemple 2).
Théorème 8. L’échantillonneur indépendant est uniformément ergodique s’il existe une

constante β > 0 telle que
q(x)
≥ β, ∀x ∈ X , (2.11)
π(x)
et alors ||PN (x, ·) − π(·)|| ≤ (1 − β )N , où N est le nombre d’itérations.
D’un autre côté, si ess inf{q(x)/π(x)} = 0 en π-mesure, l’algorithme n’est même pas
géométriquement ergodique. La borne inférieure essentielle est définie comme le supre-
mum des presque minorants d’une fonction, et donc, dans ce cas, ess infπ {q(x)/π(x)} =
q(x)
sup{a ∈ R|π({x| π(x) < a}) = 0}
Démonstration. Voir Mengersen et Tweedie (1996).
Donc, afin de garantir une convergence uniformément ergodique, il suffit de trouver

une constante 0 < β ≤ 1 (la dernière borne résulte du fait que les deux fonctions sont
des densités normalisées et intègrent à 1) telle que le ratio des densités instrumentale
41
0,74
0,72
0,7
0,68
Taux
0,66
d'acceptation
0,64
0,62
0,6
0,58
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.10 – Graphique du taux d’acceptation en fonction de λ (algorithme DA,

exemple 2).
et cible en soit supérieur pour tous les éléments de l’espace. En général, X = Rn avec
mesure sous-jacente de Lebesgue et il n’existe pas d’ensemble borné A tel que q(A) = 0
et π(A) > 0, puisque cela violerait la condition de φ -irréductibilté. En plus, la démons-
tration du théorème 8 dans Mengersen et Tweedie (1996) suppose q(x) > 0, π(x) > 0
pour tout x ∈ X. Ainsi, afin de trouver la borne inférieure essentielle du ratio, il est sou-
vent nécessaire d’en comprendre le comportement limite. Ce concept sera exploré tout
au long de cette section.
En un premier lieu, en revenant à la section 2.3.1, on se rappelle que l’échantillon-
neur indépendant (IS) avec distribution instrumentale normale semblait inefficace selon
les diagnostics de convergence présentés aux figures 2.6 et 2.8. Selon le graphique 2.5
de la moyenne des degrés de liberté proposés par la méthode DA pour cet exemple, il
semble que les queues de la densité cible s’apparentent en moyenne à celles d’une den-
sité Student avec degrés de liberté entre 27 et 33. Si on considère momentanément pour
42
0,012
0,01
0,008
σ 0,006
0,004
0,002
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.11 – Graphique de l’écart-type σ de la valeur-p pour H0 : β6 = −0, 01 (algo-

rithme DA, exemple 2).
l’exemple 1 que π1 (x) ≈ Student27 , alors
Γ( 50+3
− 50+3
2 )

x0 x 2
1 + 50+3 (50 + 3)−3/2
q50 (x) π 3/2 Γ( 50
2 )
≈ − 27+3 .
π1 (x) Γ( 27+3
2 )

x0 x 2
−3/2
π Γ( 27
1 + 27+3 (27 + 3)
2 )
3/2
Il est facile de voir que la borne inférieure essentielle de ce ratio est 0. Par exemple,
sur une ligne x = 0 + t1, la limite lorsque t → ∞ est donnée par
2 ( −50−3

) 2 ( 27+3 )
lim (50 + 3 + 3t ) 2 (27 + 3 + 3t ) 2 = 0.
t→∞
Ainsi, il semble probable que la raison expliquant la piètre convergence de l’algo-

rithme IS avec densité instrumentale normale (degrés de liberté 50) soit la décroissance
rapide de ses queues. Donc, des algorithmes présentés, seuls les méthodes DA et IS avec
distribution instrumentale Cauchy pourraient potentiellement converger de façon uni-
43
52
50
48
46
Degrés de liberté
44
moyens
42
40
38
36
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Paramètre λ
Figure 2.12 – Graphique des degrés de liberté moyens en fonction de λ (algorithme

DA, exemple 2). L’écart-type type n’est pas montré puisque plus petit que la largeur des
symboles.
formément ergodique et cela semble être corroboré par les diagnostics de convergence
présentés à la section 2.3.1.
À ce point, il est temps d’examiner un inconvénient potentiel de la méthode DA. Tel
que mentionné à la section 2.1, l’algorithme approxime la densité cible en se basant sur
√
le ratio (2.6) évalué à une distance s∗ = λ n :

−1/2 prop ∗ prop ∗
π x̂ + Ĥ u j+1 · s qf prop u j+1 · s
j+1
=
π (x̂) q f prop (0)
j+1
prop
f j+1 +n
 0 − 2
prop ∗ prop ∗
u j+1 · s u j+1 · s 
= 1 + .

prop
f j+1 +n

Dans certains cas, il peut arriver que la densité cible décroisse rapidement sur une
région, s’apparentant à une densité normale, mais que les queues restent épaisses au-delà
de cette région. Dans ce contexte, la méthode DA proposera une densité sous-optimale
44
(quasi-normale) si s∗ se trouve dans la région de forte décroissance tandis qu’elle pro-

posera une densité plus optimale (plus épaisse) si s∗ se trouve dans la deuxième région.
Cette idée est illustrée à l’aide de l’exemple unidimensionnel suivant.
Soit la fonction de densité suivante

 3 (2 − x2 ), |x| < 1
16
f (x) = . (2.12)
 3 1, |x| ≥ 1
16 x 2
Pour cet exemple, la densité est représentée à la figure 2.13, illustrant la symétrie et un
mode unique à 0. De plus, cette densité est lisse dans le sens qu’aux points {−1, 1}, la
fonction est continue et la valeur des dérivées premières coïncident. La matrice hessienne
du négatif du log de la fonction au mode est donnée par
d2 d2

3 2
(− log( f (x))) = − log − log(2 − x )
dx2 x=0 dx2 16 x=0

d 2x
=
dx 2 − x2 x=0
4 + 2x2
= = 1.
(2 − x2 )2 x=0
Il n’est donc pas nécessaire d’ajuster la densité cible dans cet exemple puisque le
mode et la matrice hessienne correspondent déjà à ceux de la densité instrumentale, qui
dans ce cas sera de la forme suivante

Γ f +1
2

x2
− f +1
2
q f (x) = 1+ ( f + 1)−1/2 .
1/2
π Γ 2 f f + 1
Dans le cas unidimensionnel, la direction proposée par l’algorithme DA est choisie

√
au hasard entre -1 et 1. Le paramètre s∗ = λ 1 = λ peut être choisi près du mode, par
exemple λ = 1. D’autres valeurs de λ petites produiront des résultats semblables, mais
λ = 1 ici simplifie l’illustration analytiquement.
Ainsi, les quantités nécessaires de (2.7) pour l’évaluation des degrés de liberté pro-
posés seront simplement
n o
2 prop ∗
r = 2 log π (0) /π u j+1 · s = 2 log(2) (2.13)

Q2 = u prop
j+1 · s∗
u prop ∗
j+1 · s = 1. (2.14)
45
0,4
0,35
0,3
0,25
0,2
0,15
0,1
0,05
0
-6 -4 -2 0 2 4 6
Figure 2.13 – Graphique de la fonction de densité (2.12).
L’équation (2.7) ne sera jamais vraie pour aucune valeur de f puisque la partie
gauche est tout au plus 1 (selon une des définitions de la constante e = limn→∞ (1 +
1/n)n ). Cependant, la partie gauche croît avec f et par conséquent l’algorithme propo-
sera la plus grande valeur possible f = 50.
Toutefois, les queues de la densité cible sont dérivées de densités épaisses de Pareto
et donc il semble qu’une densité instrumentale presque normale ne soit pas adéquate.
Selon le théorème 8, il serait nécessaire de trouver la borne inférieure essentielle du ratio
afin de juger de la convergence de l’algorithme.
Le ratio des deux distributions est donné par :

Γ f +1 − f +1
2 x2 2

f
1+ f +1 ( f + 1)−1/2
q(x) π 1/2 Γ 2
= 3 1
π(x) 16 x2
f +1
f +1
16Γ ( f + 1) 2
2 x2
= f +1 ,
3π 1/2 Γ 2f ( f + 1)1/2 ( f + 1 + x2 ) 2
46
pour |x| ≥ 1. En prenant la limite, il en résulte que :

!
x2

q(x) q(x)
lim = lim = lim C f +1
x→−∞ π(x) x→∞ π(x) x→∞
( f + 1 + x2 ) 2
 
2x
= C lim  f −1

x→∞ f +1 2
2 ( f + 1 + x ) 2x
2
 
1
= C lim  f −1

x→∞ f +1 2
2 (f +1+x ) 2
et finalement
 
32
√
C si f = 1 si f = 1

q(x) 
3π 2
lim = = .
n→∞ π(x)  0 si f > 1  0 si f > 1
En examinant le comportement limite, il est clair qu’un algorithme DA proposant

une densité instrumentale normale ne converge pas de façon uniformément ni géométri-
q(x)
quement ergodique car il n’existe pas de β > 0 tel que π(x) ≥ β pour tout x ∈ R. La seule
densité instrumentale garantissant une ergodicité uniforme est la densité Cauchy. Dans
ce cas,
q(0) 8
β= = √ ≈ 0, 6
π(0) 3π 2
et donc la vitesse de convergence devrait être, selon le théorème 8, de l’ordre de (1 −
0, 6)N = 0, 4N , où N est le nombre d’itérations.
Il peut être intéressant de simuler cet exemple et d’utiliser certaines des mesures de
convergence établies précédemment afin d’illustrer l’idée d’un point de vue pratique. La
figure 2.14 présente les degrés de liberté moyens proposés, la valeur s(x = 1) (dont la
3
valeur exacte est 16 ) obtenue et son écart-type ainsi que l’autocorrélation en fonction de
λ . De ces données, il est clair que la densité instrumentale optimale est utilisée lorsque
λ = 42, ce qui est loin des valeurs [2; 4] établies à la section précédente. Donc, un certain
degré de prudence est nécessaire lorsque la méthode DA est employée. Il est à mention-
ner, au cas où l’absence d’espérance mettrait en doute la signification pratique du dernier

47

ϲϬ Ϭ͕Ϯϴ
Ϭ͕Ϯϲ
ϱϬ
Ϭ͕Ϯϰ
ϰϬ Ϭ͕ϮϮ
ĞŐƌĠƐĚĞ
sĂůĞƵƌͲƐ Ϭ͕Ϯ
ůŝďĞƌƚĠ ϯϬ
ƉŽƵƌǆсϭ Ϭ͕ϭϴ
ƉƌŽƉŽƐĠƐ
ϮϬ Ϭ͕ϭϲ
Ϭ͕ϭϰ
ϭϬ
Ϭ͕ϭϮ
Ϭ Ϭ͕ϭ
Ϭ ϮϬ ϰϬ ϲϬ Ϭ ϮϬ ϰϬ ϲϬ
ʄ ʄ

Ȝ =1, f=50 Ȝ=41, f=2 Ȝ=42, f=1
/ŶƚĞƌǀĂůůĞ /ŶƚĞƌǀĂůůĞ /ŶƚĞƌǀĂůůĞ

Figure 2.14 – Graphique des degrés de liberté proposés, de la valeur s(x = 1) obtenue

(avec un écart-type) et de l’autocorrélation en fonction de λ (algorithme DA, exemple
2.12).
exemple, que les résultats sont sensiblement les mêmes pour une distribution possédant
une espérance ainsi qu’une variance et construite de manière identique :

3 2
16 (3 − 2x ), |x| < 1

f (x) = .
3 1

16 x4 , |x| ≥ 1
À l’inverse, il est anticipé qu’il existe également des exemples où l’algorithme DA

serait amené à proposer une densité instrumentale avec des queues épaisses alors que le
48
choix optimal se trouverait être une densité aux queues minces.

L’utilité de l’exemple précédent était d’illustrer le fait que l’algorithme DA, étant
basé seulement sur le ratio des densités par rapport au mode, peut proposer une densité
instrumentale sous-optimale. Afin de réduire le risque de telles situations, il est toujours
recommandé de simuler l’algorithme avec plusieurs valeurs de λ . Ainsi, il est possible
d’étudier le comportement des queues de la densité cible ainsi que l’optimalité de cer-
taines mesures de convergence.
Une autre manière d’éviter une densité instrumentale inefficace est de créer un al-
gorithme avec un mélange entre une densité instrumentale Cauchy et celle proposée par
l’algorithme DA. La densité Cauchy possède des queues épaisses, plus épaisses qu’une
vaste majorité de densités. En utilisant un tel mélange, il est plus probable d’obtenir
une convergence uniformément ergodique de l’algorithme. Par exemple, si on désigne la
densité instrumentale globale induite par la méthode DA par q f (x) et la densité Cauchy
par q1 (x) alors un mélange q(x) peut être construit selon
q(x) = pq1 (x) + (1 − p)q f (x),
où 0 < p < 1.
Si l’algorithme avec densité instrumentale Cauchy est uniformément ergodique, il en
résulte qu’il existe un β > 0 tel que q1 (x)/π(x) ≥ β pour x ∈ Rn et alors
q(x) pq1 (x) + (1 − p)q f (x)

=
π(x) π(x)
q f (x)
≥ pβ + (1 − p)
π(x)
≥ pβ
et la vitesse de convergence de cet algorithme sera ||PN (x, ·) − π(·)|| ≤ (1 − pβ )N .

Il est clair que cette méthode est moins efficace si la méthode DA est déjà uniformé-
ment ergodique, mais selon notre expérience l’impact n’est pas énorme pour des valeurs
de p petites, typiquement p ≤ 0, 05. En fait, si la méthode DA est déjà uniformément
49
ergodique, il en resulte qu’il existe β 0 > 0 tel que q f (x)/π(x) ≥ β 0 et que la vitesse de
convergence de cet algorithme est (1 − β 0 + p(β 0 − β ))N .
En contrepartie, si l’approche DA n’est pas uniformément ergodique, l’addition de la
densité Cauchy (ou possiblement d’autres densités aux queues épaisses, par exemple les
densités de Pareto) pourrait rendre l’algorithme uniformément ergodique. Dans ce cas, la
vitesse de convergence sera plus faible que celle d’un algorithme IS avec distribution ins-
trumentale Cauchy, par un facteur de p. Cette dernière méthode est simple, rapidement
implémentable et résulte en une approche un peu plus versatile se situant entre celle du
DA et du IS. Il est clair qu’elle est utile seulement lorsque l’on est incertain de la conver-
gence de la méthode DA et elle représente en quelque sorte un compromis entre une plus
grande probabilité d’uniformité ergodique versus une pire vitesse de convergence.
En conclusion, la méthode DA est une approche intéressante qui vise à combiner la
versatilité de l’algorithme RWMH avec la haute performance de l’algorithme IS. Elle
se base sur un paramètre λ qui doit être choisi de façon éclairée. En plus, l’approche
DA peut dans certains contextes proposer une densité très sous-optimale ce qui peut être
évité par une meilleure étude de la densité cible ou l’introduction d’un mélange avec une
densité aux queues plus épaisses. En fait, il n’existe pas de méthode infaillible et même
un algorithme des plus versatiles comme celui du RWMH peut présenter des problèmes
de convergence importants comme dévoilé dans le prochain chapitre.
CHAPITRE 3
COMPARAISON ENTRE L’ALGORITHME DA ET DES MÉTHODES

LOCALES
Au dernier chapitre, il a été question de l’algorithme avec ajustement directionnel

(DA) qui visait à établir une méthode combinant la versatilité de l’algorithme Metro-
polis de type marche aléatoire (RWMH) et l’efficacité de l’échantillonneur indépendant
(IS). Les méthodes globales, c’est-à-dire les algorithmes où la distribution instrumentale
ne dépend pas de l’état présent (comme les méthodes IS et DA), nécessitent une connais-
sance plus détaillée de la densité cible. Les méthodes locales, c’est-à-dire les algorithmes
où la distribution instrumentale est dépendante de l’état présent (comme l’algorithme
RWMH), sont plus versatiles et facilement applicables dans plusieurs contextes.
Dans Bédard et Fraser (2008), l’algorithme DA a été comparé et jugé supérieur aux
algorithmes IS avec distribution instrumentale Student7 et RWMH avec distribution ins-
trumentale normale aux composantes indépendantes pour les deux exemples spécifiques
des sections 2.3.1 et 2.3.2. Cependant, comme mentionné au premier chapitre, il existe
plusieurs méthodes MCMC avancées, la plupart dérivées des algorithmes élémentaires
présentés à la section 1.2. Par conséquent, il peut être intéressant de comparer l’approche
DA avec d’autres algorithmes plus sophistiqués. L’analyse présentée dans ce chapitre
s’attardera uniquement sur des méthodes alternatives locales puisque leur versatilité est
comparable et/ou supérieure à l’algorithme DA. En effet, il est fort probable qu’il existe
des méthodes alternatives globales dont la performance est supériure à celle de la mé-
thode DA, mais elles nécessiteront considérablement plus d’informations au sujet de la
densité cible. Donc, l’idée est de déterminer si des méthodes plus facilement applicables
que celle du DA peuvent mener à des résultats semblables. Ce chapitre vise à explo-
rer cette comparaison et à cette fin, l’exemple 2 (section 2.3.2) du chapitre dernier sera
utilisé comme paradigme.
51
3.1 Algorithmes RWMH
À titre de rappel, la distribution cible est une distribution de paramètres de régression

et est donnée par
n
π2 b, a d0 db da = c ∏ h ea di0 + Xi b ea(n−r) db da.

i=1
La coefficient β6 , correspondant au nombre de réacteurs construits par un ingénieur, est

toujours la variable d’intérêt. Pour échantillonner de cette distribution, un choix de mé-
thode locale des plus simples est l’algorithme RWMH avec distribution instrumentale
normale aux composantes indépendantes. Cependant, appliqué à cet exemple, cet algo-
rithme présente des problèmes de convergence assez importants. En effet, la figure 3.1
tirée de Bédard et Fraser (2008) illustre ce fait.
Ce graphique présente les valeurs-p obtenues par différentes méthodes en fonction
des hypothèses H0 sur la variable β6 . Les astérisques décrivent la courbe obtenue par
une approximation fréquentiste, soit la valeur-p de troisième ordre. Les détails de cette
technique ne sont pas abordés dans le présent mémoire, mais le lecteur intéressé est in-
vité à consulter Fraser et Reid (1993). Dans le contexte présent, ces approximations sont
de l’ordre de O(n−3/2 ), soit de O(32−3/2 ), et sont incluses à titre de comparaison. La
ligne continue correspondant à la méthode DA démontre que cette dernière produit des
résultats qui concordent bien avec les valeurs de troisième ordre. D’un autre coté, les
courbes pointillées illustrent deux simulations indépendantes de l’approche RWMH qui
mènent à des résultats très différents. La différence marquée entre des simulations indé-
pendantes et la discordance entre la méthode RWMH et les deux précédentes indiquent
que cet algorithme versatile possède une vitesse de convergence médiocre.
Il est généralement accepté que le choix de la distribution instrumentale pour un
algorithme RWMH est crucial. Dans l’article Bédard et Fraser (2008), la matrice de
52
1.0
DA algorithm
0.8
RWM algorithm
Third order
0.6
p−value
0.4
0.2
0.0
−0.20 −0.15 −0.10 −0.05 0.00 0.05 0.10
H0
Figure 3.1 – Graphique des valeurs-p obtenues par la méthode DA (ligne pleine), la
méthode RWMH (lignes pointillées) ainsi que les approximations de troisième ordre
(astérisques) en fonction des hypothèses H0 sur β6 .
53
covariance de la distribution instrumentale prenait la forme suivante

 
1 ... 0
 . . 
 
0, 0001  .. . . . ..  .
 
0 ... 1
La matrice identité multipliée par un facteur est souvent utilisée comme première
approche dans un tel contexte dû à la simplicité de générer des valeurs aléatoires de
la distribution instrumentale qui s’en suit. Cependant, cette méthode ne converge pas
de façon satisfaisante comme démontré par le graphique précédent. À titre illustratif,
la figure 3.2 démontre l’autocorrélation des composantes, obtenue grâce à la fonction
CODA du progiciel R. Une période de burn-in de 10 000 a été utilisée et 40 000 itérations
ont ensuite été analysées.
Figure 3.2 – Graphique de l’autocorrélation des valeurs générées par la méthode RWMH.
Donc, l’autocorrélation de chaque composante ne s’estompe pas au fur et à mesure

que l’intervalle augmente, ce qui corrobore les conclusions précédentes au sujet de la
piètre convergence de cette méthode.
54
Une manière couramment employée pour tenter d’améliorer la vitesse de conver-

gence est d’optimiser le facteur multipliant la matrice de covariance de la densité ins-
trumentale. Une notion déjà utilisée auparavant, la distance de sauts carrée moyenne
(DSCM) peut servir de critère d’optimisation. L’idée est de déterminer la valeur de la
variance instrumentale maximisant la DSCM et ensuite d’utiliser cette variance opti-
male dans les simulations subséquentes. La relation entre la DSCM et la variance instru-
mentale est illustrée dans la figure 3.3 et la relation similaire du taux d’acceptation est
également présentée dans la figure 3.4.
0,012
0,01
0,008
DMSC 0,006
0,004
0,002
0
0,05
0,1
0,15
0,2
0,25
0,3
0,35
0,4
0,45
0,5
0,55
0,6
0,65
0,7
0,75
0,8
0,85
0,9
0,95
1,00E-07
6,00E-07
2,00E-06
7,00E-06
3,00E-05
8,00E-05
4,00E-04
9,00E-04
5,00E-03
1,00E-02
Facteur multiplicatif
Figure 3.3 – Graphique de la DSCM en fonction du facteur multiplicatif de la matrice de

covariance (algorithme RWMH, exemple 2).
La distance de saut carrée moyenne semble être maximale sur l’intervalle

[10−4 ; 10−3 ] et la valeur initialement choisie dans l’article Bédard et Fraser (2008) ap-
partient à cet intervalle. Le taux d’acceptation correspondant semble raisonnable, soit
entre [0, 15; 0, 35]. Il est à noter que le taux d’acceptation tend vers 1 lorsque le fac-
teur multiplicatif tend vers 0. Tel que mentionné au chapitre 1, la méthode RWMH avec
55
1,00E+00
9,00E-01
8,00E-01
7,00E-01
6,00E-01
Taux
5,00E-01
d'acceptation
4,00E-01
3,00E-01
2,00E-01
1,00E-01
0,00E+00
1,00E-07
7,00E-07
4,00E-06
1,00E-05
7,00E-05
4,00E-04
1,00E-03
7,00E-03
0,03
0,09
0,15
0,21
0,27
0,33
0,39
0,45
0,51
0,57
0,63
0,69
0,75
0,81
0,87
0,93
Facteur multiplicatif
Figure 3.4 – Graphique du taux d’acceptation en fonction du facteur multiplicatif de la

matrice de covariance (algorithme RWMH, exemple 2).
distribution instrumentale normale est un cas d’algorithme symétrique et la probabilité

d’acceptation se simplifie à
π(y)
min 1, .
π(x)
Par conséquent, plus la variance de la distribution instrumentale est petite, plus un
π(y)
état proposé y a tendance à être proche d’un état actuel x et donc π(x) ≈ 1. En d’autres
mots, une faible variance de la distribution instrumentale mène à un fort taux d’accepta-
tion de sauts qui sont en moyenne très petits. Par conséquent, il est attendu que le gra-
phique du taux d’acceptation ne possède pas de maximum, puisque cette quantité peut
être arbitrairement proche de 1. D’un autre côté, une variance trop forte mène au rejet de
la plupart des sauts proposés puisqu’ils se trouveront généralement dans une région de
faible densité de la distribution cible. Donc, pour de grandes variances, le taux d’accepta-
tion tend vers 0. Dans le cas d’un algorithme RWMH avec densité instrumentale normale
et densité cible n-dimensionnelle avec composantes indépendantes et identiquement dis-
56
tribuées (i.i.d.), il a été montré que le taux d’acceptation optimal asymptotique (AOAR)
est de 0,234 lorsque n → ∞ (Roberts et al. (1997)). Bien que prouvé en présumant la
condition i.i.d., il apparaît que ce taux optimal est assez robuste et peut s’appliquer dans
des contextes déviant de cette supposition (Roberts et Rosenthal (2001)). En général,
il est donc courant d’ajuster la variance instrumentale afin d’obtenir un taux d’accepta-
tion d’environ 0,25. Pour l’exemple en cours, l’utilisation d’une variance maximisant la
DSCM mène à un taux d’acceptation raisonnablement près de 0,25, mais l’algorithme
reste sous-optimal.
Une autre manière d’améliorer la convergence d’un algorithme RWMH est de modi-
fier la forme de la matrice de covariance instrumentale afin de la faire correspondre à la
forme de la matrice de covariance cible. La motivation derrière cette technique est à la
fois empirique (voir Tierney (1994)) et théorique puisqu’il a été montré que la matrice de
covariance instrumentale asymptotiquement optimale est proportionnelle à la matrice de
covariance cible au moins dans le cas où les densités cible et instrumentale sont normales
(voir Roberts et al. (1997)). La seule difficulté réside dans le fait qu’aucune information
au sujet de la forme de cette matrice n’est connue à priori.
Une méthode facilement implémentable afin de remédier à ce problème est d’utiliser
l’estimateur empirique de la variance pour un échantillon généré par l’algortihme précé-
dent. Bien que l’approche RWMH avec densité instrumentale normale aux composantes
indépendentes ne donne pas de résultats convaincants, il est possible que la variance em-
pirique des valeurs générées soit tout de même suffisamment près de la variance cible
pour améliorer la convergence d’un algorithme subséquent. Par conséquent, dix essais
indépendants ont été simulés avec une période de chauffe de 10 000 et avec un nombre de
valeurs générées de 40 000 par essai. Les matrices de variance/covariance obtenues pour
chaque simulation ont été comparées entre elles et une forme approximative moyenne
57
de la variance, qu’on peut appeler A, a été construite :

 
1 0 0 0 0 0 0 0
 
 0 0, 01 0 0 0 0 0 0
 

 
 
 0 0 0, 5 0 0 0 0 0 
 
 
 0 0 0 1 0 0 0 0 
A= 
.

 0 0 0 0 2 0 0 0 
 
 
 0 0 0 0 0 1 0 0 
 
 
 0 0 0 0 0 0 1 0 
 
0 0 0 0 0 0 0 0, 1
Cette matrice a ensuite été employée comme matrice de covariance pour la distri-
bution instrumentale normale d’un nouvel algorithme RWMH. Le facteur multiplicatif
a été optimisé par un exercice identique à celui appliqué à la première méthode de ce
chapitre. Cette nouvelle approche a été simulée avec les mêmes conditions décrites au
chapitre 2, c’est-à-dire avec une période de chauffe de 10 000 et 4 000 000 valeurs gé-
nérées divisées en séries de 950. Le tableau 3.I rapporte aussi les mêmes mesures de
convergence qui ont été présentées à la section 2.3.
Le nouvel algorithme augmente la distance de saut carrée moyenne de plus de 50
fois et semble diminuer l’écart-type de la valeur-p pour l’hypothèse β6 = −0.1. Pour des
raisons de cohérence, cette hypothèse sera maintenue à travers ce chapitre pour toutes
les méthodes qui seront présentées.
En somme, cet ajustement de la matrice de covariance a permis d’améliorer quelque
peu la convergence de l’algorithme, mesurée selon les critères cités. Cependant, la
valeur-p demeure éloignée de l’approximation de troisième ordre, soit 0,75283. À pré-
sent, il reste à savoir s’il est possible de faire encore mieux et surtout à quel coût.
58
3.2 Algorithme Metropolis adaptatif
Une approche qui repose sur les mêmes fondements que l’idée précédente est l’algo-
rithme Metropolis adaptatif (AM) développé par Haario et al. (2001). La particularité de
cette méthode est le fait que la matrice de covariance est mise à jour à chaque temps (ou
saut) de la chaîne de Markov. Pour l’exemple en cours, à chaque temps j, une nouvelle
valeur y j+1 est proposée d’une distribution instrumentale N(x j ,C j+1 ), où C j+1 est calcu-
lée à partir des valeurs antécédentes x0 , x1 , x2 , . . . , x j . Cette nouvelle valeur est acceptée
avec probabilité :
π(y j+1 )
α(x j , y j+1 ) = min 1, .
π(x j )
Il est à noter que cet algorithme, même si sa probabilité d’acceptation est évocatrice
de certains algorithmes RWMH, n’est pas symétrique ni reversible. Toutefois, il existe
des résultats théoriques qui garantissent une convergence vers la distribution cible sous
certaines conditions qui seront mentionnées un peu plus loin. En général, il est préférable
de générer les premières j0 valeurs à l’aide d’une distribution instrumentale fixe et les
valeurs subséquentes à l’aide d’une distribution mise à jour selon la méthode décrite.
Donc, la forme de la matrice de covariance est donnée par

C,
0 j ≤ j0
Cj = .
 s Cov(x , . . . , x ) + s εI , j > j
n 0 j−1 n n 0
Le terme C0 représente la connaissance à priori de la matrice de covariance et j0

indique le temps à partir duquel l’algorithme adaptatif est appliqué. Le paramètre sn est
un facteur multiplicatif qui dépend de la dimension et qui peut être optimisé. Un élément
ε > 0 est introduit afin d’assurer la non singularité de la matrice C j et afin d’assurer la
convergence théorique de l’algorithme. En effet, ε > 0 ainsi qu’un support mesurable
borné S ⊂ Rn pour π sont nécessaires afin de démontrer, pour une fonction f mesurable
et bornée, que
1
Z
lim ( f (x0 ) + f (x1 ) + . . . + f (xN )) = f (x)π(dx) ,
N→∞ N + 1 S
59
presque toujours (voir Haario et al. (2001) pour de plus amples détails).
D’un côté pratique, ces deux conditions ne semblent pas être toujours cruciales
puisque l’algorithme converge de façon similaire sans ces restrictions dans tous les
exemples testés par les auteurs. Toutefois, afin de concilier théorie et pratique, il est
toujours possible de choisir S arbitrairement grand et ε arbitrairement petit par rapport à
S.
Finalement, il vaut la peine de mentionner que le calcul de la covariance ne requiert
pas autant de ressources que l’on pourrait l’imaginer à priori. La matrice de covariance
et la moyenne satisfont les formules récursives suivantes
j−1 sn
jx j−1 x0j−1 − ( j + 1)x j x0j + x j x0j + εIn

C j+1 = Cj +
j j
et
1
x j+1 = jx j + x j+1 .
j+1
L’emploi de ces équations élimine la nécessité de stocker les valeurs générées et per-
met de diminuer considérablement le coût de simulation. Les résultats de cette méthode
sont rapportés au tableau 3.I ; ceux-ci sont beaucoup plus convaincants que ceux obtenus
à la section 3.1, la valeur-p concordant avec l’approximation de troisième ordre. De plus,
la DSCM est augmentée de près de 100 000 fois et l’écart-type de la valeur-p est réduit
de moitié. Le paramètre sn a été choisi tel qu’il serait asymptotiquement optimal dans un
contexte où la densité cible est normale, c’est-à-dire égal à 2, 382 /8, tel que justifié dans
Tierney (1994) et Roberts et al. (1997).
Aussi, on remarque que la matrice de covariance (3.1) calculée par l’algorithme est
60
très différente de la matrice A utilisée précédemment :

 
675, 85 −8, 74 −12, 26 −5, 16 −2, 21 5, 75 −15, 11 −0, 53
 
 −8, 74 0, 11 0, 10 0, 07 0, 03 −0, 07 0, 20 0, 01 
 
 
 −12, 26 0, 05 −0, 01 −0, 11
 
0, 10 0, 74 0, 17 0, 02 
 
 −5, 16 0, 30 −0, 02 −0, 08
 
0, 07 0, 05 0, 20 0, 00 
  . (3.1)
 −2, 21 0, 03 −0, 01 −0, 02 0, 17 −0, 01
 
0, 08 0, 00 
 
5, 75 −0, 07 −0, 11 −0, 08 −0, 01 0, 10 −0, 16 −0, 01 
 

 
 −15, 11 0, 08 −0, 16
 
0, 20 0, 17 0, 20 0, 56 0, 01 
 
−0, 53 0, 01 0, 02 0, 00 0, 00 −0, 01 0, 01 0, 02
On constate, entre autres, la variance très grande de la première composante ainsi

qu’une covariance considérable entre la première et la plupart des autres composantes.
En rétrospective, il n’est pas surprenant que les deux premiers algorithmes démontrent
une convergence lente puisque la variance de la distribution instrumentale est vastement
différente de celle de la distribution cible.
3.3 Algorithmes Metropolis avec essais multiples
Bien que la méthode AM produise des résultats à priori satisfaisants, il peut être inté-
ressant de poursuivre l’analyse de cet exemple et d’examiner la performance de versions
spécifiques d’autres algorithmes locaux. Les méthodes suivantes seront computationnel-
lement plus intenses et demanderont davantage de ressources que celle du RWMH avec
densité instrumentale aux composantes indépendantes. Cependant, il serait intéressant
de voir si un meilleur compromis entre efficacité et effort computationnel peut être at-
teint. En plus, il serait intéressant de déterminer si d’autres algorithmes peuvent agir à
titre de remplacement aux méthodes adaptatives et si ces dernières sont véritablement
nécessaires dans ce contexte.
61
3.3.1 Algorithme MTM
Un algorithme de type marche aléatoire à essais multiples (Multiple-Try Metropolis

Algorithm (MTM)) sera la première approche utilisée. Cet algorithme a été développé
par Liu et al. (2000) et continue d’être un sujet d’intérêt en recherche MCMC. L’idée
fondamentale de cette méthode est d’élargir l’ensemble des valeurs proposées, c’est-à-
dire proposer plus d’un candidat par itération, afin d’améliorer l’exploration de l’espace.
Étant donné la valeur présente x, l’approche MTM génère k valeurs, {y1 , . . . , yk },
de la distribution instrumentale T (x, ·). Ensuite, un candidat Y = y est choisi parmi les
k valeurs proposées avec probabilité proportionnelle à un poids w(yl , x). Ce poids est
défini en général par
w(x, y) = π(x)T (x, y)λ (x, y) .
La densité instrumentale T (x, y) ne doit pas nécessairement être symétrique, mais

la condition T (x, y) > 0 ⇔ T (y, x) > 0 doit être remplie. La fonction λ (x, y) est une
fonction symétrique et nonnégative choisie par l’utilisateur avec condition T (x, y) >
0 ⇒ λ (x, y) > 0.
Subséquemment, k − 1 valeurs {x∗1 , . . . , x∗k−1 } sont générées de la distribution T (y, ·)
et x∗k = x. Finalement, y est accepté avec probabilité

∗ ∗ w(y1 , x) + . . . + w(yk , x)
α(x1 , . . . , xk , y1 , . . . , yk ) = min 1, .
w(x∗1 , y) + . . . w(x∗k , y)
Il est démontré dans Liu et al. (2000) que cette méthode avec les conditions citées
plus haut produit une chaîne de Markov réversible par rapport à la distribution cible.
Il existe plusieurs options pour le choix de λ (x, y) et cela donne naissance à plusieurs
types d’algorithmes MTM.
Une méthode MTM populaire se base sur le fait que si T (x, y) est symétrique, alors
λ (x, y) = 1/T (x, y) l’est aussi et donc que w(yl , x) = π(yl ) est un choix valide. Cela
donne lieu à la probabilité d’acceptation suivante

∗ ∗ π(y1 ) + . . . + π(yk )
α(x1 , . . . , xk , y1 , . . . , yk ) = min 1,
π(x∗1 ) + . . . + π(x∗k )
62
Cette dernière approche a été appliquée à l’exemple en cours avec k = 2 valeurs propo-
sées. Le tableau 3.I résume les résultats obtenus. Il apparaît que cet algorithme n’apporte
qu’une très légère amélioration par rapport à l’algorithme initial RWMH (x j , 0, 0001I8 ).
La DSCM est augmentée de trois fois et l’écart-type de la valeur-p est diminué de près
de 15%. Cependant, la valeur-p globale reste éloignée de l’approximation de troisième
ordre.
Il a été montré que le taux d’acceptation asymptotiquement optimal pour cet algo-
rithme se situe à 0,32 sous des conditions similaires à celles du taux optimal asympto-
tique de la méthode RWMH (voir Bédard et al. (2010b)). Dans le contexte présent, la
variance instrumentale a donc été ajustée afin de produire un taux d’acceptation global
se situant près de cette valeur.
3.3.2 Algorithme MTM hit-and-run
Il existe plusieurs autres types d’algorithmes MTM, dont une variante qui permet de
proposer plusieurs valeurs dépendantes à l’intérieur d’une même itération. En effet, il
peut arriver que la distribution instrumentale soit dépendante d’une variable auxiliaire,
R
disons e, telle que T (x, y) = Te (x, y) fx (e)de où fx (e) est la densité de cette variable.
Dans un tel cas, il est possible de générer un ensemble de k valeurs {y1 , . . . , yk } de la
distribution instrumentale Te (x, ·). Le candidat Y = y est choisi avec probabilité propor-
tionnelle à un poids we (yl , x). Ce poids est défini en général tel que
we (x, y) = π(x) fx (e)Te (x, y)λe (x, y) ,
où λe (x, y) est une fonction positive et symétrique.

Ensuite, k − 1 valeurs {x∗1 , . . . , x∗k−1 } sont générées à partir de la distribution Te (y, ·).
Le candidat y est accepté avec probabilité
( )
k−1
we (y, x) + ∑l=1 we (yl , x)
αe (x, x∗1 , . . . , x∗k−1 , y1 , . . . , yk ) = min 1, k−1
.
we (x, y) + ∑l=1 we (x∗l , y)
Cette forme du poids donné aux valeurs générées permet d’assurer la réversibilité de la
méthode (voir Liu et al. (2000) pour de plus amples détails).
63
Un algorithme appelé MTM hit-and-run (MTMHR) repose sur l’idée précédente et

tente d’améliorer l’exploration de l’espace de manière directionnelle. La première étape
est la génération d’une direction, c’est-à-dire un vecteur unité e à partir d’une distribution
f (e). Ensuite, un choix populaire est de générer (r1 , . . . , rk ) ∼ N(0, σ 2 ) et de poser yl =
x + rl e. Un candidat Y = y est choisi et enfin, k − 1 valeurs {x∗1 , . . . , x∗k−1 } sont générées
à partir de Te (y, ·) et l’algorithme est exécuté selon les étapes précédentes.
Il est à noter que si fx (e) = fy (e) = f (e), c’est-à-dire que la distribution du vecteur
est indépendante de x ou y, alors f (e) étant constante peut être exclue de la fonction
de poids. Aussi, si Te (x, ·) est symétrique alors λe (x, y) peut être choisie 1/Te (x, ·) et la
forme de la probabilité d’acceptation devient alors :
( )
k−1
π(y) + ∑l=1 π(yl )
α(x, x∗1 , . . . , x∗k−1 , y1 , . . . , yk ) = min 1, .
π(x) + ∑k−1 ∗
l=1 π(xl )
Cette derniere version de l’algorithme MTMHR avec k = 2 a été appliquée à

l’exemple en cours et la variance a été ajustée afin de produire un taux d’acceptation près
de la valeur asymptotiquement optimale de 0,46 présentée dans Bédard et al. (2010b).
Les résultats sont indiqués dans le tableau 3.I. Il apparaît que cet algorithme est pres-
qu’équivalent à l’algorithme initial RWMH(x j , 0, 0001I8 ). La DSCM est augmentée de
deux fois, mais l’écart-type de la valeur-p obtenue est identique sinon légèrement supé-
rieur. La valeur-p globale reste encore éloignée de l’approximation de troisième ordre.
3.3.3 Algorithme Metropolis avec rejet différé
Enfin, une méthode alternative aux approches de type MTM est l’algorithme avec
rejet différé (Delayed Rejection Metropolis-Hastings (DR)). Ce type d’algorithme a été
développé par Mira (2001) et repose sur un concept similaire à la méthode MTM. En
effet, l’objectif est toujours une meilleure exploration de l’espace par une génération
de valeurs multiples à l’intérieur d’une même itération. Cependant, les candidats sont
maintenant proposés de manière consécutive plutôt que simultanément. L’idée se base
sur le fait qu’un rejet suggère que le candidat proposé se trouve dans une région où π(x)
64
est faible. En utilisant cette information, une deuxième valeur est proposée, en général
d’une distribution moins étendue, avant d’incrémenter le temps. Dans ce type d’algo-
rithme (tout comme les méthodes MTM), il est clair que le taux d’acceptation doit tou-
jours être analysé avec prudence puisqu’il sera artificiellement élevé dû aux nombreuses
valeurs proposées avant un rejet.
Si l’on suppose qu’au temps j, X j = x, la première étape de la méthode DR est de
proposer une nouvelle valeur y1 à partir d’une distribution instrumentale q1 (x, dy1 ) et de
l’accepter avec la probabilité habituelle

π(y1 )q1 (y1 , x) N1
α1 (x, y1 ) = min 1, = min 1, .
π(x)q1 (x, y1 ) D1
Si la valeur y1 est acceptée, le temps est incrémenté et X j+1 = y1 . Si elle est rejetée,
une nouvelle valeur y2 est proposée d’une distribution q2 (x, y1 , dy2 ) et elle est acceptée
avec probabilité

π(y2 )q1 (y2 , y1 )q2 (y2 , y1 , x)[1 − α1 (y2 , y1 )]
α2 (x, y1 , y2 ) = min 1,
π(x)q1 (x, y1 )q2 (x, y1 , y2 )[1 − α1 (x, y1 )]

N2
= min 1, .
D2
Cette forme de la probabilité d’acceptation est une manière possible de préserver la

réversibilité de la chaîne (voir Tierney et Mira (1999) pour plus de détails). On remarque
que si y1 est rejeté cela implique que N1 < D1 et donc α1 (x, y1 ) = N1 /D1 ; la deuxième
probabilité d’acceptation se simplifie donc à

N2
α2 (x, y1 , y2 ) = min 1,
q (x, y1 , y2 )[π(x)q1 (x, y1 ) − π(y1 )q1 (y1 , x)]
2
N2
= min 1, .
q2 (x, y1 , y2 )[D1 − N1 ]
De façon générale, la i-ème valeur yi est proposée de la distribution qi (x, y1 , . . . , dyi )

si la valeur précédente yi−1 est rejetée. La probabilité d’acceptation pour la nouvelle
65
valeur devient

π(yi )q1 (yi , yi−1 )q2 (yi , yi−1 , yi−2 ) · · · qi (yi , yi−1 , . . . , x)
αi (x, y1 , . . . , yi ) = min 1,
π(x)q1 (x, y1 )q2 (x, y1 , y2 ) · · · qi (x, y1 , . . . , yi )

[1 − α1 (yi , yi−1 )][1 − α2 (yi , yi−1 , yi−2 )] · · · [1 − αi−1 (yi , . . . , y1 )]
[1 − α1 (x, y1 )][1 − α2 (x, y1 , y2 )] · · · [1 − αi−1 (x, y1 , · · · , yi−1 )]

Ni
= min 1,
Di
et
Di = qi (x, y1 , . . . , yi )(Di−1 − Ni−1 ) .
La grande utilité de ces formules récursives est naturellement de réduire le coût compu-
tationnel de la méthode.
Une manière simple d’appliquer l’algorithme DR est de proposer des valeurs in-
dépendantes des précédentes à chaque étape, c’est-à-dire qi (x, y1 , . . . , dyi ) = qi (x, dyi ).
Dans un tel cas, la variance de chaque distribution instrumentale pourrait être sélec-
tionnée telle que σ1 > σ2 > σ3 > . . . > σk , où k est le nombre d’essais permis. Bien
que relativement facile à implémenter, cette méthode a été montrée équivalente à un al-
gorithme RWMH lorsque la dimension de la distribution cible est grande (voir Bédard
et al. (2010a)). De façon intuitive, il apparaît qu’un tel algorithme n’utilise pas toute l’in-
formation contenue dans le rejet des variables passées. Une variante de cet algorithme
permet justement de corriger cet aspect par une approche directionnelle. La méthode
DR avec distribution instrumentale antithétique (delayed rejection with antithetic pro-
posal (DR anti)) propose un maximum de deux candidats dans des directions opposées.
En d’autres mots, le rejet du premier candidat est considéré comme une indication qu’il
se situait dans une région de faible densité de π(x) et qu’une région de forte densité
se trouve fort probablement dans la direction opposée. Bien que plusieurs choix de dis-
tributions instrumentales sont possibles, la présente analyse considérera la distribution
instrumentale normale seulement.
Comme d’habitude, on suppose que X = x est la valeur présente. La première étape
propose une direction Z j+1 = z provenant d’une distribution N(0, In ). Par la suite, un pre-
66
(1)
mier candidat Y j+1 = y(1) est obtenu en posant y(1) = x j +σ1 z, où σ1 > 0. La probabilité
d’acceptation associée est identique à celle d’un algorithme symétrique, c’est-à-dire
 
π y(1)
α1 x, y(1) = min 1, .
π (x)
(2)
En cas de rejet seulement, une deuxième valeur Y j+1 = y(2) est proposée, y(2) =
x j − σ2 z où σ2 > 0. Donc, le deuxième candidat se trouve effectivement dans la direc-
tion opposée au premier. Ensuite, afin de préserver la réversibilité de la chaîne, cette
deuxième valeur est acceptée avec probabilité
 h (2)
i 
2 )y −σ1 x/σ2 )
π(y(2) ) 1 − π((1+σ1 /σπ(y (2) )
α2 x, y(2) = min 1,
 h i +
,
π((1+σ1 /σ2 )x−σ1 y(2) /σ2 )
π(x) 1 − π(x) +
où [x]+ signifie max(0, x).
Les résultats de cette méthode appliquée à l’exemple en cours sont illustrés au tableau
3.I. La variance a été ajustée similairement aux méthodes précédentes afin de produire un
taux d’acceptation global près de la valeur 0,39, démontrée asymptotiquement optimale
dans Bédard et al. (2010a). Il apparaît que cet algorithme est aussi presqu’équivalent
à l’algorithme initial RWMH(x j , 0, 0001I8 ). La DSCM est augmentée de deux fois,
mais l’écart-type de la valeur-p obtenue en est légèrement supérieur. D’un autre côté,
la valeur-p globale se rapproche un peu plus de l’approximation de troisième ordre.
À première vue, il semble que les algorithmes MTM, MTMHR et DR anti n’offrent
pas une véritable amélioration par rapport à l’algorithme RWMH(x j , 0, 0001I8 ) initial.
Il est aussi nécessaire de se rappeler que les méthodes MTM et MTMHR nécessitent à
peu près deux fois plus de ressources puisqu’il faut générer deux valeurs par itération
et évaluer la densité cible à chacune d’entre elles. Si on suit la convention de diviser les
mesures de convergence par le nombre d’essais, elles représentent en fait des méthodes
dont la performance est inférieure à celle de la méthode RWMH pour cet exemple. D’un
autre côté, la méthode DR anti ne requiert qu’une simulation par itération et un facteur
67
Tableau 3.I – Résumé des résultats de méthodes locales (exemple 2) (Période de chauffe
10 000, Itérations 4 000 000).
Algorithme valeur-p pour β6 = −0, 1
RWMH - Instrumentale N(x j , 0, 0001I8 ) 0,89338

{écart-type de la valeur-p} (0, 12942)
{DSCM} 0,000247
{Taux d’acceptation} 0,345
RWMH - Instrumentale N(x j , 0, 005A) 0,84919

{DSCM} 0,011210
AM - Instrumentale N(x j ,C j )(ε = 0, 001, j0 = 40 000, sn = 0, 7) 0,75339

{DSCM} 29,97311
MTM 2 essais - Instrumentale N(x j , 0, 0003I8 ) 0,90904

{DSCM} 0,000737
MTMHR 2 essais - (r1 , r2 ∼ N(0, 0, 036) et fe (e) uniforme) 0,85552

{DSCM} 0,000571
DR anti - (Z j+1 ∼ N(0, I8 ), σ1 = σ2 = 0, 013) 0,79176

{DSCM} 0,000488
DA 0,75683
{DSCM} 1322,69
Approximation de trosième ordre 0,75283

68
de 2 semble assez exagéré dans cette situation, quoiqu’une probabilité d’acceptation plus
complexe doit être évaluée en cas de rejet. De toute manière, ces trois algorithmes locaux
semblent être moins performants que la méthode AM et encore moins performants que
la méthode DA.
3.4 Conclusion et analyse
Pour des raisons illustratives, il peut être intéressant de revenir à l’algorithme global
DA dont la performance semblait très satisfaisante. Il est possible d’obtenir une matrice
de covariance empirique à partir de valeurs générées par cette méthode. En effet, en
utilisant 40 000 itérations, un calcul empirique mène à une matrice B :
 
978, 32 −12, 66 −17, 65 −7, 33 −3, 23 8, 24 −21, 86 −0, 69
 
 −12, 66 0, 17 0, 15 0, 10 0, 04 −0, 10 0, 29 0, 01 
 
 
 −17, 65 0 −0, 16
 
0, 15 1, 04 0, 07 0, 25 0, 03 
 
 −7, 33 0, 42 −0, 04 −0, 11
 
0, 10 0, 07 0, 28 0 
B=   .(3.2)
 −3, 23 0 −0, 04 0, 24 −0, 01

0, 04 0, 11 0 
 
8, 24 −0, 10 −0, 16 −0, 11 −0, 01 0, 14 −0, 23 −0, 01 
 

 
 −21, 8 0, 11 −0, 23
 
0, 29 0, 25 0, 28 0, 81 0, 02 
 
−0, 69 0, 01 0, 03 0 0 −0, 01 0, 02 0, 03
On remarque tout d’abord que la forme de la matrice de covariance empirique gé-

nérée par l’algorithme AM, la seule méthode utilisée qui possédait une convergence
raisonnable, est comparable à B. Par conséquent, il est fort probable qu’aucune méthode
locale ne converge de manière satisfaisante sans une étape adaptative. En d’autres mots,
afin d’utiliser une méthode locale quelconque, il apparaît nécessaire pour cet exemple
d’avoir au préalable une estimation fiable de la matrice de covariance. En retrospective,
la piètre convergence des algorithmes locaux utilisés à la section 3.3 n’est donc pas sur-
prenante, car pour cet exemple, la matrice de covariance de la distribution cible B est
vastement différente de I8 .
69
Le tableau 3.II présente un résumé de certaines des méthodes locales présentées plus
tôt, cette fois avec matrice de covariance instrumentale B et ajustées afin d’approxima-
tivement obtenir un taux d’acceptation asymptotiquement optimal. On note, pour toutes
les méthodes, une amélioration remarquable des mesures de convergence utilisées. La
valeur-p globale est très près de l’approximation de troisième ordre, la DSCM est aug-
mentée de près de 700 000 fois et l’écart-type de la valeur-p est diminué de près de 6
fois. L’algorithme MTM, après l’ajustement de ses mesures d’efficacité par un facteur
de 2, semble identique à l’algorithme RWMH. La méthode DR anti apparaît légèrement
supérieure quant à elle, puisque dans ce cas, un facteur de 2 semble démesuré. L’effi-
cacité des méthodes locales semble rester inférieure à celle de la méthode DA, même
après l’ajustement par la matrice de covariance B. Cependant, l’écart paraît petit à toutes
fins pratiques et il est anticipé qu’un nombre légèrement plus élevé d’itérations pour la
méthode RWMH ou une méthode avec plus de 2 essais multiples mènent à des mesures
de convergence semblables à la méthode DA.
En conclusion, cet exemple démontre bien l’importance du choix de la distribution
instrumentale pour une méthode locale. Bien que versatiles, ces méthodes peuvent être
très sous-optimales surtout quand la variance de la distribution instrumentale est loin
de celle de la distribution cible. Pour cet exemple, l’idée initiale était d’examiner si des
méthodes locales pouvaient être aussi performantes que la méthode DA. En réalité, il
apparaît que la méthode RWMH(x j , 0, 0001I8 ), une des méthodes les plus versatiles,
ne converge pas de façon satisfaisante. Une légère modification telle que l’algorithme
RWMH(x j , 0, 005A) ou des méthodes plus avancées telles que le MTM, MTMHR ou
DR anti n’améliorent pas grandement la performance. L’algorithme AM est le seul qui,
grâce à sa nature adaptative, semble être en voie de convergence. À la lumière de ces
résultats, il semble qu’une étape adaptative est cruciale dans cet exemple, un fait qui est
démontré par la matrice de covariance (3.2) obtenue à partir de valeurs générées par l’al-
gorithme convergent DA. Dans cette situation, une estimation efficace (contrairement à
la matrice A qui a été obtenue par une méthode sous-optimale) de la matrice de cova-
70
Tableau 3.II – Résumé des résultats de méthodes locales (exemple 2) avec ajustement de
la matrice de covariance (Période de chauffe 10 000, Itérations 4 000 000).
Algorithme valeur-p pour β6 = −0, 1
RWMH - Instrumentale N(x j , 0, 7B) 0,75502

{DSCM} 139,17
MTM 2 essais - Instrumentale N(x j , B) 0,75675

{DSCM} 224,86
DR anti - (Z j+1 ∼ N(0, B), σ1 = σ2 = 0, 95) 0,75573

{DSCM} 270,73
DA 0,75683
{DSCM} 1322,69
Approximation de trosième ordre 0,75283

71
riance semble nécessaire en un premier lieu. Une fois cette estimation obtenue, des mé-
thodes locales RWMH, MTM et DR anti convergent de façon raisonnable. D’entre elles,
l’algorithme DR anti semble la méthode la plus performante et mène à des mesures de
convergence près de celles de la méthode DA. En plus d’explorer une comparaison entre
la méthode DA et certaines méthodes locales, ce chapitre visait aussi à présenter une
explication pour la piètre performance de l’algorithme RWMH(x j , 0, 0001I8 ) dans cet
exemple, un problème qui avait été laissé ouvert dans Bédard et Fraser (2008).
Une question reste toutefois sans réponse, soit comment obtenir une estimation fiable
de la matrice de covariance de la distribution cible. Précédemment, B a été calculée à
des fins d’illustration et il est clair qu’en pratique on ne dispose pas d’un algorithme déjà
convergent pour calculer cette matrice. Une réponse à cette question peut être d’utiliser
l’inverse de la matrice hessienne comme approximation de par sa relation à la matrice
d’information. Cette méthode est bien sûr vastement dépendante de la densité cible.
Une autre réponse peut être l’algorithme AM qui souvent peut donner rapidement une
estimation suffisamment fiable de la matrice de covariance. Il est possible qu’ensuite cela
soit plus avantageux, en termes de coût computationnel, d’utiliser cette estimation dans
un algorithme de type RWMH que de laisser poursuivre l’algorithme AM. Une autre
possibilité est de créer un algorithme qui combine l’étape adaptative à une méthode DR
ou MTM.
De toute manière, une méthode locale semble ardue dans le contexte de cet exemple.
Pour y arriver, un coût computationnel et analytique considérable est nécessaire. Par
conséquent, il est plus avantageux d’utiliser une méthode globale telle que le DA.
D’un autre côté, dans un contexte où une distribution instrumentale générique (telle
que N(0, In )) est déjà suffisamment adéquate et qu’aucune estimation de la matrice de
covariance est nécessaire, une méthode locale peut être à meilleur marché.
Par exemple, le tableau 3.III résume les mêmes mesures d’efficacité dans le contexte
de l’exemple de la section 2.3.1. Dans ce cas, les méthodes locales se comportent déjà
de façon suffisamment satisfaisante en utilisant une matrice de covariance instrumentale
72
Tableau 3.III – Résumé des résultats de méthodes locales (exemple 1) (Période de

chauffe 10 000, Itérations 4 000 000).
Algorithme valeur-s pour β = 1
RWMH - Instrumentale N(x j , 0, 3I3 ) 0,10794

{écart-type de la valeur-s} (0, 01483)
{DSCM} 0,123851
MTM 2 essais - Instrumentale N(x j , 0, 4I3 ) 0,10837

{DSCM} 0,194470
DR anti - (Z j+1 ∼ N(0, I3 ), σ1 = σ2 = 0, 6) 0,10851

{DSCM} 0,233924
DA 0,10764
{DSCM} 0,744361
73
I3 . Elles sont légèrement moins efficaces que la méthode DA et d’une certaine manière,
cette situation est semblable à celle décrite dans les tableau 3.II. Les méthodes globales
peuvent généralement posséder une vitesse de convergence supérieure aux méthodes
locales (qui en fait peuvent être tout au plus géométriquement ergodiques et ce seulement
sous certaines conditions particulières (voir Mengersen et Tweedie (1996))). Toutefois,
la facilité d’application de ces méthodes locales et leur versatilité sont souvent leurs
principaux attraits.
Dans le cas de l’exemple de la section 2.3.1, il peut être avantageux d’utiliser une mé-
thode locale, puisque le temps d’exécution ainsi que le temps d’implémentation (temps
de codage) en seront inférieurs à ceux de la méthode DA et donc, pour une précision don-
née, un nombre d’itérations suffisamment élevé sera rapidement atteignable. D’un côté
pratique, il peut être préférable d’opter pour un algorithme légèrement moins efficace
mais plus facilement implémentable puisque le temps de simulation est généralement à
meilleur marché que le temps de l’utilisateur.
De toute manière, la distribution cible dictera dans chaque cas la méthode appropriée,
mais il est clair que l’algorithme DA tout comme une méthode locale peuvent constituer
un choix optimal selon le contexte.
CONCLUSION
En première partie, ce mémoire a effectué un survol des algorithmes MCMC ainsi

que des grands théorèmes et résultats portant sur leur convergence. Il est à espérer que
ce chapitre a permis au lecteur d’apprécier l’importance de ces méthodes, leur facilité
d’application ainsi que les raisons de leur popularité dans la communauté scientifique et
surtout bayésienne.
Au deuxième chapitre, il a été question d’un algorithme Metropolis avec ajuste-
ment directionnel (DA) récemment développé dans Bédard et Fraser (2008). L’idée de
cette méthode était de combiner la versatilité de l’approche Metropolis-Hastings de type
marche aléatoire (RWMH) avec la performance supérieure de l’échantillonneur indé-
pendant (IS). Pour ce faire, la méthode DA employait une densité instrumentale Student
centrée au mode de la densité cible avec des degrés de liberté distincts selon la direc-
tion. Ces degrés de liberté étaient choisis tels que les rapports à un point s∗ et au mode
des densités instrumentale et cible coïncidaient. En plus, les auteurs de la méthode DA
avaient choisi de faire correspondre la courbure des deux densités au mode par l’emploi
d’une reparamétrisation menant à des matrices hesiennes identiques.
De par sa nature, l’algorithme DA a un rendement efficace dans le cas de distribu-
tions lisses et unimodales. À priori, cette condition est peu restrictive puisqu’en pratique
plusieurs problèmes de ce type existent. Toutefois, une utilisation consciencieuse de la
méthode DA nécessite une première étape d’optimisation afin d’identifier le mode et de
s’assurer qu’il est unique. À cette fin, il existe des algorithmes préétablis dans la plupart
des progiciels de statistique, comme par exemple la fonction nlm ou optim du logiciel
R. Les deux exemples étudiés plus tôt dans ce mémoire portaient sur des problèmes de
régression. Dans ces cas, l’utilisateur dispose souvent de la valeur des estimateurs par
moindres carrés. Cela peut servir de point de départ dans un algorithme d’optimisation
et, cette valeur étant généralement près du maximum, la convergence est la plupart du
temps rapide. Cependant, lorsque ce type d’estimateur n’est pas disponible ou si la fonc-
75
tion de densité est complexe ou possède des extremums locaux, le problème devient
beaucoup plus difficile à cerner. Dans ces cas, il est possiblement plus avantageux de
considérer un algorithme local comme le RWMH.
La deuxième étape de l’algorithme DA, c’est-à-dire l’évaluation de la matrice hes-
sienne au mode, représente une difficulté tout aussi importante. Dans le logiciel R, il est
possible d’évaluer numériquement cette matrice à l’aide de la fonction fdHess. Selon
le contexte, cette technique peut être très instable. Par exemple, l’évaluation de la ma-
trice hessienne pour l’exemple de la section 2.3.2 a nécessité plusieurs ajustements avant
d’obtenir des résultats fiables. D’un autre côté, il est possible d’utiliser un logiciel ana-
lytique comme Maple afin d’éviter les difficultés de la différenciation numérique. Enfin,
même si l’ajustement des courbures semble bénéfique à priori, son impact véritable n’est
pas tout à fait clair. De futures études sont nécessaires afin de déterminer l’importance
de l’utilisation de matrices hessiennes identiques et surtout si les avantages conférés en
termes de convergence l’emportent sur les difficultés numériques potentielles.
Quoiqu’il en soit, un des objectifs du chapitre 2 était d’étudier l’évolution de cer-
√
taines mesures de convergence empiriques par rapport au point s∗ = λ n . En utilisant
deux exemples tirés de Bédard et Fraser (2008), l’efficacité de l’algorithme semblait
être optimisée pour des valeurs de λ entre 2 et 4. Sur cet intervalle, les mesures empi-
riques indiquaient une amélioration modeste mais soutenue de la performance. Ces deux
exemples semblaient indiquer que le choix de λ était important mais non pas capital. Ce-
pendant, le dernier exemple unidimensionnel a démontré que l’algorithme DA pouvait
dans certaines situations mener à une convergence fortement dépendante de λ . Dans ce
cas, un choix de λ = 1 résultait en un algorithme très sous-optimal tandis qu’un choix
de λ ≥ 42 représentait la situation idéale. En pratique, la meilleure approche est donc
de simuler l’algorithme DA en employant différentes valeurs de λ tout en surveillant le
comportement des critères de convergence.
Finalement, le troisième chapitre visait à comparer l’approche DA avec des algo-
rithmes locaux beaucoup plus versatiles et facilement implémentables. L’exemple de la
76
section 2.3.2 a été utilisé comme paradigme pour la majeure partie de la discussion. Dans
ce contexte, plusieurs méthodes locales ont été trouvées inefficaces, un fait qui avait été
démontré dans Bédard et Fraser (2008) pour l’algorithme RWMH. Pour cet exemple,
une étape adaptative semblait nécessaire avant toute méthode locale et donc l’approche
DA était préférable. D’un autre côté, dans le contexte de l’exemple de la section 2.3.1,
les méthodes locales semblaient répondre au problème de façon raisonnable sans ajuste-
ment de la matrice de covariance. Dans ce cas, une méthode locale pouvait constituer un
choix moins dispendieux en termes de temps d’exécution et d’implémentation.
En conclusion, l’algorithme DA est une alternative attrayante aux algorithmes
MCMC traditionnels. Bien qu’une attention particulière doit être portée quant au choix
du paramètre λ , cette méthode semble efficace, surtout en grande dimension. Elle com-
porte certains désavantages comme un temps d’exécution considérable et des difficul-
tés numériques potentielles, mais elle peut constituer un choix optimal dans certains
contextes.
BIBLIOGRAPHIE
Bédard, M., Douc, R. et Moulines, E. (2010a). Scaling analysis of delayed rejection

MCMC methods. Article soumis.
Bédard, M., Douc, R. et Moulines, E. (2010b). Scaling analysis of multiple-try MCMC

methods. Stochastic Processes and their Applications. À paraître.
Bédard, M. et Fraser, D. A. S. (2008). On a directionally adjusted Metropolis-Hastings

algorithm. International Journal Of Statistical Sciences, 9 (Special Issue), 33–57.
Bédard, M., Fraser, D. A. S. et Wong, A. (2007). Higher accuracy for bayesian and
frequentist inference : Large sample theory for small sample likelihood. Statistical
Science, 22, 301–321.
Chan, K. et Geyer, C. (1994). Discussion de l’article Tierney (1994). Annals of Statistics,

22, 1747–1758.
Cogburn, R. (1972). The central limit theorem for Markov processes. Dans Proceedings
of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, 485–512.
University of California Press.
Cowles, M. et Carlin, B. (1996). Markov chain Monte Carlo convergence diagnostics :a

comparative review. Journal of the American Statistical Association, 91, 883–904.
Cox, D. R. et Snell, E. J. (1981). Applied Statistics : Principles and Examples. Chapman

and Hall. 192p.
Damerdji, H. (1994). Strong consistency of the variance estimator in steady-state simu-

lation output analysis. Mathematics of Operations Research, 19, 494–512.
Flegal, J. M. et Jones, G. L. (2010). Batch means and spectral variance estimators in

Markov chain Monte Carlo. Annals of Statistics, 38, 1034–1070.
78
Fraser, D. et Reid, N. (1993). Ancillaries and third order significance. Rapport technique,
Université de Toronto.
Gelman, A., Carlin, J. B., Stern, H. S. et Rubin, D. B. (1995). Bayesian Data Analysis.
New York : Chapman.
Gelman, A. et Rubin, D. B. (1992). Inference from iterative simulation using multiple

sequences. Statistical Science, 7, 473–483.
Genz, A. (1972). An adaptive multidimensional quadrature procedure. Computer Phy-

sics Communications, 4, 11–15.
Geweke, J. (1992). Evaluating the accuracy of sampling-based approaches to the calcu-

lation of posterior moments. Dans Bayesian Statistics 4, 169–193. Oxford University
Press.
Geyer, C. J. (1992). Practical Markov chain Monte Carlo (with discussion). Statistical
Science, 7, 473–511.
Haario, H., Saksman, E. et Tamminen, J. (2001). An adaptive Metropolis algorithm.

Bernoulli, 7, 223–242.
Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their
application. Biometrika, 57, 97–109.
Kipnis, C. et Varadhan, S. R. S. (1986). Central limit theorem for additive functionals of

reversible Markov processes and applications to simple exclusions. Communications
in Mathematical Physics, 104, 1–19.
Liu, J. S., Liang, F. et Wong, W. H. (2000). The multiple-try method and local optimi-
zation in Metropolis sampling. Journal of the American Statistical Association, 95,
121–134.
79
Liu, J. S., Wong, W. H. et Kong, A. (1995). Covariance structure and convergence rate of
the Gibbs sampler with various scans. Journal of the Royal Statistical Society : Series
B, 57, 157–169.
Matthews, P. (1993). A slowly mixing Markov chain with implications for Gibbs sam-
pling. Statistics and Probability Letters, 17, 231–236.
Mengersen, K. et Tweedie, R. L. (1996). Rates of convergence of the Hastings and Me-

tropolis algorithms. Annals of Statistics, 24, 101–121.
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. et Teller, E. (1953).

Equations of state calculations by fast computing machines. Journal of Chemical
Physics, 21, 1087–1092.
Mira, A. (2001). On Metropolis-Hastings algorithms with delayed rejection. Metron,

LIX, 231–241.
Raftery, A. et Lewis, S. (1992). How many iterations in the Gibbs sampler ? Dans Baye-
sian Statistics 4, 763–773. Oxford University Press.
Roberts, G. O. (1999). A note on acceptance rate criteria for CLTs for Metropolis-
Hastings algorithms. Journal of Applied Probability, 36, 1210–1217.
Roberts, G. O., Gelman, A. et Gilks, W. R. (1997). Weak convergence and optimal sca-
ling of random walk Metropolis algorithms. Annals of Applied Probability, 7, 110–
120.
Roberts, G. O. et Rosenthal, J. S. (1997). Geometric ergodicity and hybrid Markov

chains. Electronic Communications in Probability, 2, 13–25.
Roberts, G. O. et Rosenthal, J. S. (2001). Optimal scaling for various Metropolis-

Hastings algorithms. Statistical Science, 16, 351–367.
80
Roberts, G. O. et Rosenthal, J. S. (2004). General state space Markov chains and MCMC
algorithms. Probability Surveys, 1, 20–71.
Tierney, L. (1994). Markov chains for exploring posterior disiributions (with discussion).
Annals of Statistics, 22, 1701–1762.
Tierney, L. et Mira, A. (1999). Some adaptive Monte Carlo methods for Bayesian infe-
rence. Statistics in Medicine, 8, 2507–2515.
Annexe I
Programmes R
# Programme R p o u r l a m é t h o d e DA , p o u r s∗=1 . . . 20 ( e x e m p l e 1 )
ADR <− f u n c t i o n ( BurnIN , I t n s , s1 , s2 , d ) {
# P é r i o d e de c h a u f f e =10000 , i t é r a t i o n s =4M, s∗ minimum= s q r t ( 3 ) , s∗maximum=20∗ s q r t ( 3 ) , d= s q r t ( 3 )
S l<−s e q ( s1 , s2 , d )
# D é f i n i t i o n de l a d e n s i t é
P i i <− f u n c t i o n ( par ) {
x<−c ( −3 , −2 , −1 ,0 ,1 ,2 ,3)
y<−c ( − 2 . 6 8 , −4.02 , −2.91 , 0 . 2 2 , 0 . 3 8 , −0.28 , 0 . 0 3 )
exp (−7∗par [ 3 ] ) ∗prod ( c ( ( 1 + ( y−par [1] − par [ 2 ] ∗ x ) ^ 2 / 7 / exp ( 2 ∗par [ 3 ] ) ) ^ − 4 ) ) }
# D é f i n i t i o n du l o g de l a d e n s i t é
P i i l <− f u n c t i o n ( par ) {
x<−c ( −3 , −2 , −1 ,0 ,1 ,2 ,3)
y<−c ( − 2 . 6 8 , −4.02 , −2.91 , 0 . 2 2 , 0 . 3 8 , −0.28 , 0 . 0 3 )
l o g ( exp (−7∗par [ 3 ] ) ∗prod ( c ( ( 1 + ( y−par [1] − par [ 2 ] ∗ x ) ^ 2 / 7 / exp ( 2 ∗par [ 3 ] ) ) ^ − 4 ) ) ) }
# D é t e r m i n a t i o n du mode
MM<−optim ( c ( − 1 . 3 2 , 0 . 6 7 , 1 ) , P i i , method = "BFGS" , c o n t r o l = l i s t ( f n s c a l e = −1))

MMM<−a s . f u n c t i o n (MM[ 1 ] )
Xhat<−MMM( )
# C a l c u l de l a m a t r i c e h e s s i e n n e , de s a r a c i n e c a r r é e e t s a r a c i n e c a r r é e i n v e r s e
G<−f d H e s s ( Xhat , P i i l )
T<−a s . f u n c t i o n (G [ 3 ] )
H<−T ( )
pd1 <− pdSymm ( d i a g ( 3 ) )
m a t r i x ( pd1 ) <− −H
Hf<− p d M a t r i x ( pd1 , 0 . 5 )
H f i n v<− s o l v e ( Hf )
# I n i t i a l i s a t i o n de v a r i a b l e s
p<− ( 1 )
AvDOF<− ( 1 )
xvii
S v a l u e<− ( 1 )
D3d<− ( 0 )
D1d<− ( 0 )
ACCEP<− ( 1 )
StDOF<− ( 1 )
S t S v a l<− ( 1 )
StD3<− ( 1 )
StD1<− ( 1 )
# B o u c l e p o u r un s∗ donné
w h i l e ( s 2 +1>= s 1 ) {
p r i n t ( s1 )
# V a l e u r i n i t i a l e e t i n i t i a l i s a t i o n de v a r i a b l e s
VI<−c ( 1 , 1 , 1 )
j <− ( 1 )
SS<−rep ( 0 , I t n s )
F<− ( 1 : I t n s )
# S t a n d a r d i s a t i o n de l a v a l e u r i n i t i a l e e t c a l c u l du d e g r é de l i b e r t é o p t i m a l
X s t a r<−Hf %∗%( VI−Xhat )

U s t a r<− X s t a r / ( s q r t ( sum ( X s t a r ^ 2 ) ) )
R s q u a r e<−rep ( 1 , 5 0 )
Q s q u a r e<−rep ( 1 , 5 0 )
R s q u a r e<− R s q u a r e∗2∗ l o g ( ( P i i ( Xhat ) ) / ( P i i ( Xhat + H f i n v%∗% ( U s t a r ∗s 1 ) ) ) )
Q s q u a r e<− Q s q u a r e∗ ( t ( U s t a r ∗s 1 ) )%∗%( U s t a r ∗s 1 )
Min<− ( 1 : 5 0 )
Min<−abs ( ( ( Min +3 )∗ l o g ( 1+ Q s q u a r e / ( Min +3))) − R s q u a r e )
Minpt<−which . min ( Min )
F [ j ]<−Minpt
# I n i t i a l i s a t i o n d’ autres variables
Di3d<− ( 0 )
Di1d<− ( 0 )
ACPT<− ( 0 )
# Boucle pour l e s i t é r a t i o n s
w h i l e ( j < I t n s +BurnIN + 1 ) {
VIS<−VI
# P r o p o s i t i o n de v a l e u r e t c a l c u l du d e g r é de l i b e r t é o p t i m a l
R<− c ( rnorm ( 1 ) , rnorm ( 1 ) , rnorm ( 1 ) )

U s t a r N<− R/ s q r t ( sum (R^ 2 ) )
RsquareN<−rep ( 1 , 5 0 )
QsquareN<−rep ( 1 , 5 0 )
RsquareN<− RsquareN∗2∗ l o g ( ( P i i ( Xhat ) ) / ( P i i ( Xhat + H f i n v%∗% ( U s t a r N∗s 1 ) ) ) )

QsquareN<− QsquareN∗ ( t ( U s t a r N∗s 1 ) )%∗%( U s t a r N∗s 1 )
xviii
MinN<− ( 1 : 5 0 )
MinN<−abs ( ( ( MinN +3 )∗ l o g ( 1 + QsquareN / ( MinN+3))) − RsquareN )
MinptN<−which . min ( MinN )
F [ j + 1]<− ( MinptN )
# P r o p o s i t i o n d ’ une d i s t a n c e r a d i a l e
S<− r c h i s q ( 1 , MinptN , ncp = 0 )

SN<− ( MinptN + 3 ) ^ 0 . 5 ∗ s q r t ( sum (R^ 2 ) ) / s q r t ( S )
Y s t a r<−U s t a r N∗SN
# C a l c u l de a l p h a
Td<− f u n c t i o n ( x , y ) ( gamma ( ( y + 3) / 2 ) / p i ^ 1 . 5 / gamma ( y / 2 ) ) ∗ (

( 1 + ( sum ( x ^ 2 ) / ( y + 3 ) ) ) ^ −(( y + 3 ) / 2 ) ) ∗ ( y +3)^ −1.5
Num<− P i i ( Xhat + H f i n v%∗%Y s t a r ) ∗Td ( X s t a r , Minpt )
Den<− P i i ( Xhat + H f i n v%∗%X s t a r ) ∗Td ( Y s t a r , MinptN )
a l p h a<−min ( 1 ,Num / Den )
# A c c e p t a t i o n ou r e f u s
T<− r u n i f ( 1 )
i f ( T< a l p h a )
{ VI<−Xhat + ( H f i n v%∗%Y s t a r )
Minpt<−MinptN
X s t a r<− Y s t a r
ACPT<−ACPT+1
}
# D i s t a n c e e t v a l e u r −s
Di3d<− ( sum ( ( VIS−VI ) ^ 2 ) ) + Di3d

i f ( j >BurnIN ) {
i f ( VI [ 2 ] >=1) SS [ j −BurnIN ]<− ( 1 ) }
j <− j +1
}
# Fin b o u c l e i t é r a t i o n s
FN<− ( 0 )
SSN<− ( 0 )
D i s 3 d<− ( 0 )
D i s 1 d<− ( 0 )
JK<− I t n s / 1000
f o r ( g i n 1 : JK ) {
# S é p a r a t i o n en s é r i e s de 950
l o<−50+1000∗ ( g −1)
h i<−1000∗g−1
FN [ g ]<−mean ( F [ l o : h i ] )
SSN [ g ]<−mean ( SS [ l o : h i ] )
}
xix
# E n r e g i s t r e m e n t de s o r t i e s
AvDOF[ p ]<−mean ( FN )
S v a l u e [ p ]<−mean ( SSN )
StDOF [ p ]<−sd ( FN )
S t S v a l [ p ]<−sd ( SSN )
D3d [ p ]<− Di3d
ACCEP[ p ]<−ACPT
p<−p+1
s 1<−s 1 +d
}
# F i n b o u c l e s∗
# Sorties globales
par ( mfrow= c ( 2 , 2 ) )
p l o t ( Sl , AvDOF)
f o r ( mn i n 1 : l e n g t h ( S l ) ) {
s e g m e n t s ( S l [ mn ] , AvDOF[ mn] + StDOF [ mn ] , S l [ mn ] , AvDOF[ mn]−StDOF [ mn ] ) }
p l o t ( Sl , S v a l u e )
f o r ( mn i n 1 : l e n g t h ( S l ) ) {
s e g m e n t s ( S l [ mn ] , S v a l u e [ mn] + S t S v a l [ mn ] , S l [ mn ] , S v a l u e [ mn]− S t S v a l [ mn ] ) }
p l o t ( Sl , D3d / I t n s )
print ( " Svalue " )
print ( Svalue )
p r i n t ( " Dis3 " )
p r i n t ( D3d / I t n s )
p r i n t ( "AvDOF" )
p r i n t (AvDOF)
p r i n t ( " StDOF " )
p r i n t ( StDOF )
print ( " StSval " )
print ( StSval )
p r i n t ( "ACEPT" )
p r i n t (ACCEP)
}
xx
# Programme R p o u r l a m é t h o d e RWMH a v e c N( x j , 0 , 0 0 0 1 I 8 )
RWMH <− f u n c t i o n ( BurnIN , I t n s , Fac1 , AC) {
# P é r i o d e de c h a u f f e =10000 , I t n s =4M, Fac1 = 0 , 0 0 0 1 ,AC=d i a g ( 1 , 8 )

# I n i t i a l i s a t i o n de v a r i a b l e s
ACPT<−0
DIS3d<−0
A<− I t n s / 1000
Accep<−0
D i s<−0
j <−1
P<−0
SN<−0
# F a c t e u r d ’ a j u s t e m e n t ∗ M a t r i c e de c o v a r i a n c e
Sigma<−Fac1∗AC
# D é f i n i t i o n de l a d e n s i t é c i b l e
P i i <− f u n c t i o n ( par ) {
Func<− ( 0 )
P a r 2<−par [ 1 : 7 ]
for ( i in 1 :32) {
Func [ i ]<− dt ( exp ( par [ 8 ] ) ∗D[ i ] + XX[ i , ]%∗%P a r 2 , 4 )
}
prod ( Func ) ∗ exp ( 2 5∗par [ 8 ] )
}
# Valeur i n i t i a l e
VI<−c ( b0 , l o g ( s ) )
while ( j < I t n s + BurnIN + 1 ) {

VIS<−VI
# Saut proposé
Yn<− rmvnorm ( 1 , mean=VI , Sigma )
Num<− P i i ( Yn ) # N u m é r a t e u r du r a t i o a l p h a
Den<− P i i ( VI ) # D é n o m i n a t e u r du r a t i o a l p h a
a l p h a<− min ( 1 , ( Num / Den ) )
# A c c e p t a t i o n ou r e f u s
C<− r u n i f ( 1 ) ;
i f (C < a l p h a ) {
VI<−Yn
Accep<−Accep +1 }
xxi
# C a l c u l d i s t a n c e e t v a l e u r −p s i l ’ i t é r a t i o n e s t au−d e l à de l a p é r i o d e de c h a u f f e
i f ( j >BurnIN ) {
Tk<− ( VI [ 6 ] ) / ( s q r t (CK [ 6 , 6 ] ) ∗exp ( VI [ 8 ] ) )
TO<−( −0.0878 + 0 . 1 ) / ( s q r t (CK [ 6 , 6 ] ) ∗s )
i f ( Tk<TO) {P [ j −BurnIN ]<−1} e l s e {P [ j −BurnIN ]<−0}
D i s<−D i s + ( sum ( ( VI−VIS ) ^ 2 ) )
}
j <− j +1
}
# D i v i s i o n d e s s o r t i e s en s é r i e s de 950
f o r ( i i n 1 :A) {
SN [ i ]<−sum ( P [ ( 1 0 0 0 ∗ ( i −1) +51 ) : ( 1 0 0 0∗ i ) ] ) / 950

}
# C a l c u l v a l e u r −p moyenne e t é c a r t −t y p e
MN<−mean ( SN )
STDEV<−sd ( SN )
# Sorties globales
p r i n t ( " p r o b moyene " )

p r i n t (MN)
print ( " standard deviation " )
p r i n t (STDEV)
print ( " Accepations t o t a l " )
p r i n t ( Accep / ( I t n s +BurnIN ) )
print ( " distance t o t a l e " )
p r i n t ( Dis / I t n s )

Mireuta Matei 2011 Memoire

Transféré par

Droits d'auteur :

Formats disponibles

Mireuta Matei 2011 Memoire

Transféré par

Informations du document

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Mireuta Matei 2011 Memoire

Transféré par

Droits d'auteur :

Formats disponibles

Université de Montréal

Étude de la performance d’un algorithme Metropolis-Hastings avec ajustement

Département de mathématiques et de statistique

Mémoire présenté à la Faculté des études supérieures

c Matei Mireuta, 2011.

Étude de la performance d’un algorithme Metropolis-Hastings avec ajustement

a été évalué par un jury composé des personnes suivantes:

Pierre Lafaye de Micheaux, président-rapporteur

Mémoire accepté le: . . . . . . . . . . . . . . . . . . . . . . . . . .

TABLE DES MATIÈRES . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

LISTE DES TABLEAUX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

LISTE DES FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

LISTE DES ANNEXES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x

CHAPITRE 1 : INTRODUCTION AUX ALGORITHMES MCMC . . . . 1

1.6.5 Statistiques de convergence . . . . . . . . . . . . . . . . . . . 20

CHAPITRE 2 : ALGORITHME DE TYPE METROPOLIS-HASTINGS AVEC

CHAPITRE 3 : COMPARAISON ENTRE L’ALGORITHME DA ET DES

3.I Résumé des résultats de méthodes locales (exemple 2) (Période de

1.1 Graphique de la trace de paramètres en fonction des itérations. . . 16

2.1 Graphique de la densité cible (ligne pleine) et instrumentale (ligne

2.11 Graphique de l’écart-type σ de la valeur-p pour H0 : β6 = −0, 01

3.1 Graphique des valeurs-p obtenues par la méthode DA (ligne pleine),

Annexe I : Programmes R . . . . . . . . . . . . . . . . . . . . . . . . xvi

considérable d’algorithmes de simulation, comme entre autres, la méthode du rejet ou la

nues. Notre travail tentera donc d’élucider quelques aspects de la performance et de la

INTRODUCTION AUX ALGORITHMES MCMC

Les algorithmes MCMC représentent une méthode alternative pour l’échantillon-

1.1 Construction des algorithmes MCMC

Un espace X et une fonction de densité πu (x), possiblement non normalisée, mais

Conséquemment, si la chaîne est simulée assez longtemps, la distribution du j-ième

mentionnées à la section 1.3) et une première observation aléatoire Z1 = X j est alors

π(dx)P(x, dy) = π(dy)P(y, dx) ∀x, y ∈ X. (1.2)

= π(dy) P(y, dx) = π(dy) . (1.3)

Si B = 1, l’état y est accepté et devient le nouvel état de la chaîne (X j+1 = y). Si B = 0,

π(dx)P(x, dy) = c−1 πu (x)dx q(x, y)α(x, y)dy

qui est symétrique en x et y (voir Metropolis et al. (1953)).

ainsi qu’à la distribution instrumentale au-delà de ce chapitre introductoire.

1.2 Algorithmes MCMC élémentaires

Il existe une multitude de façons de construire la chaîne de Markov instrumentale,

(c) Échantillonneur Metropolis-Hastings indépendant (Independent Sampler (IS))

dans le calcul de valeurs-p.

L’algorithme GS peut s’implémenter de façon déterministe comme l’exemple ci-haut

1.3 Convergence des algorithmes MCMC

Ensemble, ces deux conditions garantissent la convergence d’une chaîne de Markov

Théorème 1. Une chaîne de Markov φ -irréductible et apériodique avec distribution

lim ||P j (x, ·) − π(·)|| = 0,

Ici, || · || dénote la distance de variation totale entre deux mesures quelconques

presque toujours si π(|h|) < ∞.

1.4 Taux de convergence des algorithmes MCMC

Bien que le théorème 1 garantisse la convergence asymptotique des chaînes MCMC,

||P j (x, ·) − π(·)|| ≤ Mρ j , j = 1, 2, 3, . . . ∀x ∈ X.

Définition 5. Un sous-ensemble C ⊆ X est dit petit s’il existe un entier positif n0 ∈ N,

Pn0 (x, ·) ≥ εν(·) x ∈ C

et donc, Pn0 (x, A) ≥ εν(A) pour tout x ∈ C.

De façon intuitive, cette définition implique que les probabilités de transition en n0

Les probabilités de transition satisfont

Théorème 3. Si une chaîne de Markov possède une distribution stationnaire π(·) et

Démonstration. Voir Roberts et Rosenthal (2004).

Dans ce mémoire, la propriété de l’ergodicité géométrique restera qualitative, mais