Répliquer une expérience sur GMRF (champ aléatoire de Markov gaussien)

Jessaie de comprendre une expérience de ceci papier , en particulier la section 5.2.

Dans larticle, ils proposent un nouvel algorithme pour calculer le log-déterminant des matrices creuses, et dans la section 5 ils le testent sur un ensemble de données quils génèrent.

Ils veulent le tester sur un jeu de données synthétique, ils créent donc une matrice clairsemée de taille 5000×5000 dont la matrice de précision (linverse de la matrice de covariance) est paramétrée par $ \ rho = -0,22 $ . Selon larticle, chaque nœud a 4 voisins avec une corrélation partielle de $ \ rho $ . Ensuite, ils utilisent un échantillonneur Gibbs pour prélever un échantillon de la distribution gaussienne multivariée qui est décrite par la matrice J. Sur ce exemple, je calcule le log-vraisemblance comme suit: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . et je trace les valeurs comme dans la figure 2.

Si ma compréhension est correcte, ils évaluent la log-vraisemblance pour un seul échantillon? Je comprends que le graphique de la figure 2 est la formule ci-dessus, qui est calculée uniquement pour un échantillon? Je calcule généralement le log-vraisemblance sur un ensemble de données, pas seulement sur un seul échantillon.

Ma question est la suivante: quelle est exactement la signification $ \ rho $ , et comment créer $ J (\ rho) $ et échantillonner ? (cest-à-dire avec un package python? Sinon, quel est lalgorithme?)?

Je pense que lhypothèse sous-jacente est que le $ \ log \ det J (\ rho ) $ pour deux échantillons différents de $ J (\ rho) $ est le même, pourquoi?

En fait, je suis allé chercher au livre référencé , qui est un très bon livre sur GMRF, mais je nai trouvé aucun lien clair entre un seul paramètre $ \ rho $ et la matrice quils génèrent. Le paramétrage des GMRF est décrit dans la section 2.7, page 87. Là, un seul paramètre nest jamais utilisé, et lespace des paramètres est en fait décrit par un vecteur à 2 dimensions $ \ Theta $ :

$$ \ pi (x | \ Theta) \ propto exp (\ frac {- \ theta_1} {2} \ sum_ {i \ approx j} (x_i – x_j) ^ 2 – \ frac {\ theta_2} {2} \ sum_i x_i ^ 2) $$ Mais ils font probablement référence à une matrice différente.

Mise à jour En fait, je pense que la matrice de précision $ J (\ rho) $ qui décrivent linteraction entre 4 voisins est juste une matrice de bande cest-à-dire une matrice à diagonales multiples. Dans ce cas (jimagine) avec 2 diagonales supérieures et 2 inférieures, toutes remplies avec $ – 0.22 $ , et juste $ 1 $ sur la diagonale principale.

Pourtant, comment puis-je échantillonner à partir de la distribution décrite par la matrice de précision? Dois-je linverser et obtenir la matrice de covariance des données, puis échantillonner à partir de celle-ci? Si oui, vous trouverez ci-dessous le code pour le faire. Il pourrait être utile pour quelquun de voir le code que nous pouvons utiliser pour échantillonner à partir de ce GMRF, en supposant que $ \ vec (0) $ signifie et une dimension de matrice de 30.

import numpy as np def precision(k, numero): return np.diag(np.repeat(1, k)) + np.diag(np.repeat(numero, k-1), -1) + np.diag(np.repeat(numero, k-2), -2) + np.diag(np.repeat(numero, k-1), 1) + np.diag(np.repeat(numero, k-2), 2) J = precision(30, -0.22) Sigma = np.linalg.inv(J) np.random.multivariate_normal(np.repeat(0, 30), Sigma)