Replicando un experimento en GMRF (Gaussian Markov Random Field)

Estoy tratando de entender un experimento de este paper , específicamente la Sección 5.2.

En el artículo, proponen un nuevo algoritmo para calcular el determinante logarítmico de matrices dispersas, y en la sección 5 lo prueban en un conjunto de datos que generan.

Quieren probarlo en un conjunto de datos sintéticos, por lo que crean una matriz dispersa de tamaño 5000×5000 cuya matriz de precisión (la inversa de la matriz de covarianza) está parametrizada por $ \ rho = -0,22 $ . Según el documento, cada nodo tiene 4 vecinos con una correlación parcial de $ \ rho $ . Luego, usan un muestreador de Gibbs para tomar una muestra de la distribución gaussiana multivariante que se describe en la matriz J. En este muestra, calculo la probabilidad logarítmica como: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . y trazo los valores como en la Figura 2.

Si mi comprensión es correcta, ¿evalúan la probabilidad logarítmica dada solo una muestra? Entiendo que el gráfico de la figura 2 es la siguiente fórmula anterior, que se calcula solo para una muestra. Por lo general, calculo la probabilidad de registro en un conjunto de datos, no solo en una sola muestra.

Mi pregunta es la siguiente: cuál es exactamente el significado $ \ rho $ y cómo creo $ J (\ rho) $ y muestra de ella? (es decir, ¿con un paquete de Python? De lo contrario, ¿cuál es el algoritmo?)?

Creo que la suposición subyacente es que el $ \ log \ det J (\ rho ) $ para dos muestras diferentes de $ J (\ rho) $ es lo mismo, ¿por qué?

De hecho, fui a buscar al muy citado libro referenciado , que es un muy buen libro sobre GMRF, pero no he encontrado ningún vínculo claro entre un solo parámetro $ \ rho $ y la matriz que generan. La parametrización de GMRFs se describe en la Sección 2.7, página 87. Allí, un solo parámetro nunca se usa, y el espacio de parámetros se describe en realidad mediante un vector bidimensional $ \ 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) $$ Pero probablemente se estén refiriendo a una matriz diferente.

Actualización En realidad, creo que la matriz de precisión $ J (\ rho) $ que describe la interacción entre 4 vecinos es solo una matriz de bandas , es decir, una matriz con múltiples diagonales. En este caso (imagino) con 2 diagonales superiores y 2 inferiores, todas llenas de $ – 0.22 $ , y solo $ 1 $ en la diagonal principal.

Aún así, ¿cómo puedo tomar muestras de la distribución descrita por la matriz de precisión? ¿Debo invertirlo y obtener la matriz de covarianza de los datos y luego tomar una muestra de ella? Si es así, a continuación se muestra el código para hacerlo. Puede ser útil que alguien vea el código que podemos usar para muestrear de este GMRF, asumiendo que $ \ vec (0) $ significa y una dimensión de matriz de 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)