Replicando um experimento em GMRF (Gaussian Markov Random Field)

Estou tentando entender um experimento deste papel , especificamente a Seção 5.2.

No artigo, eles propõem um novo algoritmo para calcular o log-determinante de matrizes esparsas e, na seção 5, eles o testam em um conjunto de dados que geram.

Eles querem testá-lo em um conjunto de dados sintético, então eles criam uma matriz esparsa de tamanho 5000×5000 cuja matriz de precisão (o inverso da matriz de covariância) é parametrizada por $ \ rho = -0,22 $ . De acordo com o artigo, cada nó possui 4 vizinhos com correlação parcial de $ \ rho $ . Em seguida, eles usam um amostrador de Gibbs para obter uma amostra da distribuição gaussiana multivariada que é descrita pela matriz J. amostra, eu calculo a probabilidade de log como: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . e ploto os valores como na Figura 2.

Se meu entendimento estiver correto, eles avaliam o log da probabilidade dada apenas uma amostra? Eu entendo que o gráfico na figura 2 é a seguinte fórmula acima, que é calculada apenas para uma amostra? Eu geralmente calculo a probabilidade de log em um conjunto de dados, não apenas em uma única amostra.

Minha pergunta é a seguinte: qual é exatamente o significado $ \ rho $ , e como faço para criar $ J (\ rho) $ e amostra dele? (ou seja, com um pacote python? caso contrário, qual é o algoritmo?)?

Acho que a suposição subjacente é que o $ \ log \ det J (\ rho ) $ para duas amostras diferentes de $ J (\ rho) $ é o mesmo, por quê?

Na verdade, fui olhar ao muito citado livro referenciado , que é um livro muito bom sobre GMRF, mas não encontrei nenhum link claro entre um único parâmetro $ \ rho $ e a matriz que eles geram. A parametrização de GMRFs é descrita na Seção 2.7, página 87. Lá, um único parâmetro nunca é usado, e o espaço de parâmetro é realmente descrito por um vetor bidimensional $ \ Theta $ :

$$ \ pi (x | \ Theta) \ propto exp (\ frac {- \ theta_1} {2} \ sum_ {i \ aprox j} (x_i – x_j) ^ 2 – \ frac {\ theta_2} {2} \ sum_i x_i ^ 2) $$ Mas provavelmente eles estão se referindo a uma matriz diferente.

Atualizar Na verdade, acho que a matriz de precisão $ J (\ rho) $ que descreve a interação entre 4 vizinhos é apenas uma matriz de banda ou seja, uma matriz com múltiplas diagonais. Neste caso (eu imagino) com 2 diagonais superiores e 2 inferiores, todas preenchidas com $ – 0,22 $ e apenas $ 1 $ na diagonal principal.

Ainda assim, como posso obter uma amostra da distribuição descrita pela matriz de precisão? Devo inverter e obter a matriz de covariância dos dados e, em seguida, fazer uma amostragem a partir dela? Se sim, abaixo está o código para fazer isso. Pode ser útil para alguém ver o código que podemos usar para obter uma amostra deste GMRF, assumindo $ \ vec (0) $ média e uma dimensão 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)