Hatena::ブログ(Diary)

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

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 

2015-03-17

復元抽出のアルゴリズム

| 00:50 | 復元抽出のアルゴリズムを含むブックマーク

もくてき

粒子フィルタ(パーティクル・フィルタ)を実行する際には、粒子のウェイト(weight)に比例する確率でリサンプリングを実行する必要がある。そのためのアルゴリズムとコードを考えたい。ここでは手元にある各粒子のウェイトはK個の要素からなるベクトルだと仮定して、さらにそれを復元を許してN個リサンプリングするという状況で考える。

こんな面倒な状況を考えなくても、要するにこれは

  • いろんな色の球が入ってる壺から、1個適当に球を取り出して、その色をメモって、球を戻す

を複数回繰り返すことと同じで、高校生で習う確率の範囲で理解できる計算なわけだ。

アルゴリズム1(逆変換法)

逆変換法のアイディアを使って以下のようにするのが素朴なアイディアでコードも短い(後述)。

1: weightを確率に直し、その累積確率を計算し、これを{Qk, k=1,2,...K}とする

2: [0, 1]の実数乱数rを一個生成する

3: rを、k=1,2,..., Kまで順番に累積確率Qkと比較していき、r<Qkとなる最小のkを答えとし、weight[k]を返す

4: 2-3をN回繰り返す

アルゴリズム2(ソートを使う方法)

ちょいと頭を使うといかのようなアルゴリズムもできる。

1: weightを確率に直し、その累積確率を計算し、これを{Qk, k=1,2,...K}とする

2: [0, 1]の実数乱数rをN個生成し、{rn, n=1,2,...N}とする

3: {rn}∪{Qk}として、これをソートし、配列unionとする

4: ソートした配列unionの間に挟まってる"乱数rn"の個数を数える

5: 4で数えた個数ぶんだけ、該当するweightの要素を返却する


速度比較

上述のアルゴリズムのオーダーは、

なので、データ(K) or リサンプリング(N)数が多い場合には、アルゴリズム2の方が効率的であると考えられるので、それを確かめてみる。

Rで書くと以下のような感じか。乱数は一度に生成するなどの計算時間節約は入れてしまっている。

# アルゴリズム1
resampling1 <- function(weight, size){
  prob_cum <- cumsum(weight)/sum(weight)
  index <- sapply(runif(size), function(x)which(x < prob_cum)[1])
  weight[index]
}
# アルゴリズム2
resampling2 <- function(weight, size){
  prob_cum <- cumsum(weight)/sum(weight)
  size_weight <- length(weight)
  union <- c(prob_cum, runif(size))
  union_indexes <- order(union, decreasing=TRUE)
  index <- numeric(size)
  value <- union_indexes[1]
  counter <- 1
  for(union_index in union_indexes[-1]){
    if(union_index > size_weight){
      index[counter] <- value
      counter <- counter + 1
    }else{
      value <- union_index
    }
  }
  weight[index]
}

速度の比較結果は

> x <- sample(1:10^4)
> system.time(resampling1(x, 10^5))
   ユーザ   システム       経過  
      9.35       0.00       9.36 
> system.time(resampling2(x, 10^5))
   ユーザ   システム       経過  
      0.25       0.00       0.25 

というわけで、コードの長さとは裏腹にアルゴリズム2の方が遥かに速い。

・・・で、ここまで書いてから、Rに組み込みのsample関数のprob引数を与えれば、同じ処理ができることに気がついた。なので、アルゴリズム2とsample関数の計算速度を比べてみる。

> system.time(resampling2(x, 10^5))
   ユーザ   システム       経過  
      0.25       0.00       0.24 
> system.time(sample(x, 10^5, replace=TRUE, prob=1:10^4))
   ユーザ   システム       経過  
         0          0          0 

・・・Rの実装の方が圧倒的に速い!?!?!?!?!?

Rcppで書きなおす

きぃぃい、悔しい!!!負けるのは悔しいのでRcppでアルゴリズム2を書きなおしてみた。RcppとRcppArmadilloパッケージは

install.packages("Rcpp")
install.packages("RcppArmadillo")

として突っ込んでおいておく。

library(Rcpp)
sourceCpp(code='
#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
NumericVector resamplingCpp(arma::vec weight, int size){
  RNGScope scope;
  const unsigned int size_weight = weight.size();  
  const arma::vec prob_cum = cumsum(weight)/sum(weight);
  arma::vec unionset = arma::join_cols(prob_cum, as<arma::vec>(runif(size)));
  arma::uvec union_indexes = arma::sort_index(unionset, "descend");  
  arma::uvec index = arma::uvec(size);
  double value = union_indexes[0];
  int counter = 0;
  const arma::uvec::iterator union_indexes_begin = union_indexes.begin()+1;
  const arma::uvec::iterator union_indexes_end   = union_indexes.end();
  for(arma::uvec::iterator union_iterator=union_indexes_begin; union_iterator!=union_indexes_end; ++union_iterator){
    if(*union_iterator > (size_weight-1)){
      index[counter] = value;
      counter++;
    }else{
      value = *union_iterator;
    }
  }
  arma::vec result = weight.elem(index);
  return NumericVector(result.begin(), result.end());
}')

そして、このresamplingCpp関数と、アルゴリズム2(R実装)・sample関数の速度比較してみる。

> system.time(resampling2(x, 10^6))
   ユーザ   システム       経過  
      2.90       0.00       2.93 
> system.time(sample(x, 10^6, replace=TRUE, prob=1:10^4))
   ユーザ   システム       経過  
      0.03       0.02       0.04 
> system.time(resamplingCpp(x, 10^6))
   ユーザ   システム       経過  
      0.22       0.00       0.22 

うおお、なんてこった。元のコード(アルゴリズム2)よりも計算速度が10倍以上速くなったRcppのコードでも負けた…R恐るべし。俺のC++実装の問題な気もするのが、生のCっぽく書かくのはしんどいので、ここでおしまい。sample関数使えばいいや、もう。

答えの確認

ここで俺俺実装した結果が、ちゃんと同じ結果を返すことを確認しておく。

> weight <- sample(1:3)
> set.seed(100)
> x1 <- resampling1(weight, 10^2)
> set.seed(100)
> x2 <- resampling2(weight, 10^2)
> set.seed(100)
> x3 <- resamplingCpp(weight, 10^2)
> all(x1[order(x1)]==x2[order(x2)])
[1] TRUE
> all(x1[order(x1)]==x3[order(x3)])
[1] TRUE

上のように、乱数のシードを揃えれば確かに全ての結果が一致する。

sample関数については、実際に復元抽出した結果から確認しておく。

> val1 <- table(sample(1:5, 10^5, replace=TRUE, prob=1:5))
> val1/min(val1)
       1        2        3        4        5 
1.000000 2.026021 3.031165 4.059304 5.012103 

確かに、モンテカルロ法の誤差のレベルで、指定したウェイト(1,2,3,4,5)に比例した答えになる。

おまけ

ここで述べた復元抽出アルゴリズムの他にも、「Walkerのエイリアス法」や

の「第7章 乱数生成:不確実性をつくる」の7.3に載っている、ルーレットを使った方法があるらしい。こちらは気が向いたらそのうち。

そして、Rのsample関数はまさに「Walkerのエイリアス法」を用いている模様。

"If replace is true, Walker's alias method (Ripley, 1987) is used"

参考

2015-03-12

クロス集計〜公式:dplyr + tidyr = (xtabs|(f)table)

| 06:27 | クロス集計〜公式:dplyr + tidyr = (xtabs|(f)table)を含むブックマーク

特にやりたくはないんだけど、クロス集計をしなければならない状況がある。そんなときExcelを使ってもいいんだろうけど、レポーティングまで含めてRでやってしまいたい、あると思います。そんな時どうするかって話。

Rには既に

というクロス集計してくれる関数があるんだけど、こいつらの使い方覚えるのめんどくさい(毎度ググってる)し、特に他のデータハンドリング系関数との相性もよくはないので、全てdplyr&tidyrパッケージで済ませたい。

なので、以下、xtabs&table関数の結果と同じになるように、dplyr&tidyrで書いてみる。

使用するパッケージ

当然、dplyr&tidyrは使う。

library(dplyr)
library(tidyr)

dplyr, tidyrの使い方については

がよくまとまっていて、大変ありがたい。

使用するデータ

なんかよくわからんけどxtabsの例に載ってたデータを使ってみることした。適当なデータフレームですね。

> df <- as.data.frame(UCBAdmissions)
> head(df)
     Admit Gender Dept Freq
1 Admitted   Male    A  512
2 Rejected   Male    A  313
3 Admitted Female    A   89
4 Rejected Female    A   19
5 Admitted   Male    B  353
6 Rejected   Male    B  207

公式1:xtabs(~X+Y) = group_by(X, Y) %>% tally %>% spread(X, n)

クロス集計において、その”キーの組み合わせ数”をカウントするコードは以下。

> df %>%
+   xtabs(~Gender+Admit, .)
        Admit
Gender   Admitted Rejected
  Male          6        6
  Female        6        6
> 
> df %>% 
+   group_by(Gender, Admit) %>%
+   tally %>%
+   spread(Admit, n)
Source: local data frame [2 x 3]

  Gender Admitted Rejected
1   Male        6        6
2 Female        6        6
公式2: xtabs(Z~X+Y) = group_by(X, Y) %>% summarize(n=sum(Z)) %>% spread(X, n)

クロス集計において、”各キーの○○要素の合算値”を計算するコードは以下。

> df %>%
+   xtabs(Freq~Gender+Admit, .)
        Admit
Gender   Admitted Rejected
  Male       1198     1493
  Female      557     1278
>
> df %>% 
+   group_by(Gender, Admit) %>%
+   summarize(n=sum(Freq)) %>%
+   spread(Admit, n)
Source: local data frame [2 x 3]

  Gender Admitted Rejected
1   Male     1198     1493
2 Female      557     1278
公式3: table(X, Y)= group_by(X, Y) %>% tally %>% spread(Y, n)

マニュアルにあるtable関数の例に合わせて、ここでは、別なデータに対して実行。出力結果がちょいと違うけど、まぁいいか&NAはあとで0埋めしよう。

> b <- factor(rep(LETTERS[1:3], 10))
> d <- factor(rep(LETTERS[1:3], 10), levels=LETTERS[1:5])
> b
 [1] A B C A B C A B C A B C A B C A B C A B C A B C A B C A B C
Levels: A B C
> d
 [1] A B C A B C A B C A B C A B C A B C A B C A B C A B C A B C
Levels: A B C D E
> table(b, d)
   d
b    A  B  C  D  E
  A 10  0  0  0  0
  B  0 10  0  0  0
  C  0  0 10  0  0
> df <- data.frame(b, d)
> df %>%
+   group_by(b, d) %>%
+   tally %>%
+   spread(b, n)
Source: local data frame [3 x 4]

  d  A  B  C
1 A 10 NA NA
2 B NA 10 NA
3 C NA NA 10

結論

総じてdplyr&tidyrの方がタイプ数が多くなるけど、余計なクラス構造&返り値を覚えなくていいので、私はdplyr&tidyrでいきたい。

hoxo_mhoxo_m 2015/03/12 13:42 tally よりも count を使うほうが楽かもしれません。
また、spread には fill という引数があり、NA に対するデフォルト値を指定できます。

3つ目の例はこうなります。

>|r|
df %>%
count(b, d) %>%
spread(b, n, fill=0)
||<

>||
d A B C
1 A 10 0 0
2 B 0 10 0
3 C 0 0 10
||<

enjoy!

teramonagiteramonagi 2015/03/12 18:35 おお、より楽になりそうですね、ありがとうございます!

hoxo_mhoxo_m 2015/03/13 10:25 spread には drop という引数もありますね。

df %>%
count(b, d) %>%
spread(d, n, drop = FALSE, fill = 0)

b A B C D E
1 A 10 0 0 0 0
2 B 0 10 0 0 0
3 C 0 0 10 0 0

2015-02-28

RmdファイルからRコードを抜くのはpurl関数

| 10:18 | RmdファイルからRコードを抜くのはpurl関数を含むブックマーク

例えば、以下のような内容のR Markdownファイル(hoge.Rmd)を用意していたとする。

げっへっへ、本文だぜー!!
```{r, cache=TRUE}
#Rのコード
x <- 1:10
```
また、本文だぜー!

ここから、Rのコードだけを抽出するためにはknitrパッケージにあるpurl関数を用いると良い。documentation引数に応じて出力が変わる(チャンク・オプション、本文のコメント化など)ので、以下にその結果をメモる。

documentation=0

コマンド

purl("hoge.Rmd", documentation=0)

結果


#Rのコード
x <- 1:10

documentation=1(デフォルト)

コマンド

purl("hoge.Rmd", documentation=1)

結果


## ----, cache=TRUE--------------------------------------------------------
#Rのコード
x <- 1:10

documentation=2

コマンド

purl("hoge.Rmd", documentation=2)

結果

#' げっへっへ、本文だぜー!!
## ----, cache=TRUE--------------------------------------------------------
#Rのコード
x <- 1:10

#' また、本文だぜー!