Replikering af et eksperiment på GMRF (Gaussian Markov Random Field)

Jeg prøver at forstå et eksperiment fra dette papir , specifikt afsnit 5.2.

I papiret foreslår de en ny algoritme til beregning af log-determinanten af sparsomme matricer, og i afsnit 5 tester de den på et datasæt, som de genererer.

De vil teste det på et syntetisk datasæt, så de opretter en sparsom matrix i størrelsen 5000×5000, hvis præcisionsmatrix (den omvendte af kovariansmatrixen) er parametreret af $ \ rho = -0.22 $ . Ifølge papiret har hver knude 4 naboer med delvis korrelation på $ \ rho $ . Derefter bruger de en Gibbs sampler til at tage en prøve fra den multivariate gaussiske fordeling, der er beskrevet af matrixen J. På denne prøve, jeg beregner log-sandsynligheden som: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . og jeg tegner værdierne som i figur 2.

Hvis min forståelse er korrekt, evaluerer de log-sandsynligheden kun givet en prøve? Jeg forstår plottet i figur 2 er følgende formel ovenfor, som kun beregnes for en prøve? Jeg beregner normalt log-sandsynligheden på et datasæt, ikke kun på en enkelt prøve.

Mit spørgsmål er følgende: hvad er præcis betydningen $ \ rho $ , og hvordan opretter jeg $ J (\ rho) $ og prøve fra det? (dvs. med en pythonpakke? ellers, hvad er algoritmen?)?

Jeg tror, at den underliggende antagelse er, at $ \ log \ det J (\ rho ) $ for to forskellige prøver af $ J (\ rho) $ er den samme, hvorfor?

Jeg gik faktisk for at se til den meget citerede refererede bog , som er en meget god bog om GMRF, men jeg har ikke fundet nogen klar sammenhæng mellem en enkelt parameter $ \ rho $ og den matrix, de genererer. Parameteriseringen af GMRFer er beskrevet i afsnit 2.7, side 87. Der bruges aldrig en enkelt parameter, og parameterrummet er faktisk beskrevet af en 2-dimensionel 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 sandsynligvis henviser de til en anden matrix.

Opdater Faktisk tror jeg, at præcisionsmatricen $ J (\ rho) $ som beskriver interaktionen mellem 4 naboer er bare en båndmatrix dvs. en matrix med flere diagonaler. I dette tilfælde (forestiller jeg mig) med 2 øvre og 2 nederste diagonaler, alle fyldt med $ – 0.22 $ , og bare $ 1 $ på hoveddiagonalen.

Stadig, hvordan kan jeg prøve fra den distribution, der er beskrevet i præcisionsmatrixen? Skal jeg invertere det og hente dataens kovariansmatrix og derefter prøve det? Hvis ja, nedenfor er koden til at gøre det. Det kan være nyttigt for nogen at se koden, vi kan bruge til at prøve fra denne GMRF, forudsat at $ \ vec (0) $ betyder og en matrixdimension 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)