M. Khaled [KM1, KM2, 2008] établit une récurrence originale permettant un calcul récursif des lois marginales d'une loi temporelle jointe π de Z=(Z₁,Z₂,⋯,Z_{T})∈E^{T}, loi utilisée en économétrie pour modéliser un régime caché [G. Lovinson, 2006]. Une conséquence intéressante est de faciliter le calcul de la constante de normalisation de π. Nous généralisons ces résultats en observant que la loi π envisagée étant une loi de Gibbs sur T={1,2,⋯,T}, π est aussi un champ de Markov sur T dont il est aisé de manipuler les lois conditionnelles. Cette approche permet de généraliser la récurrence fondamentale de M. Khaled au cas d'un espace d'état E fini quelconque et d'une loi π de Gibbs générale, temporelle ou spatiale. Cette récurrence permet un calcul simplifié de la constante de normalisation de π et la performance numérique de cet algorithme est évaluée pour π un modèle d'Ising spatial

Guyon, Xavier

Hardouin, Cécile

HAL-Paris1

Chaine et champ de markov : quelques re´currences etcalcul d’une constante de normalisationXavier Guyon, Ce´cile HardouinTo cite this version:Xavier Guyon, Ce´cile Hardouin. Chaine et champ de markov : quelques re´currences et calculd’une constante de normalisation. 2009. <hal-00397201>HAL Id: hal-00397201https://hal.archives-ouvertes.fr/hal-00397201Submitted on 19 Jun 2009HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.Chaîne et champ de Markov : quelques récurrences et calculd’une constante de normalisationX. Guyon et C. HardouinSAMOS-MATISSE, Université Paris 1June 19, 2009AbstractM. Khaled [KM1, KM2, 2008] établit une récurrence originale permettant un calcul récur-sif des lois marginales d’une loi temporelle jointe π de Z = (Z1, Z2, · · · , ZT ) ∈ ET , loi utiliséeen économétrie pour modéliser un régime caché [G. Lovinson, 2006]. Une conséquence in-téressante est de faciliter le calcul de la constante de normalisation de π.Nous généralisons ces résultats en observant que la loi π envisagée étant une loi de Gibbssur T = {1, 2, · · · , T}, π est aussi un champ de Markov sur T dont il est aisé de manipulerles lois conditionnelles. Cette approche permet de généraliser la récurrence fondamentale deM. Khaled au cas d’un espace d’état E fini quelconque et d’une loi π de Gibbs générale,temporelle ou spatiale. Cette récurrence permet un calcul simplifié de la constante de nor-malisation de π et la performance numérique de cet algorithme est évaluée pour π un modèled’Ising spatial.Mots clés : modèle de Gibbs, chaîne de Markov, champ de Markov, loi marginale, con-stante de normalisation, calcul récursif, modèle d’Ising.Classification AMS : 60G60, 62M401 IntroductionSoit T un entier > 0 fixé, E = {e1, e2, · · · , eN} un espace d’états fini à N éléments, Z =(Z1, Z2, · · · , ZT ) une suite temporelle de loi jointe π sur ET . On suppose que π est une loi deGibbs d’énergie UT dérivant de potentiels de singletons (θs)s=1,T et de paires (Ψs)s=2,T , c’est àdire, notant z(t) = (z1, z2, · · · , zt),π(z(T )) = C expUT (z(T )) avec C−1 =Xz(T )∈ETexpUT (z(T )) oùUt(z(t)) =Xs=1,tθs(zs) +Xs=2,tΨs(zs−1, zs) si 2 ≤ t ≤ T ,et U1(z1) = θ1(z1). Généralement, le calcul de la constante de normalisation C est impossible dufait de la grande cardinalité de ET . L’objectif ici est d’obtenir des formules explicites permettantle calcul des lois marginales πt(z(t)) = πt(z1, z2, · · · , zt), 1 ≤ t ≤ T (avec π = πT ) ainsi que lecalcul de la constante de normalisation C.Le graphe de voisinage sur T = {1, 2, · · · , T} associé aux potentiels de π étant celui aux 2-plus proches voisins (2-ppv), π est aussi un champ de Markov bilatéral aux 2-ppv ([Kindermannet Snell, 1980];[Guyon, 1995];[Lauritzen, 1996]) : si 1 < t < T ,π(zt | zs, 1 ≤ s ≤ T et s 6= t) = π(zt | zt−1, zt+1) (1)1En eﬀet :π(zt | zs, 1 ≤ s ≤ T et s 6= t) = π(z1, z2, · · · , zT )Pu∈E π(z1, z2, · · · , zt−1, u, zt+1, ..zT )=exp{θt(zt) +Ψt(zt−1, zt) +Ψt+1(zt, zt+1)}Pu∈E exp{θt(u) +Ψt(zt−1, u) +Ψt+1(u, zt+1)}= π(zt | zt−1, zt+1).La loi conditionnelle bilatérale π(zt | zt−1, zt+1) s’explicite donc facilement dès que N , le cardinalde E, n’est pas trop grand.2 Propriété de chaîne de MarkovProposition 1 Z est une chaîne de Markov : π(zt | zs, s ≤ t− 1) = π(zt | zt−1) si 1 < t ≤ T .Preuve : π(zt | zs, s ≤ t− 1) = πt(z1,z2,··· ,zt)πt−1(z1,z2,··· ,zt−1) . Identifions πt(z1, z2, · · · , zt). Pour 1 ≤ s < t ≤T , notant uts = (us, us+1, · · · , ut), on a :πt(z1, z2, · · · , zt) =XuTt+1∈ET−tπ(z1, z2, · · · , zt, uTt+1)= C exp{Xs=1,tθs(zs) +Xs=2,tΨs(zs−1, zs)} × exp{θ∗t (zt)} oùexp{θ∗t (zt)} =XuTt+1exp{TXs=t+1θs(us) +Ψt+1(zt, ut+1) +TXt+2Ψs(us−1, us)}.On en déduit que la loi de zt conditionnelle au passé vaut :π(zt | zs, s ≤ t− 1) = exp{θt(zt) + θ∗t (zt)− θ∗t−1(zt−1) +Ψt(zt−1, zt)} = π(zt | zt−1). (2)¥Remarques :1. Dans la preuve, on peut encore écrire la loi marginale πt(z1, z2, · · · , zt) sous la formeπt(z(t)) = C exp{[Xs=1,t−1θs(zs)] + [θt(zt) + θ∗t (zt)] +Xs=2,tΨs(zs−1, zs)}.Ainsi, le champ marginal (Z1, · · · , Zt) est encore un champ de Gibbs aux 2-ppv pourles mêmes potentiels à l’exception du potentiel de singleton terminal qui vaut eθt(zt) =θt(zt) + θ∗t (zt).2. Une diﬀérence importante apparaît entre les formules (2) et (1) : (1) est calculable alorsque (2) ne l’est pas du fait de la sommation en uTt+1 de complexité NT−t.3 Récursions sur les lois marginales3.1 Contribution conditionnelle au futur Γt(z(t))Pour t ≤ T − 1, la loi conditionnelle au futur, π(z1, z2, · · · , zt | zt+1, zt+2, · · · , zT ), ne dépendque de zt+1. En eﬀet,π(z1, z2, · · · , zt | zt+1, zt+2, · · · , zT ) = π(z1, z2, · · · , zT )Put1∈Et π(ut1, zt+1, ..zT )=exp {Ut(z1, ..., zt) +Ψt+1(zt, zt+1)}Put1∈Et expnPs=1,t θs(us) +Ψt+1(ut, zt+1) +Ps=2,tΨs(us−1, us)o= π(z1, z2, · · · , zt | zt+1).2On peut encore écrire cette loi conditionnelle sous la forme suivante :π(z1, z2, · · · , zt | zt+1) = exp {Ut(z1, ..., zt) +Ψt+1(zt, zt+1)}Put1∈Et exp {Ut(u1, ..., ut) +Ψt+1(ut, zt+1)}, c-à-dπ(z1, z2, · · · , zt | zt+1) = Ct(zt+1) expU∗t (z1, z2, · · · , zt; zt+1)où U∗t est l’énergie conditionnelle au futur définie par :U∗t (z1, z2, · · · , zt; zt+1) = Ut(z1, z2, · · · , zt) +Ψt+1(zt, zt+1). (3)et Ct+1(zt+1)−1 =Put1∈Et exp {U∗t (u1, ..., ut; zt+1)}. Donc pour i = 1, N :π(z1, z2, · · · , zt | zt+1 = ei) = Ct(ei)γt(z1, z2, · · · , zt; ei) où γt(z(t); ei) = expU∗t (z(t); ei).Ainsi, γt(z(t); ei) est la contribution à la loi πt de Z(t) conditionnelle au futur zt+1 = ei.Definition 1 Pour t ≤ T − 1, le vecteur Γt(z(t)) des contributions conditionnelles aufutur est le vecteur de RN de i-ième coordonnée, 1 ≤ i ≤ N :(Γt(z(t)))i = γt(z(t); ei)Pour t = T , il n’y a pas de conditionnement et ΓT (z(T )) est le vecteur constant de coor-données γT (z(T )) = expUT (z(T )). La définition de ΓT (z(T )) est analogue à celle de Γt pourt ≤ T − 1 sous la convention que ΨT+1 ≡ 0. Avec cette même convention, définissons pour1 ≤ t ≤ T la matrice At de taille N ×N de terme général :At(i, j) = exp{θt(ej) +Ψt+1(ej , ei)}, et i, j = 1, N. (4)Notons que AT est une matrice à colonne constante égale à (exp{θT (ej)}, j = 1, N). On a alorsla récurrence fondamentale suivante :Proposition 2 Pour tout 1 ≤ t ≤ T , z(t) = (z1, z2, · · · , zt) ∈ Et et ei ∈ E, on a :γt(z(t− 1), ej ; ei) = At(i, j)× γt−1(z(t− 1); ej). (5)et Xzt∈EΓt(z(t− 1), zt) = AtΓt−1(z(t− 1)) (6)Preuve : (i) Examinons d’abord le cas où 2 ≤ t ≤ T − 1. L’énergie Ut vérifiant Ut(z(t− 1), zt) =Ut−1(z(t− 1)) + θt(zt) +Ψt(zt−1, zt), (3) donne pour tout (zt, zt+1) = (a, b) ∈ E2 :U∗t (z(t− 1), a; b) = Ut−1(z(t− 1)) + θt(a) +Ψt(zt−1, a) +Ψt+1(a, b)= U∗t−1(z(t− 1); a) + {θt(a) +Ψt+1(a, b)}.La récurrence (5) en découle en posant (ej , ei) = (a, b) = (zt, zt+1). La sommation en zt = ejdonne alors la coordonnée i du terme de gauche de (6). D’où le résultat.(ii) Pour t = T , puisque ΨT+1 ≡ 0, U∗T (z(T − 1), a) = U∗T−1(z(T − 1); a)+ θT (a). On obtientalors (6) pour AT (i, j) = exp θT (ej). ¥33.2 Calcul récursif des lois marginales et de la constante de normalisationCommençons par donner quelques notations sur les vecteurs lignes 1 × N suivants : soientE1 = BT = (1, 0, · · · , 0) et la suite (Bt)t=T définie par la récursion Bt−1 = BtAt si t ≤ T ; onnote également K1 =Pz1∈E Γ1(z1) ∈ RN . Le calcul de K1 est facile si N n’est pas trop grandpuisque sa i-ème coordonnée vaut K1i =Pz1∈E exp{θ1(z1) +Ψ2(z1, ei)}.Proposition 3 Calcul des lois marginales πt et de la constante de normalisation C.(1) Pour 1 ≤ t ≤ T :πt(z(t)) = C ×BtΓt(z(t)). (7)(2) La constante de normalisation C de la loi globale π vérifie :C−1 = E1ATAT−1 · · ·A2K1. (8)Preuve :(1) Vérifions (7) par récurrence. Pour t = T , la relation est bien vérifiée puisque,π(z1, z2, · · · , zT ) = πT (z(T )) = C expUT (z(T )) = C ×E1ΓT (z(T )).Supposons (7) vérifiée à l’ordre t, 2 ≤ t ≤ T . Utilisant (6) :πt−1(z(t− 1)) =Xzt∈Eπt(z(t− 1), zt) = C ×Bt{Xzt∈EΓt(z(t− 1), zt)}= C ×BtAtΓt−1(z(t− 1)) = C ×Bt−1Γt−1(z(t− 1)).(2) D’après (7), π1(z1) = C × B1 Γ1(z1). Sommant en z1, on obtient C × B1K1 = 1. Lerésultat résulte de l’égalité B1 = E1ATAT−1 · · ·A2. ¥La formule (8) permettant le calcul de C se simplifie si les potentiels de π sont invariants dansle temps. En eﬀet, dans ce cas, At ≡ A pour 1 ≤ t ≤ T −1, avec A(i, j) = exp{θ(ej)+Ψ(ej , ei)},AT (i, j) = exp θ(ej) et C−1 = E1ATAT−2K1. Si la taille N de E permet la diagonalisation deA, le calcul de C est possible indépendemment de la dimension temporelle T .******************************Exemples1. J’ai fait la vérification de la formule donnant C−1 pour trois situations, les plus simples :(i) un processus à , ceci pour T = 2 et T = 3. Ca marche.(ii) pour E = {0, 1}2 (N = 4 états), les vecteurs E1,K1 de dimension 4 et les matricesA sont 4 × 4 avec T = 2 et le modèle, posant z = (x, y), et le modèle type Ising nonisotropique :π2((x1, y1), (x2, y2)) = C exp{αx1 + βy1 + γx1y1 + αx2 + βy2 + γx2y2 + δ(x1x2 + y1y2)}.C’est mieux de ne pas mettre d’isotropie, on controle mieux les calculs avec les 4 paramètres.Ici il suﬃt de calculer A2, en fait sa première ligne A2(1, ·) = E1A2, et K1. La for-mule n’est autre que C−1 = hA2(1, ·),K1i, soit 16 termes. Et on retrouve bien C−1 =Pz1,z2 expU(z1, z2) !2. Après quoi, il faudra s’assurer de la forme des potentiels θt et Ψt associés au modèle d’Isinget voir le calcul eﬀectif. Dans ce calcul eﬀectif, par exemple sur une image 10×10, il n’y apas lieu de diagonaliser. Il est plus facile de faire le produit de matrices A10×A9×· · ·×A2avec une matrice de “bord” A10 et les autres égales à A2. K1 ne pose pas de problème.43. La diagonalisation vaudra la peine si T est grand, et ce serait bien de valoriser un telexemple. Une chose envisageable, qui nécessite du travail numérique, est de faire del’estimation par MV pour une observation Ising sur un rectangle S = {1,M}×{1, T} d’unmodèle d’Ising. On simule par échantillonneur de Gibbs un tel Ising et on examine lesrésultats asymptotiques pour T grand de comparaison du PMV, Codage et MV. Là, onpeut explicitement programmer la maximisation de la vraie vraisemblance.******************************J’ai suivi ton exple 1 maisje reprendrai l’exple du ising classique avec {-1,+1}La programmation (en matlab) donne un résultat immédiat pour les calculs de C- sous la forme explicite en fonction des valeurs propres λ1 et λ2, voir ci-dessous Meth1- sous la forme C−1 = E1ATAT−2K1 Meth3- sous la forme C−1 = E1ATPDT−2P−1K1 Meth2 et Meth2bis ( selon qu’onutilise les matrices de passage explicites ci-dessous (Meth2) ou celles automatiquementcalculées par matlab (Meth2bis))et un calcul un peu plus long sous la forme de sommation. Meth4Voici le tabelau des temps de caclul, en secondesα = 1, β = 0.7 Meth1 Meth2 Meth2b Meth3 Meth4 Valeur C−1T = 10 0 0 0 0 0.0150 3.6420e+ 007T = 20 0 0 0 0 41.22 2.2611e+ 015Pour un T plus grand, j’ai des pbs de mémoire avec matlab, il n’est pas content d’avoir unematrice de taille 2 puissance T ! il va falloir que je fasse autrement pour stocker.1. *****************************************************Exemple 1Prenons le modèle à N = 2 états avec E = {0, 1} = {e1, e2}, de potentiels de singletonsθt(zt) = αzt et de paires Ψt+1(zt, zt+1) = βztzt+1 pour t ≤ T − 1.Dans ce cas, on obtient At = A =µ1 eα1 eα+β¶pour t = 1, T − 1, AT =µ1 eα1 eα¶, E1 =¡1 0¢et K1 =µ1 + eα1 + eα+β¶.Diagonalisation de A;det(A− λI) = λ2 − λ(1 + eα+β) + eα+β − eα∆ = 1 + e2(α+β) − 2eα+β + 4eα = (1− eα+β)2 + 4eα > 0Deux valeurs propres λ1 = 1+eα+β−√∆2 et λ2 =1+eα+β+√∆2Vecteurs propres associés u1 =µeαλ1 − 1¶et u2 =µeαλ2 − 1¶d’où P =µeα eαλ1 − 1 λ2 − 1¶et P−1 = 1√∆e−αµλ2 − 1 −eα1− λ1 eα¶C−1 = E1ATPDT−2P−1K1C−1 = 1√∆h(λT−11 (λ2 − 1) + λT−12 (1− λ1))(1 + eα) + eα(1 + eα+β)(λT−12 − λT−11 )i***********************RéférencesKindermann R. et Snell J.L. (1980), Markov random fields and their applications, Contemp.Maths.Guyon X. (1995), Random fields on a network : modeling, statistics and applications,Springer5Lauritzen, Graphical modelling (1996), Clarendon Press Oxford.Lovinson G. (2006), A matrix-valued Bernoulli distribution, Journal of Multivariate Analysis,97:1573-1585.Khaled M. (2008), A multivariate generalization of a Markov switching model, working paper,C.E.S., Université Paris 1.Khaled M. (2008), Estimation bayésienne de modèles espace-état non linéaires, Thèse del’Université Paris1.6

Chaine et champ de markov : quelques récurrences et calcul d'une constante de normalisation

https://hal.archives-ouvertes.fr/hal-00397201

Chaine et champ de markov : quelques récurrences et calcul d'une constante de normalisation

Abstract

Similar works

Full text

Available Versions

HAL-Paris1