GMRF（ガウスマルコフ確率場）での実験の複製

から実験を理解しようとしています。紙、特にセクション5.2。

この論文では、スパース行列の対数行列式を計算するための新しいアルゴリズムを提案し、セクション5では、生成したデータセットでそれをテストします。

合成データセットでテストしたいので、サイズ5000×5000のスパース行列を作成します。この行列の精度行列（共分散行列の逆行列）は、によってパラメーター化されます。 $ \ rho = -0.22 $ 。論文によると、各ノードには $ \ rho $ の偏相関を持つ4つの隣接ノードがあります。次に、ギブスサンプラーを使用して、行列Jで記述される多変量ガウス分布から1つのサンプルを取得します。サンプルでは、対数尤度を次のように計算します。 $$ \ log（x | \ rho）= \ log \ det J（\ rho）-x ^ TJ（\ rho）x + G $$ 。そして、図2のように値をプロットします。

私の理解が正しければ、1つのサンプルのみが与えられた場合の対数尤度を評価しますか？図2のプロットは、上記の次の式であると理解しています。これは、 1つのサンプルに対してのみ計算されますか？私は通常、単一のサンプルだけでなく、データセットで対数尤度を計算します。

私の質問は次のとおりです。 $ \ rho $の意味は正確には何ですか。 、および create $ J（\ rho）の方法$ とサンプルから？ （つまり、Pythonパッケージを使用しますか？それ以外の場合、アルゴリズムは何ですか？）？

基本的な仮定は、 $ \ log \ det J（\ rho ） $ J（\ rho）$ の2つの異なるサンプルの$ は同じですが、なぜですか？

実際に調べに行きましたよく引用されている参照されている本へ。これはGMRFに関する非常に優れた本ですが、単一のパラメーター $ \ rho $ とそれらが生成する行列。 GMRFのパラメーター化については、87ページのセクション2.7で説明しています。そこでは、単一のパラメーターが使用されることはなく、パラメーター空間は実際には2次元ベクトル $ \ 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）$$ しかし、おそらくそれらは別の行列を参照しています。

更新実際、精度行列 $ 4つの隣接行列間の相互作用を表すJ（\ rho）$ は、単なる帯行列、つまり複数の対角線を持つ行列です。この場合（私が想像する）、2つの上部対角線と2つの下部対角線があり、すべて $ -0.22 $ で埋められ、 $ 1だけで埋められています。主対角線上の$ 。

それでも、精度マトリックスによって記述された分布からどのようにサンプリングできますか？それを反転してデータの共分散行列を取得し、それからサンプリングする必要がありますか？はいの場合、以下はそれを行うためのコードです。 $ \ vec（0）$ の平均と行列の次元が30であると仮定して、このGMRFからサンプリングするために使用できるコードを誰かが見ると役立つ場合があります。

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)