Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード Twitter

2015-04-15

tidyrでダミー変数を作成する(赤ペン先生希望)

| 19:11 | tidyrでダミー変数を作成する(赤ペン先生希望)を含むブックマーク

とりあえず必要なパッケージを読み込んでおく。

library(tidyr)
library(dplyr)
library(pipeR)

例えば、以下のようなデータを考える。

> df <- data.frame(id=1:10, feature=c(1,NA,2,2,2,NA,NA,3,3,1))
> df
   id feature
1   1       1
2   2      NA
3   3       2
4   4       2
5   5       2
6   6      NA
7   7      NA
8   8       3
9   9       3
10 10       1

このデータのfeature列の値を回帰分析なりにつっこみたいので、feature列の値毎にダミー変数(列)を作る必要がある。model.matrix関数なんかでイケるのかなと思ったけど、どうもそうではないっぽいので自分でやってみようということです。

とりあえず何も考えずにtidyrパッケージのspreadをすれば良さそうだぞと。NAはうざいのでとりあえず0埋めで。

> df %>>% spread(feature, feature, fill=0)
   id 1 2 3 NA
1   1 1 0 0  0
2   2 0 0 0  0
3   3 0 2 0  0
4   4 0 2 0  0
5   5 0 2 0  0
6   6 0 0 0  0
7   7 0 0 0  0
8   8 0 0 3  0
9   9 0 0 3  0
10 10 1 0 0  0

さらに新しく出来たNA列がうざいので、こいつを削除したい。したいんだが・・・

> df %>>% 
+   spread(feature, feature, fill=0) %>>%
+   select(-contains("NA"))
data frame with 0 columns and 10 rows
> df %>>% 
+   spread(feature, feature, fill=0) %>>%
+   select(-NA)
   id 1 2 3 NA
1   1 1 0 0  0
2   2 0 0 0  0
3   3 0 2 0  0
4   4 0 2 0  0
5   5 0 2 0  0
6   6 0 0 0  0
7   7 0 0 0  0
8   8 0 0 3  0
9   9 0 0 3  0
10 10 1 0 0  0

・・・のように、うまくいかなかった。

〜追記〜

教えてもらったのだが、いったんdata.frameをかましてNA列をNA.列にするか、あるいはNAをバッククオート``で囲めばOKとのこと。

> df %>% 
+   spread(feature, feature, fill=0) %>%
+   data.frame %>%
+   select(-NA.)
   id X1 X2 X3
1   1  1  0  0
2   2  0  0  0
3   3  0  2  0
4   4  0  2  0
5   5  0  2  0
6   6  0  0  0
7   7  0  0  0
8   8  0  0  3
9   9  0  0  3
10 10  1  0  0
> 
> df %>% 
+   spread(feature, feature, fill=0) %>%
+   select(-`NA`)
   id 1 2 3
1   1 1 0 0
2   2 0 0 0
3   3 0 2 0
4   4 0 2 0
5   5 0 2 0
6   6 0 0 0
7   7 0 0 0
8   8 0 0 3
9   9 0 0 3
10 10 1 0 0

〜〜

しょうがなく列名を変更してから削除

> df %>>% 
+   spread(feature, feature, fill=0) %>>% 
+   setNames(c("id", paste0("X", colnames(.)[-1]))) %>>%
+   select(-XNA)
   id X1 X2 X3
1   1  1  0  0
2   2  0  0  0
3   3  0  2  0
4   4  0  2  0
5   5  0  2  0
6   6  0  0  0
7   7  0  0  0
8   8  0  0  3
9   9  0  0  3
10 10  1  0  0

それから、各列の値を全部1に持っていきたいので、以下にifelseで書く&その処理はid列には非適用

> df %>>% 
+   spread(feature, feature, fill=0) %>>% 
+   setNames(c("id", paste0("X", colnames(.)[-1]))) %>>%
+   select(-XNA) %>>%
+   mutate_each(funs(ifelse(.>=1, 1, 0)), -id)
   id X1 X2 X3
1   1  1  0  0
2   2  0  0  0
3   3  0  1  0
4   4  0  1  0
5   5  0  1  0
6   6  0  0  0
7   7  0  0  0
8   8  0  0  1
9   9  0  0  1
10 10  1  0  0

・・・で、御望みの形になった。なったはいいけど、もっといいやり方があるはずなので、赤ペン先生希望します。

追記

コメント欄含めアドバイス多数頂戴した、ありがたいことです。

NA非対応だけどこんなかんじ。はじめから0埋めしておけばいいのか。

> model.matrix(~feature-1, df %>% mutate(feature = factor(feature)))
   feature1 feature2 feature3
1         1        0        0
3         0        1        0
4         0        1        0
5         0        1        0
8         0        0        1
9         0        0        1
10        1        0        0

若干遠回りしてる感あるけど、これでもいいのか。

> library(caret)
> feat <- df %>% mutate(feature = factor(feature))
> obj <- dummyVars(~feature, feat)
> predict(obj, feat)
   feature.1 feature.2 feature.3
1          1         0         0
2         NA        NA        NA
3          0         1         0
4          0         1         0
5          0         1         0
6         NA        NA        NA
7         NA        NA        NA
8          0         0         1
9          0         0         1
10         1         0         0

こちらのスライドの中で紹介されているdummiesパッケージを使うのもよさげ。dummy, dummy.data.frame関数でサクッといく。

> library(dummies)
> dummy(df$feature)
      feature1 feature2 feature3 featureNA
 [1,]        1        0        0         0
 [2,]        0        0        0         1
 [3,]        0        1        0         0
 [4,]        0        1        0         0
 [5,]        0        1        0         0
 [6,]        0        0        0         1
 [7,]        0        0        0         1
 [8,]        0        0        1         0
 [9,]        0        0        1         0
[10,]        1        0        0         0
> dummy.data.frame(df %>% mutate(feature=factor(feature)))
   id feature1 feature2 feature3 featureNA
1   1        1        0        0         0
2   2        0        0        0         1
3   3        0        1        0         0
4   4        0        1        0         0
5   5        0        1        0         0
6   6        0        0        0         1
7   7        0        0        0         1
8   8        0        0        1         0
9   9        0        0        1         0
10 10        1        0        0         0

qdapToolsパッケージというのもあるらしい。こいつのmtabulate関数でもイケる。

> qdapTools::mtabulate(df$feature)
   1 2 3
1  1 0 0
2  0 0 0
3  0 1 0
4  0 1 0
5  0 1 0
6  0 0 0
7  0 0 0
8  0 0 1
9  0 0 1
10 1 0 0

tidyrパッケージの前身であるreshape2でも出来る。出来るっちゃ出来るが抱える問題もtidyrとほぼ同じか。

> reshape2::dcast(df, id~feature, value.var = "feature", length)
   id 1 2 3 NA
1   1 1 0 0  0
2   2 0 0 0  1
3   3 0 1 0  0
4   4 0 1 0  0
5   5 0 1 0  0
6   6 0 0 0  1
7   7 0 0 0  1
8   8 0 0 1  0
9   9 0 0 1  0
10 10 1 0 0  0

通りすがり通りすがり 2015/04/15 19:44 tidyrではなくて恐縮ですがqdapTools::mtabulate(df$feature)

通りすがり通りすがり 2015/04/15 21:01 reshape2::dcast(df, id~feature, value.var = "feature", length)

teramonagiteramonagi 2015/04/15 22:58 いずれにせよ、ありがとうございます!

2015-04-11

ガウシアンを一様乱数でモンテカルロ積分した際の一様乱数の幅依存性

| 00:19 | ガウシアンを一様乱数でモンテカルロ積分した際の一様乱数の幅依存性を含むブックマーク

ガウシアンを、一様乱数モンテカルロ積分した際の結果の、一様乱数の範囲依存性がみたい。

こんな感じか。

f <- function(x){exp(-0.5*x^2)}
N <- 10^3
L <- c(1:10, 20, 50, 100)
y <- numeric(length(L))
for(i in 1:length(L)){
  y[i] <- L[i]/N*sum(f(runif(N, min=-0.5*L[i], max=0.5*L[i])))
}
plot(L, y, type="b", col=2)
abline(h=sqrt(2*pi))

横軸が一様乱数の幅で、縦軸が評価値。黒線が答えで赤線モンテカルロ積分評価した値。3σもあれば十分かなと思っていたが、6σくらいないと結果が収束してないように見える。また、幅がでかすぎてもいらないところを評価しまくることになるので、推定精度は下がる。

f:id:teramonagi:20150412002650p:image

ありがたい事に添削してもらった。RMSE(=sqrt(M)×エラーバー)も算出してもらっているので、こっちのほうがいいね&可視化ggplot2苦手勢にはありがたい。

2015-04-10

時系列データに対するブートストラップ法(ブロック・ブートストラップ法)について

| 19:42 | 時系列データに対するブートストラップ法(ブロック・ブートストラップ法)についてを含むブックマーク

あたまだし

検定やクロスバリデーション等への応用を企図した、サンプル数を水増しするための手法としてブートストラップ法がある。これをRで実行するにはsample関数を使って自分でリサンプリングするコードを実装するか、あるいはbootパッケージのboot関数を用いればいい。

ただ、通常このようなリサンプリングにおいては、例えば、データのレコードの行番号を一様にリサンプリングするなど、"データの順序"を考慮したものとはなっておらず、これはデータの順序に意味があるデータ、特に(時)系列データに対して問題となってくるので、通常のブートストラップ法を適用することはできない。

時系列データに対するブートストラップ法に関しては、まず、大枠としてのブロック・リサンプリング法があり、その構成要素としてブロック・ジャックナイフ法、ブロック・ブートストラップ法が研究されてきた。ブロック・ジャックナイフ法はさほどメジャーではないということなので、ここではブロック・ブートストラップ法のみを取り上げる。また、時系列データには定常性を仮定する。

詳しくは

を見るとよい。論文へのリファレンスも多くオススメ。本記事の内容は当該書籍のRのコード部分の要約に近い。

ブロック・ブートストラップ法

ブロック・ブートストラップ法のアイディアは非常に単純で、

  • データの順序に意味があるのならば、データを1つだけではなく、塊(ブロック)としてリサンプリングすればよかろう

というものである。そうすることで、大局的なデータの相関構造は壊れるが、ブロック内での局所的なデータの相関構造は担保されるだろうというものである。

ブロック・ブートストラップ法の中でも、ここでは

  • 移動ブロック・ブートストラップ法(Moving Block Bootstrap, 以下MBB法)
  • 円形ブロック・ブートストラップ法(Circular Block Bootstrap, 以下CBB法)
  • 定常ブートストラップ法(Stationary Bootstrap, 以下SB法)

の3つを取り上げる。上述の書籍ではこれに加え「非重複ブロック・ブートストラップ法」も取り喘げているが、私的にあまりイケてない感があったので、こいつは端折る。

以下ではbootパッケージのtsboot関数を使うので、これはインストールして読み込んでおく。

install.packages("boot")
library("boot")

移動ブロック・ブートストラップ法(Moving Block Bootstrap, MBB法)

MBB法は例えば、

> x <- 1:5
> x
[1] 1 2 3 4 5

というデータがあった時に、ブロックサイズを3とすると

c(1:3)
c(2:4)
c(3:5)

という3つのブロックにデータ分割し、このブロック自身をリサンプリングするものである。なお、ブロックは元のデータと同じ数だけデータがリサンプリングされるよう、リサンプリングされるその個数が決まる。今回の場合は2個ブロックがあればデータが6個できるので、ブロックは2個リサンプリングされ、はみ出した1個のデータは捨てられる。

tsboot関数を使ってMBB法を行うためには、引数として

  • sim="fixed"
  • endcorr=FALSE

と設定する。その他の引数については、第一引数(tseries)にデータを、第二引数(statistic)にブートストラップ法で計算したい統計量を算出したい関数を指定する。

ここで第二引数にはMBB法を使って計算したい統計量を算出する関数を渡すが、統計量ではなく、リサンプルされたサンプルそのものが欲しい場合には、引数をそのまま返す関数を指定すればよい。Rには引数をそのまま返す関数、force関数があるので、これを指定している。force関数の実装は要するに

force <- function(x){x}

ということだ。

実際にMBB法を使って見ると

> # ブロックの長さを3としたMBB法によりテストデータと同じ長さのブートストラップ標本を4回発生させる
> tsboot(x, force, R=4, l=3, sim="fixed", endcorr=FALSE)$t
[,1] [,2] [,3] [,4] [,5]
[1,]    2    3    4    2    3
[2,]    2    3    4    1    2
[3,]    2    3    4    1    2
[4,]    1    2    3    1    2

となる。使用した引数の意味を改めて書いておくと

  • tseries: (時系列)データ
  • stat: ブートストラップ法により計算したい統計
  • R: リサンプリング回数
  • l: ブロックの長さ
  • sim: リサンプルの発生方法
  • endcorr: ブロックの端点補正有無

となる。

tsboot関数の返り値自体はbootクラスのオブジェクトになっており、リサンプリングの際の条件なども持っているが、ここでは特段いらないので、その返却値に対して$tと付けることで、リサンプリングの結果だけを取得している。

円形ブロック・ブートストラップ法(Circular Block Bootstrap, CBB法)

これは物理でいう”周期境界条件”をつけたものと考えると良い。具体的には、先ほど同様のデータ

> x
[1] 1 2 3 4 5

に対して、ブロックサイズを3とすると

c(1:3)
c(2:4)
c(3:5)
c(4,5,1)
c(5,1,2)

という5つのブロックにデータ分割し、このブロック自身をリサンプリングするものである。こうすることで先ほどのデータでいうと1や5という数値が2,3に比べて出にくいというMBB法に存在した、"リサンプリングデータに対する偏り"がなくなる。

RでCBB法を実行する上で、MBB法と違う点はendcorr=TRUEとするだけである。また、このendcorr引数デフォルトでTRUEとなっているので、実際には書かなくていい。

# ブロックの長さを3としたCBB法によりテストデータと同じ長さのブートストラップ標本を4回発生させる
> tsboot(x, force, R=4, l=3, sim="fixed")$t
[,1] [,2] [,3] [,4] [,5]
[1,]    3    4    5    3    4
[2,]    3    4    5    2    3
[3,]    3    4    5    4    5
[4,]    3    4    5    1    2
定常ブートストラップ法(Stationary Bootstrap, SB法)

理論的な背景を調べきれていないのだが、元の時系列が定常であっても、MBB・CBBでリサンプリングした場合には定常性が崩れる場合がある(らしい)。そういう状況を避けるためには、SB法を使うとよい。RでSB法を実行する上で、CBB法と違う点はsim="geom"とするだけである。こうすることで、平均l(ブロックサイズ)になる幾何分布からブロックサイズをサンプリングするようになる。なので、SB法の場合、ブロック数は一定とならない。

# ブロックの長さを3としたSB法によりテストデータと同じ長さのブートストラップ標本を4回発生させる
> tsboot(x, force, R=4, l=3, sim="geom")$t
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    4    5    5    1
[2,]    5    1    1    1    2
[3,]    1    2    3    4    5
[4,]    1    2    3    4    5

ブロックサイズの選び方、

上述の参考書籍に「サブサンプリング法に基づくブロック長さの選択法」なるものが載っていた。疲れたから省略。

まとめ

結局定常ブートストラップ法使っておけばいいのか?なんかタラタラ書いてたら堅い日本語になっちゃった・・・

MSYS2のbash上にgitのブランチ(branch)名を表示する

| 12:56 | MSYS2のbash上にgitのブランチ(branch)名を表示するを含むブックマーク

最近のWindowsにおけるLinuxコマンド環境として(俺の中で)熱いのが、MSYS2なのですが、そいつのコンソール上で、作業中のgitのブランチ(branch)名を出すためには、MSYS2がデフォルトインストール先だとすると、.bashrcの中に

source /c/msys64/usr/share/git/completion/git-prompt.sh
source /c/msys64/usr/share/git/completion/git-completion.bash
GIT_PS1_SHOWDIRTYSTATE=true
export PS1='\[\033[32m\]\u@\h\[\033[00m\]:\[\033[1;33m\]\w\[\033[31m\]$(__git_ps1)\[\033[00m\]\$'

と書いておけばOKだ。PS1環境変数の中で指定している色については、好みに合わせて適当に変えたらいい。

2015-03-23

リサンプリング(復元抽出)で積分値を評価する

| 10:17 | リサンプリング(復元抽出)で積分値を評価するを含むブックマーク

の応用かつ、「これで間違ってねぇよな?」という自分への備忘録。インポータンスサンプリングや測度変換とは微妙に分母の形式が違ってくるので間違いやすいんだ。

数式の整理

以下の積分I

 I = ¥int_1^{L} w(x) f(x) dx

を計算することを考える。積分範囲は以下の関数例にあうように適当に設定してしまった。

これはモンテカルロ法により以下のように近似される。UNIFORMは一様分布のこと。

 I = (L-1)¥int_1^{L} ¥frac{1}{L-1} w(x) f(x) dx ¥sim ¥frac{L-1}{N} ¥sum_{i=1}^{N} w(x_i)f(x_i), ¥,¥,¥,¥,¥, x_i ¥sim UNIFORM(1, L)

ここで、積分Iの式を、以下のように式変形して考える。

 I = ¥frac{L-1}{N} sumw ¥sum_{i=1}^{N} ¥frac{w(x_i)}{sumw} f(x_i), ¥,¥,¥,¥,¥, sumw=¥sum_{i=1}^{N} w(x_i)

ここで定義した

 p(x_i)=¥frac{w(x_i)}{sumw}

は各サンプルiに対する重み(確率)とみなすことができるので、この値に比例するようにリサンプリングすると

 I = ¥frac{L-1}{N} sumw ¥frac{1}{N}¥sum_{i=1}^{N} f(x_i)=¥frac{L-1}{N} weight ¥sum_{i=1}^{N} f(x_i), ¥,¥,¥,¥,¥, x_i ¥sim p(x_i), ¥,¥,¥, weight=¥frac{sumw}{N}

として計算することができる*1


できる・・・はずだと信じてやってみる。

やってみる

#適当な関数w, f
w <- function(x){log(x)}
f <- function(x){exp(-x^2)}
#Lを4にしてどんなもんになるか描画

f:id:teramonagi:20150323101637p:image

L <- 4
x <- seq(1, L, length.out = 100)
plot(x, w(x)*f(x))
#ふつーのモンテカルロ法(1万乱数)
N <- 10^4
x <- runif(N, min=1, max=L)
(L-1)*sum(w(x)*f(x))/N
#際サンプリングして計算
index <- sample(N, replace=TRUE, prob=w(x))
weight <- sum(w(x))/N
(L-1)*sum(weight*f(x[index]))/N
#実際に積分させてみた答え
stats::integrate(function(x){w(x)*f(x)}, 1, L)

実行結果の箇所だけ示すと

> (L-1)*sum(w(x)*f(x))/N
[1] 0.03627271
> (L-1)*sum(weight*f(x[index]))/N
[1] 0.03558152
> stats::integrate(function(x){w(x)*f(x)}, 1, L)
0.03588273 with absolute error < 4.1e-12

というわけで、ちゃんと一致してそうだ。

補記

ここで定義した

 p(x_i)=¥frac{w(x_i)}{sumw}

は各サンプルiに対する重み(確率)とみなすことができるので、この値に比例するようにリサンプリングすると

・・・の流れは

 I = ¥int_1^{L} w(x) f(x) dx = (¥int_1^{L} w(x)dx) ¥int_1^{L} ¥frac{w(x)}{¥int_1^{L} w(x)dx} f(x) dx

=(¥int_1^{L} w(x)dx) ¥frac{1}{N}¥sum_i^N f(x_i), ¥,¥,¥,¥, x_i ¥sim ¥frac{w(x)}{¥int_1^{L} w(x)dx}

¥sim ¥left( ¥frac{L-1}{N} ¥sum_i^N w(z_i)¥right) ¥frac{1}{N}¥sum_i^N f(x_i), ¥,¥,¥,¥, x_i ¥sim ¥frac{w(x)}{¥left( ¥frac{L-1}{N} ¥sum_i^N w(z_i)¥right)}, z_i ¥sim UNIFORM(1, L)

¥sim ¥frac{L-1}{N} sumw ¥frac{1}{N}¥sum_i^N f(x_i), ¥,¥,¥,¥, x_i ¥sim ¥frac{w(x)}{¥frac{L-1}{N} sumw}, z_i ¥sim UNIFORM(1, L), sumw=¥sum_i^N w(z_i)

¥sim ¥frac{L-1}{N} weight ¥sum_i^N f(x_i), ¥,¥,¥,¥, x_i ¥sim ¥frac{w(x)}{¥frac{L-1}{N} sumw}, z_i ¥sim UNIFORM(1, L), sumw =¥sum_i^N w(z_i), weight = ¥frac{sumw}{N}

とした方が見通しがよさげ

*1:ここで、"二回Nで割る"操作が入っているのを一回だと勘違いしてひどい目にあった。そのためのメモ

2015-03-22

何回も使うリストの要素は変数に突っ込んでおいた方がいい

| 08:17 | 何回も使うリストの要素は変数に突っ込んでおいた方がいいを含むブックマーク

当たり前のお話ですが、何回もアクセスするなら代入コスト払っても、一時変数使った方がよさげ。

> x <- list(list(a=1:3, b=100:103), list(a=2:4, b=200:203))
> system.time(sapply(1:10^6, function(i){
+   force(x[[2]]$a[[2]])
+   force(x[[2]]$a[[2]])
+   force(x[[2]]$a[[2]])
+ }))
   user  system elapsed 
   6.95    0.00    6.95 
> system.time(sapply(1:10^6, function(i){
+   hoge <- x[[2]]$a[[2]]
+   force(hoge)
+   force(hoge)
+   force(hoge)
+ }))
   user  system elapsed 
   4.88    0.01    4.95 

リストのリストの場合、変数を2回以上使うなら代入した方が速かった。

アクセス一回のケース

> system.time(sapply(1:10^6, function(i){
+   force(x[[2]]$a[[2]])
+ }))
   user  system elapsed 
   3.12    0.00    3.13 
> system.time(sapply(1:10^6, function(i){
+   hoge <- x[[2]]$a[[2]]
+   force(hoge)
+ }))
   user  system elapsed 
   3.43    0.00    3.44 

アクセス二回のケース

> system.time(sapply(1:10^6, function(i){
+   force(x[[2]]$a[[2]])
+   force(x[[2]]$a[[2]])
+ }))
   user  system elapsed 
   4.80    0.00    4.83 
> system.time(sapply(1:10^6, function(i){
+   hoge <- x[[2]]$a[[2]]
+   force(hoge)
+   force(hoge)
+ }))
   user  system elapsed 
   3.84    0.00    3.87