# Méthode de Monte-Carlo par chaînes de Markov

Les méthodes de Monte-Carlo par chaînes de Markov, ou méthodes MCMC pour Markov chain Monte Carlo en anglais, sont une classe de méthodes d'échantillonnage à partir de distributions de probabilité. Ces méthodes de Monte-Carlo se basent sur le parcours de chaînes de Markov qui ont pour lois stationnaires les distributions à échantillonner.
Certaines méthodes utilisent des marches aléatoires sur les chaînes de Markov (algorithme de Metropolis-Hastings, échantillonnage de Gibbs), alors que d'autres algorithmes, plus complexes, introduisent des contraintes sur les parcours pour essayer d'accélérer la convergence (Monte Carlo Hybride, Surrelaxation successive).
Ces méthodes sont notamment appliquées dans le cadre de l'inférence bayésienne.

## Approche intuitive

On se place dans un espace vectoriel $\varepsilon$ de dimension finie n. On veut générer aléatoirement des vecteurs $x$ suivant une distribution de probabilité π. On veut donc avoir une suite de $N$ vecteurs $(x_{0},x_{1},\dots ,x_{N-1})$ telle que la distribution des $x_{i}$ approche π.
Les méthodes de Monte-Carlo par chaînes de Markov consistent à générer un vecteur $x_{i}$ uniquement à partir de la donnée du vecteur $x_{i-1}$ ; c'est donc un processus « sans mémoire », ce qui caractérise les chaînes de Markov. Il faut donc trouver un générateur aléatoire avec une distribution de probabilité $q_{x_{i-1}}$ qui permette de générer $x_{i}$ à partir de $x_{i-1}$. On remplace ainsi le problème de génération avec une distribution π par $N$ problèmes de génération avec des distributions $q_{x_{i}}$, que l'on espère plus simples.

## Principe

On veut simuler une loi π sur un espace d'états général (Ω ; Ɛ). Une transition sur (Ω ; Ɛ) est une application P : (Ω ; Ɛ) → [0 ; 1] telle que P(·, A) pour tout A∈Ɛ est mesurable et P(x,·) pour tout x∈Ω est une probabilité sur (Ω ; Ɛ). Soit X = (X$_k$,k∈ℕ) une chaîne de Markov homogène de transition P et de loi initiale X$_0$ ~ v, on note P$_v$ la loi de la chaîne X.
Pour simuler π, on veut savoir construire une transition P telle que ∀ v, vP$^k$ → π, avec convergence en norme en variation totale ∥μ∥ = sup$_{A∈Ɛ}$ μ(A) − inf$_{A∈Ɛ}$ μ(A). La chaîne de transition P est ergodique.
Convergence d'une chaîne de Markov — Si une transition P est π-réductible, π-invariante, apériodique, Harris-récurrente, alors ∀x, P$^k$(x,·) →$_{k→∞}$ π.
La dernière condition, délicate, est satisfaite dans les cas de l'échantillonnage de Gibbs et de l'algorithme de Metropolis-Hastings.

## Méthodes d'échantillonnage

Les méthodes d’échantillonnage sont utilisées pour estimer la distribution (complète) postérieure des paramètres d’un modèle. Parmi elles, les méthodes de Monte-Carlo sont très précises. Néanmoins, elles sont coûteuses en calculs et prennent beaucoup de temps à converger. Elles sont également limitées par la taille de l’échantillon, puisqu’elles deviennent insolubles lorsque ceux-ci sont trop grands. Même sur de petits échantillons, le calcul d'une distribution de probabilité peut s’avérer difficile. Il existe plusieurs approches à ce problème, utilisées pour obtenir un ensemble de bons réseaux d'échantillonnage à partir de la distribution postérieure.
* Algorithme de Metropolis : l'algorithme de Metropolis consiste en une marche aléatoire avec une correction afin que chaque nouvelle valeur extraite puisse être rejetée si elle correspond à une valeur de la fonction d'intégration plus petite que le point précédent, avec une probabilité égale à un moins le rapport entre ces valeurs. Le processus aléatoire se concentre ainsi dans les zones où la fonction d'intégration a des valeurs plus élevées, d'où l'efficacité. Toutefois, contrairement aux méthodes de Monte-Carlo simples, où les valeurs échantillonnées sont statistiquement indépendantes, dans l'algorithme de Metropolis (et les chaînes de Monte-Carlo de Markov en général), elles sont auto-corrélées.
* Algorithme de Metropolis–Hastings : Hastings a publié une généralisation de l'algorithme Metropolis (qui a valu le nom d'algorithme Metropolis-Hastings) qui utilise un processus aléatoire de base (distribution des propositions) qui n'est pas nécessairement du type marche aléatoire.
* Echantillonnage de Gibbs : exige que toutes les distributions conditionnelles de la distribution cible soient connues correctement et puissent être utilisées pour les échantillonner. Cette méthode est devenue populaire car, si elle peut être appliquée, elle ne nécessite pas de paramètres d'algorithme prédéfinis.
* Échantillonnage par tranches : dépend du principe selon lequel il est possible d'échantillonner une distribution uniformément à partir de la région située sous la courbe de sa fonction de densité. Cette méthode alterne entre l'échantillonnage uniforme dans la direction verticale et l'échantillonnage de la couche identifiée par la position verticale actuelle.
* Metropolis à tentatives multiples : il s'agit d'une variante de l'algorithme de Metropolis-Hastings qui permet des tentatives multiples à chaque point. Cela permet à l'algorithme de définir des étapes plus importantes sous une forme systématique, ce qui permet de s'attaquer à des problèmes présentant des dimensions intrinsèquement importantes.
* Méthode de sur-relaxation : une version Monte Carlo de cette technique peut être considérée comme une variation de l'échantillonnage Gibbs ; elle évite parfois les répétitions de parcours.
* Méthode hybride/Hamiltonienne de Monte Carlo (HMC) : Implante la dynamique Hamiltonienne dans la chaîne de Markov en augmentant l'espace étudié avec un espace auxiliaire de temps. La densité cible devient alors la fonction d'énergie potentielle et le temps est graduellement tiré d'une distribution prédéterminée conditionnée au point xi , où une trajectoire d'énergie constante est calculée et arrêtée après un temps arbitraire. À la fin de l'échantillonnage, les données sont rejetées. Le résultat final est d'allonger les étapes proposées tout en les maintenant dans les régions à haute densité.
* No U-Turn Sampling (NUTS) : Variante de la méthode HMC qui permet d'éviter les répétitions en chaîne en calculant de manière adaptative le temps des trajectoires.
* La méthode MCMC de Langevin : La méthode MCMC de Langevin et d'autres méthodes basées sur la méthode du gradient (éventuellement en dérivée seconde) du logarithme de la distribution objective, proposent des chemins qui sont de préférence dans la direction de la densité de probabilité la plus élevée.
* Méthodes basées sur le changement de dimension : La méthode du saut réversible est une variante de la méthode de Metropolis-Hastings qui permet de proposer des pas qui modifient la dimension de l'espace paramétrique sur lequel la marche aléatoire opère. Cette méthode a été proposée en 1995 par le statisticien Peter Green de l'Université de Bristol. Les méthodes MCMC qui modifient la dimensionnalité sont utilisées depuis longtemps dans les applications de physique statistique, où les distributions basées sur l'ensemble grand-canonique (par exemple, lorsque le nombre de molécules dans une boîte est variable) ont été utilisées pour divers problèmes. Une certaine variation de la méthode du saut réversible est nécessaire dans le cas de l'échantillonnage MCMC ou Gibbs sur des modèles bayésiens non paramétriques tels que ceux impliqués dans le processus de Dirichlet (un processus stochastique qui peut être considéré comme une distribution de probabilités dont le domaine est lui-même une distribution aléatoire) ou celui du processus du restaurant chinois, où le nombre de composants, de grappes, etc. est automatiquement déduit des données.

## Construction de réseaux aléatoires et recherche de motifs

La méthode de Monte Carlo par chaines de Markov est souvent utilisée dans les algorithmes traitant les réseaux (biologiques ou non). Plusieurs applications sont possibles dont 2 principales : la génération de réseaux aléatoires et la classification d’éléments dans les graphes. Les exemples choisis pour illustrer l'utilisation des MCMC par la suite seront basées sur la biologie.

### Réseaux aléatoires

La MCMC est très souvent utilisée pour générer des réseaux aléatoires servant comme modèles nuls, le plus proche possible du hasard, tout en gardant les caractéristiques de bases d’un réseau étudié afin de rester comparable. Ces modèles nuls permettent ensuite de déterminer si des caractéristiques d’un réseau étudié sont statistiquement significatives ou non.
Le déroulement de la méthode est simple : 2 arêtes sont prises en compte (A -> B ; C -> D) et un échange des nœuds de ces arêtes (A -> D ; C -> B) est ensuite considéré et accepté seulement si les nouvelles arêtes n’existent pas déjà et qu’il n’y a pas d’arête cyclique (A -> A). Le nombre d’échanges considérés suit souvent la formule Q x E, où E est le nombre d’arêtes d’un réseau réel (souvent le réseau étudié) et Q est un nombre suffisamment élevé pour s’assurer de l’aspect aléatoire du réseau produit, souvent mis a 100.
L’utilisation de la méthode de Monte Carlo par chaînes de Markov pour générer un modèle nul contre lequel on détermine la significativité d’un (ou plusieurs) caractère(s) peut être rencontré sous le nom de «switching algorithm», avec le «matching algorithm» étant l’alternatif dans la génération des réseaux aléatoires. Ce dernier, qui n’emploie pas le MCMC, présente aussi comme désavantage la présence d’un biais dans les réseaux produits. Ces algorithmes sont le plus souvent utilisés dans la biologie pour rechercher des motifs dans les réseaux, des sous-graphes composés d’un nombre limité de nœuds et qui se retrouvent dans un très grand nombre de réseaux. Parmi les outils qui utilisent le MCMC pour générer des réseaux aléatoires sont mfinder, FANMOD, KAVOSH et NetMODE.

### Recherche de motifs

En ce qui concerne l’utilisation du MCMC pour la classification d’éléments dans un graphe, on parle notamment de «Markov clustering» (MCL), une méthode de classification non-supervisée basée sur le concept de «marche aléatoire», nécessitant presque aucune information a priori afin de pouvoir classifier les éléments d’un graphe. Plus spécifiquement, en partant d’un nœud, si on se «promène» de nœud en nœud en passant par les arêtes, on a tendance à passer plus souvent entre les nœuds qui se trouvent au sein d’un même groupe que de faire un passage uniforme entre tous les nœuds. Ainsi, en augmentant l’importance des arêtes fréquemment traversées et diminuant celle des arêtes moins traversées, les groupes dans un réseau sont progressivement mises à jour.

## Application à la biologie

Il existe deux principaux types d’utilisation du MCMC dans la biologie de systèmes, à savoir les réseaux de gènes et la dynamique moléculaire pouvant se résumer à un système de molécules. Dans ces deux cas, il s’agit de comprendre les interactions complexes entre plusieurs éléments. C'est pourquoi la méthode MCMC est combinée aux réseaux bayésiens.

### Réseaux bayésiens

Les réseaux bayésiens sont couramment utilisé en biologie de systèmes et associés à la méthode MCMC. Ils offrent une représentation nette et compacte de distribution de probabilités communes de systèmes. Leur utilisation dans des modèles graphiques présentent plusieurs avantages. Dans un premier temps, ils peuvent représenter les relations/dépendances causales entre les variables et ainsi être interprétés comme un modèle causal qui génère les données. De plus, les réseaux bayésiens sont utiles pour adapter les paramètres d'un modèle en apprentissage automatique, qu'il soit utilisé pour effectuer de la prédiction ou de la classification de données. La théorie des probabilités bayésiennes s'est également révélée efficace pour décrire l'incertitude et pour adapter le nombre de paramètres à la taille des données. La représentation et l'utilisation de la théorie des probabilités dans les réseaux biologiques permet à la fois de combiner les connaissances et les données d'un domaine, d’exprimer les relations de cause à effet, d’éviter de sur-ajuster un modèle pour former des données et de tirer des enseignements d'ensembles de données incomplets.

### Réseaux bayésiens dynamiques

La rétroaction est une caractéristique essentielle de nombreux systèmes biologiques. Ce qui fait de la création de mesures expérimentales de séries temporelles est un défi particulièrement important en modélisation des systèmes biologiques.
Les réseaux bayésiens sont adaptés à la modélisation de boucles de rétroaction, ainsi qu'aux séries temporelles. Ce sont les réseaux bayésiens dynamiques qui consistent en l'utilisation de réseaux bayésiens lors de la modélisation de séries temporelles ou de boucles de rétroaction. Dans ces conditions, les variables sont indexées par le temps et reproduites dans les réseaux. Des modèles de Markov cachés ainsi que des systèmes dynamiques linéaires sont des cas particuliers des réseaux bayésiens dynamiques. Dans les modèles de Markov cachés, il y a un ensemble caché de nœuds (normalement des états discrets), un ensemble de variables observées, et les tranches du graphique n'ont pas besoin d'être temporelles. Ils sont souvent utilisés pour l'analyse des séquences notamment dans le cas de réseaux de gènes.
Les réseaux bayésiens dynamiques ont été utilisés pour déduire des interactions de régulation génétique à partir de données de puces à ADN, à partir de quelques dizaines de points temporels issus d'un cycle cellulaire. Notamment en étant combiné au MCMC, cela a permis d'accéder aux performances d'inférence du réseau avec différentes tailles d'ensembles d'entraînement, d'antécédents et de stratégies d'échantillonnage.

### Réseaux de gènes

Les réseaux de gènes sont bien adaptés à la modélisation des systèmes biologiques complexes et stochastiques. Ils représentent chaque expression de gène par une variable décrivant la régulation entre les gènes.
Dans ce type de réseaux, l'inférence de leurs structures, par exemple, l'identification des réseaux de régulation et de signalisation à partir de données, est l'aspect le plus intéressant. La distinction de corrélation permet l'élucidation de dépendances entre les variables mesurées et ainsi à l'apprentissage de relations inconnues et d'excellents modèles prédictifs. Néanmoins, l'apprentissage des structures de modèles est difficile et nécessite souvent une conception expérimentale minutieuse.

### Dynamique moléculaire

La méthode classique pour simuler l'évolution de molécules en réaction dans le temps est basée sur l'hypothèse du continuum. Lorsque le nombre de ces molécules réagissant dans un ensemble de réactions est de l'ordre du nombre d'Avogadro, on peut supposer que cette concentration (nombre de molécules dans l’ensemble) est une variable réelle continue. Dans ce cas, la cinétique classique de l'action de masse peut être utilisée pour décrire les vitesses de réaction. Cependant, lorsque le nombre de ces molécules est de l'ordre de centaines ou de milliers, nous ne pouvons plus utiliser l'hypothèse du continuum. Par conséquent, au lieu de concentrations à valeur réelle, nous devons considérer le nombre de molécules à valeur entière. Un autre effet du faible nombre de molécules est que la cinétique classique de l'action de masse n'est plus valable. Les taux de réaction ne sont plus déterministes, et une approche probabiliste est donc nécessaire. Au lieu de tenir compte de la quantité de réactifs consommés et de produits fabriqués dans un intervalle de temps, nous devons tenir compte de la probabilité qu'une réaction se produise dans un intervalle de temps. Cette approche de modélisation des réactions est connue sous le nom de cinétique stochastique.

### Exemples en biologie de systèmes

Les processus cellulaires en biologie systémique présentent souvent un caractère aléatoire$^,$$^,$causé par le faible nombre de molécules en réaction.
* Des modèles cinétiques stochastiques$^,$ont été utilisés pour expliquer comment une population de cellules d'E. coli infectées pouvait se diviser en sous-populations suivant des chemins différents.
* Il a été démontré que la bifurcation phénotypique des cellules T infectées par le VIH-1 pouvait être expliquée par la présence de faibles concentrations moléculaires et d'une boucle de rétroaction transcriptionnelle stochastique.
* Hensel et al. ont développé un modèle détaillant la croissance intracellulaire d'un virus impliqué dans des stomatites pour démontrer que l'expression stochastique des gènes contribue à la variation du rendement viral.

### Estimation des paramètres par MCMC

La cinétique stochastique est devenue un élément de base pour la modélisation de divers phénomènes en biologie de systèmes. Ces modèles, bien plus encore que les modèles déterministes, posent un problème difficile dans l'estimation des paramètres cinétiques à partir de données expérimentales. Initialement, les paramètres ont été fixés à des valeurs biologiquement plausibles, puis ajustés à l'œil nu de sorte que la simulation du modèle ressemble aux données expérimentales. Dans ce cas, les estimations des paramètres peuvent être facilement obtenues par ajustement des moindres carrés ou par estimation du maximum de vraisemblance. Cependant l'utilisation des méthodes statistiques de Monte Carlo permettent le maintien de la stochasticité du modèle lors de l'estimation de ces paramètres. Ces méthodes d'estimation peuvent être classées en deux catégories bien connues : les méthodes de maximum de vraisemblance et les méthodes d'inférence bayésiennes.
Les méthodes d'inférence bayésiennes tentent d'obtenir la distribution postérieure des paramètres, qui peut ensuite être maximisée pour obtenir des estimations maximales a posteriori. La plupart des méthodes d'inférence bayésiennes s'appuient sur les techniques MCMC. L'application à la biologie de systèmes ne fait pas exception :
* Rempala et al. utilisent l'inférence bayésienne avec MCMC pour estimer les paramètres dans un modèle spécifique de transcription de gènes.
* Wilkinson et al. ont développé une série d'algorithmes d'inférence bayésienne, qui utilisent tous des méthodes MCMC$^,$$^,$$^,$$^,$.
* Autorégulation des gènes procaryotes$^,$$^,$.
* Modèles de transcription des gènes$^,$.
* Expression de l'ARN$^,$$^,$.
* Dynamique de l'ADNmt à cellule unique.