読んでいる本(出典): Rによるベイジアン動的線形モデル (統計ライブラリー) | G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎 | 本 | Amazon.co.jp
前回:ノート3 / 次回:ノート5
目次:Rによるベイジアン動的線型モデル
今日読んだページ: 22~35ページ
以下、自分の解釈。誤っている可能性があります。
- 前回までのあらすじ:
- 線形回帰モデルに対しても共役な事前分布を選択し、回帰係数/分散共分散行列/その両方のベイズ推定が可能。
- マルコフ連鎖モンテカルロ法の導入:
- しばしばパラメータ の事後分布は陽に扱えない。→ そこでモンテカルロ法。
- 事後分布からの標本はマルコフ連鎖によりシミュレーションするのが標準的。
- 一般的なマルコフ連鎖モンテカルロ法:
- 適応棄却メトロポリス・サンプリング(ARMS):
- 棄却サンプリング: 目標分布に似た提案分布を使用し、受容確率により採否を判定。
- 適応棄却サンプリング: 棄却が発生したときに提案分布を目標分布に近づけるよう更新。
下図のようなイメージ。但し、(Cに)無駄があるような絵。
R ではパッケージ dlm に ARMS が含まれているということなのでつかってみる。
元旦にちなんだ確率分布をつくる。→ 2000点サンプリング → できた(ARMSの無駄遣いが)。
コード(どうでもいい)は以下。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))
演習問題をやりたかったができていない。