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

Rによるベイジアン動的線型モデル: ノート4

本読み DLM

読んでいる本(出典): Rによるベイジアン動的線形モデル (統計ライブラリー) | G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎 | 本 | Amazon.co.jp

前回:ノート3 / 次回:ノート5
目次:Rによるベイジアン動的線型モデル

今日読んだページ: 22~35ページ
以下、自分の解釈。誤っている可能性があります。

  • 前回までのあらすじ:
    • 線形回帰モデルに対しても共役な事前分布を選択し、回帰係数/分散共分散行列/その両方のベイズ推定が可能。
  • マルコフ連鎖モンテカルロ法の導入:
    • しばしばパラメータ  \psi の事後分布は陽に扱えない。→ そこでモンテカルロ法
    • 事後分布からの標本はマルコフ連鎖によりシミュレーションするのが標準的。
    • 一般的なマルコフ連鎖モンテカルロ法:
    • 適応棄却メトロポリス・サンプリング(ARMS):
      • 棄却サンプリング: 目標分布に似た提案分布を使用し、受容確率により採否を判定。
      • 適応棄却サンプリング: 棄却が発生したときに提案分布を目標分布に近づけるよう更新。
        下図のようなイメージ。但し、(Cに)無駄があるような絵。f:id:cookie-box:20160101171656p:plain:w420



R ではパッケージ dlm に ARMS が含まれているということなのでつかってみる。
元旦にちなんだ確率分布をつくる。→ 2000点サンプリング → できた(ARMSの無駄遣いが)。
f:id:cookie-box:20160101183451p:plain:w450
コード(どうでもいい)は以下。ARMS が対数密度を使用する都合上、"2016" の文字の外も0でない密度がある(on と off の密度差を contrast で定義 = 100.0)。

my.dens <- function(x) {
  contrast <- 100.0
  off <- 1.0
  on <- contrast
  z <- on * 40 + off * 65
  off <- off / z
  on <- on / z
  if (x[[1]] < 1 || 15 <= x[[1]] || x[[2]] < 1 || 6 <= x[[2]]) {
    return(off)
  }
  if (5 <= x[[1]] && x[[1]] < 6) { return(on) }
  if (7 <= x[[1]] && x[[1]] < 8) { return(on) }
  if (9 <= x[[1]] && x[[1]] < 10) { return(on) }
  if (11 <= x[[1]] && x[[1]] < 12) { return(on) }
  if (x[[2]] < 2 || 5 <= x[[2]]) {
    if (1 <= x[[1]] && x[[1]] < 4) { return(on) }
    if (6 <= x[[1]] && x[[1]] < 7) { return(on) }
    if (12 <= x[[1]] && x[[1]] < 14) { return(on) }
  } else if (x[[2]] < 3) {
    if (1 <= x[[1]] && x[[1]] < 2) { return(on) }
    if (13 <= x[[1]] && x[[1]] < 14) { return(on) }
  } else if (x[[2]] < 4) {
    if (1 <= x[[1]] && x[[1]] < 4) { return(on) }
    if (12 <= x[[1]] && x[[1]] < 14) { return(on) }
  } else {
    if (3 <= x[[1]] && x[[1]] < 4) { return(on) }
  }
  return(off)
}
my.log.dens <- function(x) { return(log(my.dens(x))) }
supp <- function(x) { (0 <= x[[1]] && x[[1]] < 15 && 0 <= x[[2]] && x[[2]] < 7) }
y <- arms(c(1, 1), my.log.dens, supp, 2000)
plot(y, xlim=c(0, 15), ylim=c(0, 7))



演習問題をやりたかったができていない。