IGMRFの尤度におけるrankの減少分に関するメモ

以下の書籍を読んで、IGMRF（Intrinsic Gaussian Markov Ramdom Field）の尤度に関して自分の理解をまとめたメモです。

[asin:B00YBV6YLI:detail]

この発表資料の18ページにおいて、（観測モデル部分を除いた）IGMRFの対数尤度は以下に比例すると書きました（ただしをに、をに読み替えてください）。

ここではnodeの数、はの精度行列*1でnodeのつながりの情報を反映していて、は線形制約に由来するのrankの減少分です。

1次元の1階階差のIGMRF

線状につながれたGMRFの場合、は以下になります。

このように精度行列が帯行列になってスパースになるところがGMRFの特徴です（分散共分散行列はスパースにならない）。

ここで、という線形の制約を満たすことに注意してください。に定数を足して平行移動しても、の要素の差しか尤度に関わってこないので、全体の尤度は不変であることがこの制約の由来です。数式で表現すると、をに平行移動しても、はと同じになります。を位置のインデックスとしてをに変えても同じという意味です。という制約のため、のrankは1だけ減少して、となります。よって、IGMRFの対数尤度は以下になります。

1次元の2階階差のIGMRF

線状につながれたGMRFの場合、は以下になります。

ここで、という制約を満たします。 は以下の行列です。

これはに直線を足しても差の差（例：）しか尤度に関わってこないので、直線の傾きがキャンセルされて全体の尤度は不変であることがこの制約の由来です。数式で表現すると、をに変えても、はと同じになります。ここでです。を位置のインデックスとしてをに変えても同じという意味です。という制約のため、のrankは2だけ減少してになります。よって、IGMRFの対数尤度は以下になります。

次元の階階差のIGMRF

次元の階階差のGMRFではrankはいくつになるでしょうか。結論を先に書くと、以下になります。

ここで、はnodeの数です。そこから二項係数を引いています。ここで、をすべてのnodeを1列に並べたベクトルとします。すると、をに変えても尤度は不変となります。このの列の数（＝の要素数）がrankの減少分に相当します。これまでの例ではは、だったり、だったりしました。これまでの例と同じように、 に加えても不変になる項を求めると以下になります。

ここから、の要素数は次のように数えます。の添字の数が個で、各添字の数値の範囲が～で、添字の数値の合計は以下となる制限があります。この状況は、部屋の仕切り個と個の玉の並び替えの場合の数と同じになるので、の要素数は以下になります。

例えば2次元3階階差の場合、 に加えても不変になる項は、以下になります。

　

の要素数は6個になりますので、rankの減少分は6になります。例えば、は、以下の部屋の仕切りと玉の図に相当します。

Stanによる実装例

「StanとRでベイズ統計モデリング」の12.8節の例題（2次元1階階差の地図）のStanコードは以下になります。

data {
  int N;
  vector[N] Y;
  int I;
  int<lower=1, upper=N> From[I];
  int<lower=1, upper=N> To[I];
}

parameters {
  vector[N] r;
  real<lower=0> s_r;
  real<lower=0> s_Y;
}

model {
  target += - (N-1) * log(s_r) - 0.5 * dot_self(r[To] - r[From]) * inv_square(s_r);
  Y ~ normal(r, s_Y);
}

なお、先にを頑張って構築してquad_formなどを利用してモデルを書く方法は、上記のように書く方法より遅くなりました。ベクトル化がうまく効かないからだと思います。なお、12.8節のモデルと比べるとの部分は等価ですが、外側のlog(s_r)に掛かる定数が異なります（12.8節のモデルがよりが小さくなる方へペナルティが入っています）。

*1:厳密にはrankが落ちているので精度行列そのものとは言えませんが、ここでは[4]にならって区別しないでそう呼びます。