読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

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

久保緑本11章のマルコフ場モデル(空間構造のあるベイズモデル)

Stan

空間ベイズモデルの一種であるマルコフ場モデルがやっと分かった気がしますので備忘録を書きます。 きっかけとなったのは以下の動画の復習です。0:13:50~0:16:30と1:14:50~1:32:20の部分です。統数研の伊庭先生がマルコフ場モデルについての講義をしている箇所です。

今回動画を見て、Stanで実装する上で頑張るところは「まず同時分布を考える。次に同時分布をDAGであらわす。」ところだと理解しました。

1次元の場合

上記を踏まえて久保本11章をDAGで表したStanコード(1次のマルコフ場モデル)は以下になります。

data {
   int<lower=1> N_site;
   int<lower=0> Y[N_site];
}

parameters {
   real r[N_site];
   real<lower=0> s_r;
}

model {
   for (j in 2:N_site)
      r[j] ~ normal(r[j-1], s_r);
   for (j in 1:N_site)
      Y[j] ~ poisson_log(r[j]);
}

generated quantities {
   real<lower=0> Y_mean[N_site];
   for (j in 1:N_site)
      Y_mean[j] <- exp(r[j]);
}
  • 12~13行目: 同時分布をDAGであらわす方法は1つではありませんが、1次元の場合はこう書くよりほかはないでしょう。動画で解説のあった通り、時系列の状態空間モデルと全く同じ形になります。
  • 15行目: Y[j] ~ poisson(exp(r[j])); と書くよりも計算が安定して収束しやすくなります。StanのTipsです。

2次のマルコフ場モデルのStanコードは以下のようになります。

data {
   int<lower=1> N_site;
   int<lower=0> Y[N_site];
}

parameters {
   real r[N_site];
   real<lower=0> s_r;
}

model {
   for (j in 3:N_site)
      r[j] ~ normal(2*r[j-1] - r[j-2], s_r);
   for (j in 1:N_site)
      Y[j] ~ poisson_log(r[j]);
}

generated quantities {
   real<lower=0> Y_mean[N_site];
   for (j in 1:N_site)
      Y_mean[j] <- exp(r[j]);
}
  • 12~13行目: 変化があるのはここだけです。

結果の図は以下のようになります。左図が1次・右図が2次。黒の実線はMCMCサンプルの中央値、グレイの帯は80%信用区間、黒の点線はポアソン分布の真の平均パラメータです。2次の方が信用区間が狭くなる傾向がありました。

20150127_kubo11a0-2.png 20150127_kubo11a1-2.png

2次元の場合

まず格子が3×3の合計9個の場合を考えます。 1次のマルコフ場モデルの場合は隣同士の値が似ていることを考慮した同時分布になります。2次のマルコフ場モデルの場合は傾きの値が似ていること(まっすぐになって欲しいこと)を考慮した同時分布になります。以下、模式図。左図が1次・右図が2次。

20150127_9x9.png

式は順に以下のようになります(s_r=1の場合)。 20150127_image001.png

20150127_image002.png

今回は例として以下のデモデータを作成し用いました。30×30の格子点でZの値を持つデータです。

N <- 30
x <- seq(-10, 10, length=N)
y <- x
f <- function(x, y) { r <- sqrt(x^2+y^2); 10*sin(r)/r }
z <- outer(x, y, f)
set.seed(123)
noise <- matrix(rnorm(N*N, 0, 0.5), nrow=N, ncol=N)
Z <- z + noise

png('persp-Z_data.png', pointsize=24, width=800, height=600)
persp(x, y, Z, theta=30, phi=30, expand=0.5, col=rainbow(N), border=NA)
dev.off()

png('persp-z_hidden.png', pointsize=24, width=800, height=600)
persp(x, y, z, theta=30, phi=30, expand=0.5, col=rainbow(N), border=NA)
dev.off()

グラフを描くと以下のようなデータです。ノイズを入れる前のデータはzで左図、ノイズを入れた推定に使うデータはZで右図です。

20150127_persp-z_hidden-2.png 20150127_persp-Z_data-2.png

Stanでの実装に移ります。同時分布をDAGで表す方法はいくつかありますが、今回は動画でもあったように角(1,1の点)から決めていくようなDAGにします。1次のマルコフ場モデルの場合のStanコードは以下になります。

data {
   int<lower=1> N;
   real Z[N,N];
}

parameters {
   real r[N,N];
   real<lower=0> s_r;
   real<lower=0> s_z;
}

model {
   for (i in 1:N)
      for (j in 2:N)
         r[i,j] ~ normal(r[i,j-1], s_r);
   for (j in 1:N)
      for (i in 2:N)
         r[i,j] ~ normal(r[i-1,j], s_r);
   for (i in 1:N)
      for (j in 1:N)
         Z[i,j] ~ normal(r[i,j], s_z);
}

一見r[2,2]などが「~」の左辺に2回登場するのでオヤッと思うかもしれませんが、r[2,2] ~ normal(r[1,2], s_r)lp__(log posterior)に -log(s_r)-0.5*pow((r[2,2]-r[1,2])/s_r,2) を加えていることと等価であることを考えると、これで正しく同時分布を定義できていることが分かります。

2次のマルコフ場モデルの場合のStanコードは以下になります。

data {
   int<lower=1> N;
   real Z[N,N];
}

parameters {
   real r[N,N];
   real<lower=0> s_r;
   real<lower=0> s_z;
}

model {
   for (i in 1:N)
      for (j in 3:N)
         r[i,j] ~ normal(2*r[i,j-1] - r[i,j-2], s_r);
   for (j in 1:N)
      for (i in 3:N)
         r[i,j] ~ normal(2*r[i-1,j] - r[i-2,j], s_r);
   for (i in 1:N)
      for (j in 1:N)
         Z[i,j] ~ normal(r[i,j], s_z);
}

1次のマルコフ場モデルの場合の、キックするRコードは以下になります。

library(rstan)
...略...
data <- list(N=N, Z=Z)
stan_fit <- stan(file="model1.stan", chains=0)

fit <- stan(
   fit=stan_fit,
   data=data,
   init=function() { list(r=Z, s_r=1, s_z=1) },
   iter=1200, warmup=200, thin=1,
   seed=1234, chains=3
)

la <- extract(fit)
r_median <- apply(la$r, c(2,3), median)
png('persp-M1-r-median.png', pointsize=24, width=800, height=600)
persp(x, y, r_median, theta=30, phi=30, expand=0.5, col=rainbow(N), border=NA)
dev.off()

大切なのは9行目でrの初期値をかなりちゃんと設定しないとうまく収束しませんでした。なお計算時間はSurface Pro 3(core i5)で1次のマルコフ場モデルが1chainあたり10秒ぐらい、2次のマルコフ場モデルが1chainあたり30秒ぐらいでした。

結果はパラメータrのMCMCサンプルの中央値だけとってきて、3次元でplotしたものが以下になります。左図が1次・右図が2次。この図からは分かりませんが、こちらも2次の方が信用区間が狭くなる傾向がありました。

20150127_persp-M1-r-median-2.png 20150127_persp-M2-r-median-2.png

これ以上はやりませんが、地図の場合でもインデックスをうまく持てば問題なくいけそうですね。

追記

「DAGで表している」という表現は誤解を招くのではないかという指摘がありましたので補足します。たしかにr[2,2]などは「~」の左側に2回登場しますしますので厳密なDAGではないですね…。むしろStanで「~」を使った書き方がたまたま無向グラフの同時分布を適切に表現できているという感じなので、誤解を招きやすい「~」を使わないでincrement_log_prob()を使って2次元の2次マルコフ場モデルを書きなおしたStanコードを以下に掲載しておきます。

data {
   int<lower=1> N;
   real Z[N,N];
}

parameters {
   real r[N,N];
   real<lower=0> s_r;
   real<lower=0> s_z;
}

model {
   for (i in 1:N)
      for (j in 3:N)
         increment_log_prob(normal_log(r[i,j], 2*r[i,j-1] - r[i,j-2], s_r));
   for (j in 1:N)
      for (i in 3:N)
         increment_log_prob(normal_log(r[i,j], 2*r[i-1,j] - r[i-2,j], s_r));
   for (i in 1:N)
      for (j in 1:N)
         Z[i,j] ~ normal(r[i,j], s_z);
}
  • 2016.12.25追記

CARモデルとマルコフ場モデルは、同時分布のexpの中身は同じなのですが、正規化項(大きすぎる標準偏差へのペナルティ項)がわずかに異なります。その項は対数尤度では-N*log(sigma)となりますが、CARモデルはガウス過程と同じでNはノード数になりますが、マルコフ場モデルではNはエッジ数になります。そこで、CARモデルと書かないでマルコフ場モデルという表記にしました。