Replikerer et eksperiment på GMRF (Gaussian Markov Random Field)

Jeg prøver å forstå et eksperiment fra dette papir , spesielt avsnitt 5.2.

I avisen foreslår de en ny algoritme for beregning av logg-determinanten til sparsomme matriser, og i avsnitt 5 tester de den på et datasett de genererer.

De vil teste det på et syntetisk datasett, slik at de lager en sparsom matrise med størrelse 5000×5000 hvis presisjonsmatrise (den omvendte av kovariansmatrisen) blir parametrisert av $ \ rho = -0.22 $ . Ifølge papiret har hver node fire naboer med delvis korrelasjon på $ \ rho $ . Deretter bruker de en Gibbs sampler for å ta en prøve fra den multivariate gaussiske fordelingen som er beskrevet av matrisen J. På denne prøve, jeg beregner log-sannsynligheten som: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . og jeg plotter verdiene som i figur 2.

Hvis min forståelse stemmer, vurderer de logg sannsynligheten gitt bare ett utvalg? Jeg forstår plottet i figur 2 er følgende formel ovenfor, som bare beregnes for en prøve? Jeg beregner vanligvis sannsynligheten for logg på et datasett, ikke bare på et enkelt utvalg.

Spørsmålet mitt er følgende: hva betyr egentlig $ \ rho $ , og hvordan lager jeg $ J (\ rho) $ og eksempler fra det? (dvs. med en python-pakke? ellers, hva er algoritmen?)?

Jeg tror den underliggende antagelsen er at $ \ log \ det J (\ rho ) $ for to forskjellige eksempler på $ J (\ rho) $ er det samme, hvorfor?

Jeg gikk faktisk for å se til den mye siterte refererte boken , som er en veldig god bok om GMRF, men jeg har ikke funnet noen tydelig kobling mellom en enkelt parameter $ \ rho $ og matrisen de genererer. Parameteriseringen av GMRF er beskrevet i Seksjon 2.7, side 87. Der brukes aldri en enkelt parameter, og parameterområdet er faktisk beskrevet av en todimensjonal 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 sannsynligvis refererer de til en annen matrise.

Oppdater Egentlig tror jeg presisjonsmatrisen $ J (\ rho) $ som beskriver samspillet mellom fire naboer er bare en båndmatrise dvs. en matrise med flere diagonaler. I dette tilfellet (forestiller jeg meg) med to øvre og to nedre diagonaler, alle fylt med $ – 0.22 $ , og bare $ 1 $ på hoveddiagonalen.

Likevel, hvordan kan jeg prøve fra fordelingen beskrevet av presisjonsmatrisen? Bør jeg invertere det og få tak i kovariansematrisen til dataene og deretter prøve det? Hvis ja, nedenfor er koden for å gjøre det. Det kan være nyttig for noen å se koden vi kan bruke til å prøve fra denne GMRF, forutsatt at $ \ vec (0) $ betyr og en matrisedimensjon 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)