Replizieren eines Experiments mit GMRF (Gaußsches Markov-Zufallsfeld)

Ich versuche, ein Experiment aus zu verstehen Papier , insbesondere Abschnitt 5.2.

In diesem Artikel schlagen sie einen neuen Algorithmus zur Berechnung der logarithmischen Determinante von dünn besetzten Matrizen vor und testen ihn in Abschnitt 5 an einem von ihnen generierten Datensatz.

Sie möchten es an einem synthetischen Datensatz testen und erstellen daher eine dünne Matrix der Größe 5000×5000, deren Präzisionsmatrix (die Umkehrung der Kovarianzmatrix) durch parametrisiert wird $ \ rho = -0,22 $ . Dem Papier zufolge hat jeder Knoten 4 Nachbarn mit einer partiellen Korrelation von $ \ rho $ . Dann verwenden sie einen Gibbs-Sampler , um eine Probe aus der multivariaten Gaußschen Verteilung zu entnehmen, die durch die Matrix J beschrieben wird Beispiel: Ich berechne die Log-Wahrscheinlichkeit wie folgt: $$ \ log (x | \ rho) = \ log \ det J (\ rho) – x ^ TJ (\ rho) x + G $$ . und ich zeichne die Werte wie in Abbildung 2 auf.

Wenn mein Verständnis korrekt ist, bewerten sie die Log-Wahrscheinlichkeit bei nur einer Stichprobe? Ich verstehe, dass die Darstellung in Abbildung 2 die folgende Formel oben ist, die nur für eine Stichprobe berechnet wird. Normalerweise berechne ich die Log-Wahrscheinlichkeit für einen Datensatz, nicht nur für eine einzelne Stichprobe.

Meine Frage lautet wie folgt: Was genau bedeutet $ \ rho $ und wie erstelle ich $ J (\ rho) $ und Beispiel davon? (dh mit einem Python-Paket? Andernfalls, was ist der Algorithmus?)?

Ich denke, die zugrunde liegende Annahme ist, dass der $ \ log \ det J (\ rho ) $ für zwei verschiedene Beispiele von $ J (\ rho) $ ist das gleiche, warum?

Ich habe tatsächlich nachgesehen zu dem vielzitierten referenzierten Buch , das ein sehr gutes Buch über GMRF ist, aber ich habe keine klare Verbindung zwischen einem einzelnen Parameter $ \ rho $ und die von ihnen erzeugte Matrix. Die Parametrisierung von GMRFs wird in Abschnitt 2.7, Seite 87 beschrieben. Dort wird niemals ein einzelner Parameter verwendet, und der Parameterraum wird tatsächlich durch einen zweidimensionalen Vektor $ \ Theta $ :

$$ \ pi (x | \ Theta) \ propto exp (\ frac {- \ theta_1} {2} \ sum_ {i \ ca. j} (x_i – x_j) ^ 2 – \ frac {\ theta_2} {2} \ sum_i x_i ^ 2) $$ Aber wahrscheinlich beziehen sie sich auf eine andere Matrix.

Update Eigentlich denke ich, dass die Präzisionsmatrix $ J (\ rho) $ , das die Interaktion zwischen 4 Nachbarn beschreibt, ist nur eine -Bandmatrix , dh eine Matrix mit mehreren Diagonalen. In diesem Fall (ich stelle mir vor) mit 2 oberen und 2 unteren Diagonalen, alle gefüllt mit $ – 0,22 $ und nur $ 1 $ auf der Hauptdiagonale.

Wie kann ich dennoch eine Stichprobe aus der durch die Präzisionsmatrix beschriebenen Verteilung ziehen? Sollte ich es invertieren und die Kovarianzmatrix der Daten erhalten und dann eine Stichprobe daraus erstellen? Wenn ja, finden Sie unten den entsprechenden Code. Es kann für jemanden nützlich sein, den Code zu sehen, den wir zum Abtasten aus diesem GMRF verwenden können, vorausgesetzt, der Mittelwert $ \ vec (0) $ und eine Matrixdimension von 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)