Replikerar ett experiment på GMRF (Gaussian Markov Random Field)

Jag försöker förstå ett experiment från detta papper , särskilt avsnitt 5.2.

I tidningen föreslår de en ny algoritm för beräkning av logg-determinanten för glesa matriser, och i avsnitt 5 testar de den i en dataset som de genererar.

De vill testa det på en syntetisk dataset, så de skapar en gles matris av storlek 5000×5000 vars precisionsmatris (den inversa av kovariansmatrisen) parametreras av $ \ rho = -0.22 $ . Enligt tidningen har varje nod fyra grannar med partiell korrelation på $ \ rho $ . Sedan använder de en Gibbs sampler för att ta ett prov från den multivariata gaussiska fördelningen som beskrivs av matrisen J. På denna prov, jag beräknar logg-sannolikheten som: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . och jag plottar värdena som i figur 2.

Om min förståelse stämmer, utvärderar de logg sannolikheten med endast ett prov? Jag förstår att handlingen i figur 2 är följande formel ovan, som bara beräknas för ett prov? Jag brukar beräkna log-sannolikheten för en dataset, inte bara för ett enda exempel.

Min fråga är följande: vad är meningen $ \ rho $ , och hur skapar jag $ J (\ rho) $ och exempel från det? (dvs. med ett pythonpaket? annars, vad är algoritmen?)?

Jag tror att det underliggande antagandet är att $ \ log \ det J (\ rho ) $ för två olika exempel på $ J (\ rho) $ är samma, varför?

Jag gick faktiskt för att titta till den mycket citerade refererade boken , vilket är en mycket bra bok om GMRF, men jag har inte hittat någon tydlig koppling mellan en enda parameter $ \ rho $ och matrisen de genererar. Parametreringen av GMRF: er beskrivs i avsnitt 2.7, sida 87. Där används aldrig en enda parameter och parameterutrymmet beskrivs faktiskt av en tvådimensionell vektor $ \ 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) $$ Men förmodligen hänvisar de till en annan matris.

Uppdatera Egentligen tror jag precisionsmatrisen $ J (\ rho) $ som beskriver interaktionen mellan fyra grannar är bara en bandmatris dvs en matris med flera diagonaler. I det här fallet (föreställer jag mig) med två övre och två nedre diagonaler, alla fyllda med $ – 0.22 $ , och bara $ 1 $ på huvuddiagonalen.

Hur kan jag ändå prova från den distribution som beskrivs av precisionsmatrisen? Ska jag invertera den och få kovariansmatrisen för datan och sedan prova den? Om ja, nedan är koden för att göra det. Det kan vara användbart för någon att se koden vi kan använda för att sampla från denna GMRF, förutsatt att $ \ vec (0) $ betyder och en matrisdimension på 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)