Replicând un experiment pe GMRF (Gaussian Markov Random Field)

Încerc să înțeleg un experiment din hârtie , în special secțiunea 5.2.

În lucrare, ei propun un nou algoritm pentru calcularea log-determinantului matricelor rare, iar în secțiunea 5 îl testează pe un set de date pe care îl generează.

Vor să-l testeze pe un set de date sintetic, astfel încât să creeze o matrice rară de dimensiunea 5000×5000 a cărei matrice de precizie (inversa matricei de covarianță) este parametrizată de $ \ rho = -0.22 $ . Conform lucrării, fiecare nod are 4 vecini cu corelație parțială de $ \ rho $ . Apoi, utilizează un eșantionator Gibbs pentru a prelua un eșantion din distribuția gaussiană multivariată care este descrisă de matricea J. eșantion, calculez probabilitatea jurnalului ca: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . și trasez valorile ca în Figura 2.

Dacă înțelegerea mea este corectă, ei evaluează probabilitatea jurnalului având doar un eșantion? Înțeleg că graficul din figura 2 este următoarea formulă de mai sus, care se calculează numai pentru eșantionul unul ? De obicei calculez log-probabilitatea pe un set de date, nu doar pe un singur eșantion.

Întrebarea mea este următoarea: care este exact sensul $ \ rho $ , și cum pot crea $ J (\ rho) $ și eșantion din acesta? (adică cu un pachet python? altfel, care este algoritmul?)?

Cred că presupunerea de bază este că $ \ log \ det J (\ rho ) $ pentru două eșantioane diferite de $ J (\ rho) $ este același, de ce?

De fapt, m-am dus să caut la mult citată carte de referință , care este o carte foarte bună despre GMRF, dar nu am găsit nicio legătură clară între un singur parametru $ \ rho $ și matricea pe care o generează. Parametrizarea GMRF-urilor este descrisă în secțiunea 2.7, pagina 87. Acolo, un singur parametru nu este folosit niciodată, iar spațiul parametrilor este de fapt descris de un vector 2 dimensiuni $ \ 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) $$ Dar probabil se referă la o altă matrice.

Actualizați De fapt, cred că matricea de precizie $ J (\ rho) $ care descrie interacțiunea dintre 4 vecini este doar o matrice de bandă , adică o matrice cu diagonale multiple. În acest caz (îmi imaginez) cu 2 diagonale superioare și 2 diagonale inferioare, toate umplute cu $ – 0.22 $ și doar $ 1 $ pe diagonala principală.

Totuși, cum pot preleva din distribuția descrisă de matricea de precizie? Ar trebui să-l inversez și să obțin matricea de covarianță a datelor și apoi să prelevez din ea? Dacă da, mai jos este codul pentru ao face. Ar putea fi util ca cineva să vadă codul pe care îl putem folosi pentru a preleva din acest GMRF, presupunând $ \ vec (0) $ medie și o dimensiune a matricei 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)