Replicando un esperimento su GMRF (Gaussian Markov Random Field)

Sto cercando di capire un esperimento da questo paper , in particolare la Sezione 5.2.

Nel documento, propongono un nuovo algoritmo per calcolare il determinante logaritmico di matrici sparse e nella sezione 5 lo testano su un set di dati che generano.

Vogliono testarlo su un set di dati sintetico, quindi creano una matrice sparsa di dimensioni 5000 x 5000 la cui matrice di precisione (linverso della matrice di covarianza) è parametrizzata da $ \ rho = -0,22 $ . Secondo larticolo, ogni nodo ha 4 vicini con correlazione parziale di $ \ rho $ . Quindi, utilizzano un campionatore Gibbs per prelevare un campione dalla distribuzione gaussiana multivariata descritta dalla matrice J. Su questo sample, calcolo la verosimiglianza come: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . e tracciamo i valori come nella Figura 2.

Se la mia comprensione è corretta, valutano la probabilità logaritmica dato un solo campione? Capisco che il grafico nella figura 2 sia la seguente formula sopra, che viene calcolata solo per un campione? Di solito calcolo la probabilità di log su un set di dati, non solo su un singolo campione.

La mia domanda è la seguente: qual è esattamente il significato $ \ rho $ e come faccio a creare $ J (\ rho) $ e sample da esso? (cioè con un pacchetto Python? altrimenti, qual è lalgoritmo?)?

Penso che il presupposto sottostante sia che $ \ log \ det J (\ rho ) $ per due diversi esempi di $ J (\ rho) $ è lo stesso, perché?

In realtà sono andato a guardare al tanto citato libro di riferimento , che è un ottimo libro su GMRF, ma non ho trovato alcun collegamento chiaro tra un singolo parametro $ \ rho $ e la matrice che generano. La parametrizzazione dei GMRF è descritta nella Sezione 2.7, pagina 87. Lì, un singolo parametro non viene mai utilizzato e lo spazio dei parametri è effettivamente descritto da un vettore bidimensionale $ \ 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) $$ Ma probabilmente si riferiscono a una matrice diversa.

Aggiorna In realtà, penso che la matrice di precisione $ J (\ rho) $ che descrive linterazione tra 4 vicini è solo una matrice a bande , cioè una matrice con diagonali multiple. In questo caso (immagino) con 2 diagonali superiori e 2 inferiori, tutte riempite con $ – 0,22 $ e solo $ 1 $ sulla diagonale principale.

Tuttavia, come posso campionare dalla distribuzione descritta dalla matrice di precisione? Devo invertirlo e ottenere la matrice di covarianza dei dati e quindi campionarla da essa? Se sì, di seguito è riportato il codice per farlo. Potrebbe essere utile per qualcuno vedere il codice che possiamo usare per campionare da questo GMRF, supponendo che $ \ vec (0) $ media e una dimensione di matrice di 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)