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

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

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

f:id:StatModeling:20201106165153p:plain

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

1次元の1階階差のIGMRF

線状につながれたGMRFの場合、 $\textbf{Q}$ は以下になります。

f:id:StatModeling:20201106165157p:plain

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

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

f:id:StatModeling:20201106165201p:plain

1次元の2階階差のIGMRF

線状につながれたGMRFの場合、 $\textbf{Q}$ は以下になります。

f:id:StatModeling:20201106165204p:plain

ここで、 $\textbf{Q} \textbf{S}_{1} = \boldsymbol{0}$ という制約を満たします。 $\textbf{S}_{1}$ は以下の行列です。

f:id:StatModeling:20201106165209p:plain

これは $\boldsymbol{x}$ に直線を足しても差の差（例： $(x_3 - x_2) - (x_2 - x_1)$ ）しか尤度に関わってこないので、直線の傾きがキャンセルされて全体の尤度は不変であることがこの制約の由来です。数式で表現すると、 $\boldsymbol{x}$ を $\boldsymbol{x} + \textbf{S}_1 \boldsymbol{\beta}$ に変えても、 $( \boldsymbol{x} + \textbf{S}_1 \boldsymbol{\beta})^T \textbf{Q} ( \boldsymbol{x} + \textbf{S}_1 \boldsymbol{\beta} )$ は $\boldsymbol{x}^T \textbf{Q} \boldsymbol{x}$ と同じになります。ここで $\boldsymbol{\beta} = ( \beta_0 \ \ \beta_1 )^T$ です。 $i$ を位置のインデックスとして $x_i$ を $x_i + \beta_0 + \beta_1 i$ に変えても同じという意味です。 $\textbf{Q} \textbf{S}_{1} = \boldsymbol{0}$ という制約のため、 $\textbf{Q}$ のrankは2だけ減少して $n-2$ になります。よって、IGMRFの対数尤度は以下になります。

f:id:StatModeling:20201106165212p:plain

$d$ 次元の $k$ 階階差のIGMRF

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

f:id:StatModeling:20201106165216p:plain

ここで、 $n$ はnodeの数です。そこから二項係数を引いています。ここで、 $\boldsymbol{x}$ をすべてのnodeを1列に並べたベクトルとします。すると、 $\boldsymbol{x}$ を $\boldsymbol{x} + \textbf{S} \boldsymbol{\beta}$ に変えても尤度は不変となります。この $\textbf{S}$ の列の数（＝ $\boldsymbol{\beta}$ の要素数）がrankの減少分に相当します。これまでの例では $\textbf{S}$ は、 $\boldsymbol{1}$ だったり、 $\textbf{S}_{1}$ だったりしました。これまでの例と同じように、 $(x_{i_1}, x_{i_2}, ..., x_{i_n})$ に加えても不変になる項を求めると以下になります。

f:id:StatModeling:20201106165223p:plain

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

f:id:StatModeling:20201106165223p:plain

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

　 $\beta_{00} + \beta_{10} i + \beta_{01} j + \frac{1}{2} \beta_{20} i^2 + \beta_{11} i j + \frac{1}{2} \beta_{02} j^2$

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

f:id:StatModeling:20201106165227p:plain

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);
}

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

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

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.