ともにゃん的データ分析ブログ

勉強したことの備忘録とかね

【R】ダミー変数を一度に生成する関数

R package 【R】自作関数

探せばあるんだろうけど、データフレームを引数に複数列を一度にダミー変数化する関数を作りました。
よければ使ってください。


使い方：

(0) {dummies} パッケージをインストールする

(1) 引数 data にダミー変数化したいデータフレームを入れる

(2) 回帰等で多重共線性を避けるために1列除きたい場合は is.drop = TRUE とする

convDummies <- function(data, is.drop = FALSE){
  library(dummies)
  
  N <- ncol(data)
  row_names <- names(data)
  
  names_list <- c()
  new_data <- rep(NA, nrow(data))
  for(n in 1:N){
    unique_value <- sort(unique(data[,n]))
    dummied_data <- dummy(data[,n])
    
    if(is.drop == TRUE){
      new_data <- cbind(new_data, dummied_data[,-ncol(dummied_data)])
      names_list <- c(names_list, 
                      paste(row_names[n], unique_value, sep = ".")[-ncol(dummied_data)])
    } else {
      new_data <- cbind(new_data, dummied_data)
      names_list <- c(names_list, paste(row_names[n], unique_value, sep = "."))
    }
  }
  
  new_data <- as.data.frame(new_data)
  names(new_data) <- c("temp", names_list)
  
  return(new_data[,-1])
}

pythonによる粒子フィルタの実装

python ベイズ統計学 時系列解析

2階差分トレンド+季節(週)トレンドを考慮した以下の状態空間モデル(線形ガウス状態空間モデル)

{ \displaystyle \begin{eqnarray*} y_{t}&=&\mu_{t}+s_{t}+w_{t},\qquad w_{t}\sim N(0,\sigma^{2})\\ \mu_{t}&=&2\mu_{t-1}-\mu_{t-2}+v_{t},\qquad v_{t}\sim N(0,\alpha_{\mu}^{2}\sigma^{2})\\ s_{t}&=&-\sum_{i=1}^{6}s_{t-i}+z_{t},\qquad z_{t}\sim N(0,\alpha_{s}^{2}\sigma^{2}) \end{eqnarray*}\\ }

について、粒子フィルターをpythonで実装しました。
ここでy_{t},\ \mu_{t}, s_{t}はそれぞれ時刻tにおける観測値、トレンド(平均)、季節トレンドです。
状態空間モデルの詳しい説明については

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)


が分かりやすいです。
行列とベクトルの表記を使って綺麗に状態空間モデルのモデリングや種々のアルゴリズムの導出を行っていますし、P.89には粒子フィルタの疑似コードが載っています。


ただしこの本に従って実際に実装するとなると、モデルによっては膨大な量の行列とベクトルを実装しなければなりません。
特に拡大状態ベクトルを用いて固定ラグ平滑化を導出していますが、いざ固定ラグ平滑化を実装しようとなると拡大状態ベクトルの実装が面倒だと感じます。また固定ラグ平滑化については実装例が載っていません。


そこで、固定ラグ平滑化まで含めて粒子フィルタをpythonで実装しましたので参考にしていただけたらと思います。

github.com


使用したデータは上記書籍のサポートページ
http://daweb.ism.ac.jp/yosoku/
からダウンロードできます。



固定ラグ平滑化の概要は以下であると解釈しています。予測とフィルタリングは普通の粒子フィルタです。縦に5つ並んでいる丸が粒子で、この図では5つの粒子が描かれています。x_{i,t|t-1}はt-1時点でのデータで条件付けられたt時点でのi番目の粒子ということを表しています。粒子フィルタはシステムモデルによってt-1時点のデータからt時点の値を"予測"し、実際にt時点でのデータが手に入った時点で観測モデルの尤度でもって状態xをフィルタリング(更新)します。

f:id:kefits:20160805172611p:plain


さて、t時点での最新のデータが手に入った時点で例えばt-1やt-2の過去のデータを更新するというのが平滑化という操作で、平滑化の方法はMCMCといった方法で実行することもできますが、粒子フィルターの枠組みでは固定ラグ平滑化という方法が一般に用いられます。その固定ラグ平滑化の概要は上図のように、最新のデータでもって現在の状態xをフィルタリングした後に、そのフィルタリングされた粒子xは前の時点t-1ではどの粒子から派生したものかというのを逆向きに辿ることで平滑化ができる(と私は解釈しています)。

青い四角で囲まれた数字は、前の時点のどの(上から何番目の)粒子から派生したかのインデクスを表しています。上図ではラグ幅が1の固定ラグ平滑化を実行していることになっていますが、プログラムを実装する上では上図のように各時点の粒子がどの粒子から派生したものかのインデックスを記憶しておいて、それをデータが手に入ってフィルタリングをするたびに記憶しておいたインデックスを頼りに過去に辿っていけば平滑化ができるということです。




今回は冒頭で紹介した書籍を参考に、pythonで実際に粒子フィルタを実装をしてみたわけですが、これをStanで実装した例として@berobero11さんのスライドやブログ

www.slideshare.net


statmodeling.hatenablog.com


は大変参考になります。





間違いや質問等ございましたら@kefismまで遠慮なくご連絡ください。