Replikowanie eksperymentu na GMRF (Gaussian Markov Random Field)

Próbuję zrozumieć eksperyment z tego papier , w szczególności sekcja 5.2.

W artykule zaproponowali nowy algorytm obliczania wyznacznika logarytmicznego rzadkich macierzy, aw sekcji 5 testują go na generowanym przez siebie zbiorze danych.

Chcą to przetestować na syntetycznym zbiorze danych, więc tworzą rzadką macierz o rozmiarze 5000×5000, której macierz dokładności (odwrotność macierzy kowariancji) jest sparametryzowana przez $ \ rho = -0,22 $ . Według artykułu każdy węzeł ma 4 sąsiadów z częściową korelacją $ \ rho $ . Następnie używają samplera Gibbsa , aby pobrać jedną próbkę z wielowymiarowego rozkładu gaussowskiego opisanego przez macierz J. przykład, prawdopodobieństwo dziennika obliczam jako: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . i wykreślam wartości tak, jak na rysunku 2.

Jeśli moje rozumienie jest prawidłowe, oceniają one logarytm prawdopodobieństwa dla jednej próbki? Rozumiem, że wykres na rysunku 2 to następujący wzór powyżej, który jest obliczany tylko dla jednej próbki? Zwykle obliczam prawdopodobieństwo dziennika na zbiorze danych, a nie tylko na pojedynczej próbce.

Moje pytanie jest następujące: jakie dokładnie jest znaczenie $ \ rho $ i jak utworzyć $ J (\ rho) $ i próbka z tego? (tj. z pakietem Pythona? w przeciwnym razie jaki jest algorytm?)?

Myślę, że podstawowym założeniem jest to, że $ \ log \ det J (\ rho ) $ za dwie różne próbki $ J (\ rho) $ to to samo, dlaczego?

Właściwie poszedłem sprawdzić do często cytowanej książki , która jest bardzo dobrą książką o GMRF, ale nie znalazłem żadnego wyraźnego powiązania między jednym parametrem $ \ rho $ i wygenerowaną przez nie macierz. Parametryzacja GMRF jest opisana w Sekcji 2.7, strona 87. W tym przypadku pojedynczy parametr nigdy nie jest używany, a przestrzeń parametrów jest w rzeczywistości opisana przez dwuwymiarowy wektor $ \ Theta $ :

$$ \ pi (x | \ Theta) \ propto exp (\ frac {- \ theta_1} {2} \ sum_ {i \ ok j} (x_i – x_j) ^ 2 – \ frac {\ theta_2} {2} \ sum_i x_i ^ 2) $$ Ale prawdopodobnie odnoszą się do innej macierzy.

Zaktualizuj Właściwie myślę, że macierz precyzji $ J (\ rho) $ , które opisują interakcję między 4 sąsiadami, to po prostu macierz pasmowa , czyli macierz z wieloma przekątnymi. W tym przypadku (wyobrażam sobie) z 2 górnymi i 2 dolnymi przekątnymi, wszystkie wypełnione $ – 0,22 $ i tylko 1 $ $ na głównej przekątnej.

Mimo to, jak mogę pobrać próbkę z rozkładu opisanego przez macierz dokładności? Czy powinienem go odwrócić i uzyskać macierz kowariancji danych, a następnie pobrać z niej próbkę? Jeśli tak, poniżej znajduje się kod, aby to zrobić. Ktoś mógłby zobaczyć kod, którego możemy użyć do próbkowania z tego GMRF, zakładając $ \ vec (0) $ średnią i wymiar macierzy 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)