Een experiment repliceren op GMRF (Gaussian Markov Random Field)

Ik probeer een experiment te begrijpen van dit paper , in het bijzonder Sectie 5.2.

In de paper stellen ze een nieuw algoritme voor voor het berekenen van de log-determinant van sparse matrices, en in sectie 5 testen ze het op een dataset die ze genereren.

Ze willen het testen op een synthetische dataset, dus maken ze een dunne matrix van 5000×5000, waarvan de precisiematrix (de inverse van de covariantiematrix) wordt geparametriseerd door $ \ rho = -0,22 $ . Volgens de paper heeft elk knooppunt 4 buren met een gedeeltelijke correlatie van $ \ rho $ . Vervolgens gebruiken ze een Gibbs-sampler om één steekproef te nemen uit de multivariate Gauss-verdeling die wordt beschreven door de matrix J. Hierop voorbeeld, ik bereken de log-waarschijnlijkheid als: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . en ik zet de waarden uit zoals in figuur 2.

Als ik het goed begrijp, evalueren ze de log-waarschijnlijkheid gegeven slechts één steekproef? Ik begrijp dat de plot in figuur 2 de volgende formule hierboven is, die alleen wordt berekend voor één steekproef? Meestal bereken ik de log-waarschijnlijkheid op een dataset, niet alleen op een enkele steekproef.

Mijn vraag is de volgende: wat is precies de betekenis $ \ rho $ , en hoe maak ik $ J (\ rho) $ en voorbeeld ervan? (dwz met een python-pakket? Wat is anders het algoritme?)?

Ik denk dat de onderliggende aanname is dat de $ \ log \ det J (\ rho ) $ voor twee verschillende voorbeelden van $ J (\ rho) $ is hetzelfde, waarom?

Ik ging echt kijken naar het veel geciteerde boek waarnaar wordt verwezen , dat een zeer goed boek is over GMRF, maar ik heb geen duidelijke link gevonden tussen een enkele parameter $ \ rho $ en de matrix die ze genereren. De parametrering van GMRFs wordt beschreven in Paragraaf 2.7, pagina 87. Daar wordt nooit een enkele parameter gebruikt en wordt de parameterruimte feitelijk beschreven door een tweedimensionale vector $ \ 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) $$ Maar waarschijnlijk verwijzen ze naar een andere matrix.

Update Eigenlijk denk ik dat de precisiematrix $ J (\ rho) $ die de interactie tussen 4 buren beschrijven, is slechts een bandmatrix , dwz een matrix met meerdere diagonalen. In dit geval (stel ik me voor) met 2 bovenste en 2 onderste diagonalen, allemaal gevuld met $ – 0,22 $ , en slechts $ 1 $ op de hoofddiagonaal.

Maar hoe kan ik bemonsteren uit de distributie die wordt beschreven door de precisiematrix? Moet ik het omkeren en de covariantiematrix van de gegevens verkrijgen en er vervolgens een steekproef van nemen? Zo ja, hieronder staat de code om het te doen. Het kan voor iemand nuttig zijn om de code te zien die we kunnen gebruiken om een steekproef uit deze GMRF te nemen, uitgaande van een $ \ vec (0) $ gemiddelde en een matrixdimensie van 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)