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

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

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

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

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

  • やりたいこと(前回の範囲):
    •  Y_{1:n}=y_{1:n} を観測したとき、その観測値を実現させるモデルのパラメータ  \theta を知りたい。或いは、 \theta には特に興味はなく次の観測値  Y_{n+1} そのものを予測したい。
  • やりたいことへのベイズアプローチ(前回の範囲):
    •  \thetaパラメトリックな事前分布と  \theta の下での  Y_{1:n} の尤度関数は予め何か決め打っておいて、取得した観測値  y_{1:n} の下でのパラメータ  \theta の事後分布を求める。
  • といっても観測値がとる分布への制約をどう考えていけばよいのか:
    • Ex. 条件付き独立性( \theta の下で  Y_1, \dots, Y_n {\rm i.i.d.} にしたがう)を課してしまう。
    • Ex. 交換可能性を課してしまう。これを課すとド・フィネッティの表現定理が成り立ち、「事前分布がパラメトリックであれば各  Y_t がしたがう分布は  \theta の下で  {\rm i.i.d.} である」までいえる。
      • ちなみに、9ページのド・フィネッティの表現定理では「事前分布がパラメトリック」は課していないので  \pi は経験分布関数がとりうる関数そのものの分布である(下図のようなイメージ)。f:id:cookie-box:20151230155538p:plain:w230
    • Ex. 上記のような仮定は強すぎるので各  Y_t が対応する  \theta_t の下で独立というようにゆるめる(この場合ベクトル  \theta_{1:n} がしたがう確率分布を求めることがゴールになる)。
  • では  \theta の分布まで求まったら推定値をどう考えればよいのか:
    • 考える損失関数によって事後密度/予測密度を統合することになる。
    • 観測値が  \theta の下で  {\rm i.i.d.} であるとき、 n が大きいときは  \theta の推定値は MLE を中心とする正規密度に近づく(ベイズ主義での推定と頻度主義での推定が合致する例)。
    • しかし、分散共分散行列が既知の  \sigma^2 I である3次元以上の多変量正規分布の平均ベクトルの推定にあたっては、MLE(観測値ベクトルそのものになる)は二次損失に関して最適ではなく、より適切な別の推定量が存在する(観測が互いに独立であるにもかかわらず、 t 番目の観測値の分布の平均の推定値に  j \neq t 番目の観測値がかかわってくる!?)。



二次損失、線形損失、0-1損失を実感してみるために、 \theta の事後分布が以下のような離散分布(適当)になるケースを考える。

pi <- data.frame(
  theta=c(  1,   2,    3,    4,   5,    6,    7),
  p    =c(0.1, 0.25, 0.23, 0.2, 0.1, 0.07, 0.05)
)

このとき、

  •  \theta の平均値: 3.36
  •  \theta の中央値: 3
  •  \theta の最頻値: 2

 \theta の分布と損失関数を図示すると下図のようになる。二次損失は(見た感じ)  \theta の平均値で最小に、線形損失は中央値で最小に、0-1損失は最頻値で最小になっていることがわかる。
f:id:cookie-box:20151230185050p:plain:w480

※ 実際には 0-1 損失は  \theta の取りうる値でのみ  0 でない値を取る。グラフが汚いのはサンプリングの問題。
ソースは以下。

plot(pi$theta, pi$p, xlim=c(0.5, 7.5), xlab="θ", ylab="PDF", type="h")

quadratic.loss <- function(x){sum(pi$p * (pi$theta - x) * (pi$theta - x))}
linear.loss    <- function(x){sum(pi$p * abs(x - pi$theta))}
zeroone.loss   <- function(x){sum(pi$p * ifelse(x == pi$theta, 0, 1))}

theta.target <- seq(0.5, 7.5, by=0.01)
plot(theta.target, sapply(theta.target, quadratic.loss), xlab="a", ylab="Quadratic loss", type="l")
plot(theta.target, sapply(theta.target, linear.loss   ), xlab="a", ylab="Linear loss",    type="l")
plot(theta.target, sapply(theta.target, zeroone.loss  ), xlab="a", ylab="Zero-one loss",  type="l")