適応棄却メトロポリス・サンプリング

概要

・パッケージdlmの中にARMSという関数があり、これを使うことで適応棄却メトロポリス・サンプリングが使えるらしいのでメモがでらデモ
・とりあえず動いているように見えるが詳しいロジックまでは追えていないので悪しからず

適応棄却メトロポリス・サンプリング

任意の分布\( \pi \)からのサンプリングを行う単純なアルゴリズムらしい
『正確に目標分布からの独立な抽出系列を生成することができるのに対し、メトロポリスーヘイスティング・アルゴリズムでは、シミュレートされた確率変数には一般に従属性があり、その極限においてのみ目標分布に従う分布を持つ』

系列相関を持たないサンプリングが可能ですと。

ARMSを使ってみる(1変量)

まずはお試しで1変量の簡単な分布のサンプリングをやってみる。対象とする目標分布は適当に正規分布の組み合わせ

140524_1

ARMSを使って上記分布からサンプリングするには以下。

関数fには、先ほどの関数gの対数を返すようにします。
下から2行目でarmsを使ってサンプリングをしておりまして、
第一引数:初期値
第二引数:対象分布(対数)
第三引数:サンプリング範囲?
第四引数:サンプリング個数
というような感じかと。

実際にサンプルのヒストグラムを見ると

140524_2

という感じに、うまくサンプリングされているように見えます。
また、サンプルの系列相関を見てみると、以下のようになり、系列相関はほぼ無いと判断できそうですね。

140524_3

ARMSを使ってみる(2変量)

実はこのARMSは多変量にも対応しているようでして、例えば2変量正規分布の場合

関数fは普通の2変量正規分布の密度関数の定義です。これはcontourで等高線図を書くためにとりあえず定義しています。
ARMS用の対数密度は関数gで定義しています。

140524_4

なんとなーく分布には従っていそうな感じ

ARMSを使ってみる(2変量・多峰性)

普通にギブスサンプリングを行うと多峰性の分布を対象とした場合、1個の山から抜け出すことができなくて正確にサンプリングすることが困難。一方でARMSを使えば大丈夫(らしい)。例えば以下のうような分布の場合。

140524_5

参考文献

Amazon.co.jp: Rによるベイジアン動的線形モデル (統計ライブラリー): G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎: 本
http://www.amazon.co.jp/dp/4254127960

ウィシャート分布についてメモ

概要

・タイトルどおりウィシャート分布がよくイメージ湧かないので簡単な俺様メモ
・簡単な期待値の計算など

分散共分散行列(or 精度行列)をモデル化するときに出てくる

多変量のベイズ推定なんかをしようとするとよくお目にかかる分布。よくわからん。
単変量の場合は精度(分散の逆数)をガンマ分布に当てはめると都合が良かったですが多変量の場合はウィシャート分布がそれに該当すると。

行列に対する分布なのでイマイチイメージが湧かない。
(行列の各要素に対する同時分布なんですが)

で、例えば\( y_i \)が\( m \)変量の多変量正規分布\( N_m \left( 0, \Sigma \right) \)に従うとすると、標本分散共分散行列\( W \)は
$$
W = \sum_{i=1}^{N} y_i \cdot y_i^T
$$
となり、これがウィシャート分布に従いまして
$$
W \sim W_m \left( \frac{N}{2}, \frac{1}{2} \Sigma^{-1} \right)
$$
と書かれます。
で、一方ウィシャート分布\( W_k \left( \alpha, B \right) \)に従う行列\( A \)の期待値は
$$
E \left( A \right) = \alpha B^{-1}
$$
と書けますので、先ほどの標本分散共分散行列\( W \)の期待値は
$$
E \left( W \right) = \frac{N}{2} \cdot \left( \frac{1}{2} \Sigma^{-1} \right)^{-1} = n \Sigma
$$
となります。

実際に計算してみる using MASS

多変量正規分布に従う乱数はMASSパッケージで取得します。
とりあえずということで分散共分散行列は適当においています。

まず周辺分布から確認すると

となりましてまあ平均0, 分散1に従っていそうというのは確認できます。
ちなみにプロットすると以下の様な感じ

140517_1

あとは何回かこの作業を繰り返しぶん回すことで分散共分散行列の平均を計算する

で得られたのがこちら

実際に計算してみる using dlm

実はdlmにはウィシャート分布に従う乱数を取得するrwishartという関数がありましてこちらでも確認することができます。

出てきた値がこちら。

実際の期待値

計算するまでもないですが期待値は \( n \Sigma \)ですので

で計算され

ということでまあ多分あってそうという感じでしょうか。

参考文献

Amazon.co.jp: Rによるベイジアン動的線形モデル (統計ライブラリー): G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎: 本
http://www.amazon.co.jp/dp/4254127960

Amazon.co.jp: パターン認識と機械学習 上: C.M. ビショップ, 元田 浩, 栗田 多喜夫, 樋口 知之, 松本 裕治, 村田 昇: 本
http://www.amazon.co.jp/dp/4621061224

dlmを用いて状態推定、信頼区間の表示

概要

・点推定ではなく信頼区間も計算せよとの天の声を頂いたのでdlmを使って計算してみる

使用データ

いつもNFPだと飽きるので今回は新規失業保険申請件数(Initial Claim)を対象とします。

Initial Claims – FRED – St. Louis Fed
http://research.stlouisfed.org/fred2/series/ICSA

プロットしてみるとこのような感じ(FREDのサイトで確認できますが。)

150510_1

失業保険の申請件数(週次)なもんですから、景気が安定しているときは減少トレンドを保ちつつ相対的に低い水準で推移する傾向があるようですね。ただし景気が悪くなるとかなり激しく上昇します。ボラティリティと似たような動きをしていますな

カルマンフィルタ

今回は線形成長モデルを用います。

まず最尤推定を行うことで観測方程式の分散、状態方程式の分散を推定。
推定結果は以下

信頼区間の計算

信頼区間を計算するためには、各時点での状態の分布(というか標準偏差)が欲しいわけです。
この分布の分散はdlmの中では\( C_t \)と定義されています。
ただし、dlmの中では\( C_t \)がそのまま保存されているわけではなくて特異値分解(SVD)後のU.CとD.Cという2つの要素に別れて保存されています。なのでその2つの要素から\( C_t \)を構築してあげる必要があります。

ただそのためのメソッドはdlm内にちゃんと準備してありましてdlmSvd2var(そのままなネーミング)を呼ぶことで計算できます。

普通にdlmSvd2varを実行しただけだと行列のリストが返ってきます。一方で状態(ローカルレベル)の標準偏差として本当に欲しいのは[1,1]要素だけなので5行目のように書くことで[1,1]要素だけを抽出しています。
あとは\( C_0 \)を取り除いたり後々のプロットのためにxts型にしたりとしています。

図示

データ数が結構あるので見た目的にも2013年以降を表示。
今回はカッコつけて系列名なども表示しています。

150510_2

2013年以降でみると概ね95%の信頼区間内に入っていますがそれでも時々上方にブレているのがありますな。

ちなみにリーマン・ショック期の2008~2009年を見るとこちら
140510_3

参考文献

Amazon.co.jp: Rによるベイジアン動的線形モデル (統計ライブラリー): G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎: 本
http://www.amazon.co.jp/dp/4254127960

自作関数での感度分析

概要

・mapplyの練習も兼ねて自作関数の感度分析という初心者ポスト^^
・未だにapply族の関数が使いこなせない
・for文で感度分析は書きたくない

Rコード

適当にBlack-Scholesを評価してみます。感度分析するパラメータは行使価格Kとします。

で、行使価格以外のパラメータは任意の値で固定しつつ、行使価格を80~100で移動させながら評価させていきます。
それをやっているのが下記3行目のmapply

150506_1

for文で書くよりはスッキリ書けますね^^