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

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

今日読んだページ： 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. 交換可能性を課してしまう。これを課すとド・フィネッティの表現定理が成り立ち、「事前分布がパラメトリックであれば各がしたがう分布はの下でである」までいえる。
  - ちなみに、9ページのド・フィネッティの表現定理では「事前分布がパラメトリック」は課していないので $\pi$ は経験分布関数がとりうる関数そのものの分布である（下図のようなイメージ）。
- Ex. 上記のような仮定は強すぎるので各 $Y_t$ が対応する $\theta_t$ の下で独立というようにゆるめる（この場合ベクトル $\theta_{1:n}$ がしたがう確率分布を求めることがゴールになる）。
ではの分布まで求まったら推定値をどう考えればよいのか：
- 考える損失関数によって事後密度／予測密度を統合することになる。
- 観測値が $\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")