Hatena::ブログ(Diary)

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

2016-08-17

三連ドット(..., dot-dot-dot, ellipsis)の取り扱い

| 16:22 | 三連ドット(..., dot-dot-dot, ellipsis)の取り扱いを含むブックマーク

これもr-wakalangに投げ込んで教えてもらった話なので、まずは簡単にまとめる。

基本的な使い方

適当な...を持つ関数を定義する。

f1 <- function(x, ...)
{
  dots <- list(...)
  print(dots)
}

これに対して以下の実行結果からわかるように、引数にマッチしなかったもの(ここではx以外)がリストのdots変数として関数内で使えるようになっているのがわかる。

要するに、多言語でいうところの可変長変数のようなもんだ。

list(...)という書き方の原点がどこにあるのかは不明だが、下記参考資料にあるR Language Definitionにも載ってるし、まぁこれはこういうもんかと思っておく。

> f1("a", 3)
[[1]]
[1] 3

> f1("b", a=3, b=7)
$a
[1] 3

$b
[1] 7

...のマッチに関する注意

普通のプログラミング言語だと可変長引数の後にはなんも書けない(引数を追加できない)ような気がするが、Rではできてしまう。

f2 <- function(x, ..., y)
{
  dots <- list(...)
  print(dots)
  print(paste0("y: ", y))
}

この場合、x, y以外の変数が全部...にマッチされている。

> f2("a", 3, z=7, y=3)
[[1]]
[1] 3

$z
[1] 7

[1] "y: 3"

R言語徹底解説にあったようなお話(data.frameに対するfilter操作)

R言語徹底解説のどこだかに載っていた気がする「data.frameに対するfilter操作」をdplyrぽく書く話も以下のように実現できる。

詳しくはoreorefilter関数のコメント参照。

oreorefilter <- function(df, ...)
{
  # 評価前に言語オブジェクトにしちゃう
  object <- as.list(substitute(list(...)))
  # どういう型として内部で保持されているかの確認
  print(str(object))
  # as.listしない場合の構造(つながった1リストになっちゃってる)
  print(str(substitute(list(...))))
  # 条件の順次評価(再帰で書いたほうがスマートか?)、一個目の要素はいらないのでループに入れない
  condition <- rep(TRUE, nrow(df))
  for(obj in object[-1]){
    # 条件を順次df環境で評価し、AND(&)条件として結ぶ
    condition <- condition & eval(obj, df)
  }
  df[condition, ]
}

実際に使ってみると、ちゃんとirisデータをFilteringできていることが分かる。

> oreorefilter(iris, Species=="setosa", Petal.Width==0.2, Sepal.Length==5.0)
List of 4
 $ : symbol list
 $ : language Species == "setosa"
 $ : language Petal.Width == 0.2
 $ : language Sepal.Length == 5
NULL
 language list(Species == "setosa", Petal.Width == 0.2, Sepal.Length == 5)
NULL
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5             5         3.6          1.4         0.2  setosa
8             5         3.4          1.5         0.2  setosa
26            5         3.0          1.6         0.2  setosa
36            5         3.2          1.2         0.2  setosa
50            5         3.3          1.4         0.2  setosa

参考

2016-08-16

金融市場が混乱する間隔を日経平均ボラティリティー・インデックス(日経平均VI)から推定したい

| 18:48 | 金融市場が混乱する間隔を日経平均ボラティリティー・インデックス(日経平均VI)から推定したいを含むブックマーク

日経ボラティリティー・インデックス日経平均VI)祭りはまだ続く。

やりたいことは簡単で

だと考えて、"日経VIが跳ねるタイミング"の確率分布がどのような形になるのかを推定するというお話。

個人的には指数分布なんじゃね?と思っていたが、違った。

データの準備

ほぼ

でやった内容とおなじで、日経VIが跳ねたポイント(下図赤点)を手でやるのがめんどうなので、Rに同定させるためのコード。

↑のLINKとはrollapp()関数がちょっと違ってるので注意

library("Quandl")
# Sorce: https://www.quandl.com/data/NIKKEI/VLTL-Nikkei-Stock-Average-Volatility-Index
# Download Nikkei VI
vi <- Quandl("NIKKEI/VLTL")
#Decay detect
rollapp <- function(x, k, fun)
{
  w <- k/2
  sapply(seq_len(length(x)), function(i){
    index_s <- max(i-w, 1)
    index_e <- min(i+w, length(x))
    middle <- index_s + i-w
    fun(1+i-index_s, x[index_s:index_e])
  })
}
x <- vi$Close
max_point <- rollapp(x, 74, function(i, x){x[i]==max(x)})
max_point <- replace(max_point, is.na(max_point), FALSE)
plot(x, type="l")
points(seq_len(length(x))[max_point], x[max_point], col = "red", pch=16)

f:id:teramonagi:20160816185327p:image

このデータの赤点同士の時間間隔がどういう確率分布に従っているのかを評価したい。

推定(指数分布)

指数分布で確率密度関数を推定してみる

# 適当にdata.frame化
df <- dplyr::group_by(data.frame(index=cumsum(max_point)), index) %>% summarize(n=n())
lambda <- 1/mean(df$n)
# 下記でも同じ
# library("MASS")
# exponential <- coef(fitdistr(df$n, "exponential"))
hist(df$n, freq=FALSE)
curve(dexp(x, rate=lambda),add=TRUE,lty=2)

f:id:teramonagi:20160816185326p:image

う、うーん、これではない感。

推定(ワイブル分布)

によると、ワイブル分布は「時間に対する劣化現象や寿命を統計的に記述するためにも利用される」とあるので、こいつでやってみたい。

金融市場が破裂する=寿命が尽きると考えるわけだ。

library("MASS")
weibull <- coef(fitdistr(df$n, "weibull"))
hist(df$n, freq=FALSE)
curve(dweibull(x, weibull["shape"], weibull["scale"]), add=TRUE, lty=2)

f:id:teramonagi:20160816185325p:image

これだなー!

あとはこの裏側の寿命強度過程的なものがどうなっているのかを調べれば大体OKだな。

参考

日経平均ボラティリティー・インデックス(日経平均VI)の半減期が知りたい

| 15:26 | 日経平均ボラティリティー・インデックス(日経平均VI)の半減期が知りたいを含むブックマーク

引き続き日経VIネタ。

日経VIの挙動は、株価のような対数幾何分布に従うような確率過程(そうは見えないけどそうだと思っておく!)ではなく、平均回帰的であり、上昇するも下落するも急なものであるというのが(私的に)一般的に知られた事実である。

ここでは日経VIがどのくらいの速さで平均回帰するのか?、すなわち急上昇したのちにどれくらいの日数をかけて落ちるのかについて考えたい。

物理学では元の値の半分の大きさ(すなわち、金融でいうところの半値に相当)になるまでにかかる時間(日数)を半減期(とか、特徴的な長さスケール)だなんて呼ぶのでそれを計算してやる。

これがわかるとおおよそ日経VIがピークになってからその半分の値になるまでの時間(日数)がわかるのである。

データの取得

データは従前の記事同様、Quandlパッケージを活用して取得する。

library("Quandl")
# Sorce: https://www.quandl.com/data/NIKKEI/VLTL-Nikkei-Stock-Average-Volatility-Index
# Download Nikkei VI
vi <- Quandl("NIKKEI/VLTL")

ピーク位置の同定

これが今回の一番悩んだところ。

ここでは以下のアイディアに基づいて抽出しているが、これが一般的なのかどうかはさっぱりわからない*1

後から調べた感じ、

この辺を参考にパッケージ使えばよかった…。


ここではまず、時系列データの各点を中心として半径(単系列での半径なので、ただの一次元円の半径、すなわちただの距離)k/2以内にあるデータに対してfunという処理をかませるrollapp()関数を定義する。

これは、時系列の一番後ろ(時間的に昔)からデータを舐めて処理するzooパッケージにあるrollapplyとは違い、各時点の未来・過去k/2単位(今回の場合日次データなので、単位は日)にまたがったデータに対して処理を適用する関数

#Decay detect
rollapp <- function(x, k, fun)
{
  w <- k/2
  sapply(seq_len(length(x)), function(i){
    if((i-w) <= 0 || (i+w) > length(x)){
      NA
    } else{
      fun(x[(i-w):(i+w)])
    }
  })
}

こいつを使って

  • 上値ピークポイント:直近150日(未来・過去方向それぞれ75日)分のデータの中で最も大きく、かつ、日経VIの値自体は30以上
  • 下値ピークポイント:直近125日(未来・過去方向それぞれ62日)分のデータの中で最も小さい

点を同定し、それぞれmax_point, min_pointとする。

NAはよしなに削除で。

x <- vi$Close
max_point <- rollapp(x, 150, function(x){x_middle <- x[1+length(x)/2]; (x_middle==max(x)) && (x_middle > 30)})
max_point <- replace(max_point, is.na(max_point), FALSE)
min_point <- rollapp(x, 125, function(x){x_middle <- x[1+length(x)/2]; (x_middle==min(x))})
min_point <- replace(min_point, is.na(min_point), FALSE)

ここで計算したピークポイントがどの辺にあるのかというと、下図のとおり、まぁ見た感じよさそうだぞと。

(上(下)値ピークが赤(青)点、日経VIの値は黒線で表示)

(逆に、こういう見た目になるように↑の日数とか日経VIの水準を足きりしている点に注意)

plot(x, type="l")
points(seq_len(length(x))[max_point], x[max_point], col = "red", pch=16)
points(seq_len(length(x))[min_point], x[min_point], col = "blue", pch=16)

f:id:teramonagi:20160816153707p:image

上昇→下落フェーズのデータの抽出

ここで、「定義した上値・下値の間にあるデータ」のみを抽出する。

つまり上の図でいうところの、各赤点・青点の間にある日経VIのデータを抽出するということだ。

そのためには汚いが、データを全舐めして、下記のようなコードでデータを抽出する。

index変数が0以外の値で、かつ、共通の値を持つ箇所がおなじ上昇→下落局面に対応している。

index <- numeric(length(x))
counter <- 0
mode <- "max"
for(i in seq_len(length(x))){
  if(mode == "max" && max_point[i]){
    mode <- "min"
    counter <- counter + 1
    index[i] <- counter
  } else if(mode == "min" && min_point[i]){
    mode <- "max"
  } else if(mode == "min"){
    index[i] <- counter
  }
}
xs <- lapply(seq_len(max(index)), function(i){x[which(index == i)]})

結果は、各赤点・青点の間にあるデータごとのリストになっている。

半減期の同定

上で作ったリストの各要素が上昇→下落局面に対応する

 NIKKEI-VI = a*¥exp¥{-t/¥tau¥}

という回帰式に対してτを推定することがメインであり、このτが半減期、すなわち日経VIがどの程度の速度で下落するのかの指標になっている*2

推定する部分を関数にし、元データ・τ・推定結果(データ自体をlog化し、そのデータの誤差が正規分布になってほしいなと仮定して、線形回帰に直して計算させる)を返却させるようにしよう。

# x = a*exp(-t/tau)
# log(x) = log(a) - 1/tau * t
decay <- function(x) {
  df <- data.frame(t=seq_len(length(x)), logx=log(x))
  regression <- lm(logx~t, data=df)
  tau <- -1/coef(regression)["t"]
  list(x=x, tau=tau, regression=regression)
}

これを、上で作った全データに適用してやると

result <- lapply(xs, decay)

となる。

得られたτは

> tau <- sapply(result, function(x){x$tau})
> tau
        t         t         t         t         t         t         t         t         t 
 56.60291  58.95572 165.28740 395.43998 201.99652 109.37092  40.24398  33.13756  27.06593 

となり、30くらいから400くらいまでとかなり幅広だ。

各々の結果を出したデータと回帰の当てはまり具合を下記のように目視で確認すると、長い半減期の系はだいたいフィッティングがおかしいので、無視してよさそうだ。

これは、どちらかというと、データの中に小さなピークが立っていることからもわかるように、ピークの同定がうまくないからなんだろうなと。

従って、短い方を意識して大体、半減期は10-30日程度と意識して日経VI先物投資すればよいだろう(と期待したい)。

τ=56.60291

f:id:teramonagi:20160816160309p:image

τ=58.95572

f:id:teramonagi:20160816160308p:image

τ=165.28740

f:id:teramonagi:20160816160307p:image

τ=395.43998

f:id:teramonagi:20160816160306p:image

τ=201.99652

f:id:teramonagi:20160816160304p:image

τ=109.37092

f:id:teramonagi:20160816160303p:image

τ=40.24398

f:id:teramonagi:20160816160302p:image

τ=33.13756

f:id:teramonagi:20160816160301p:image

τ=27.06593

f:id:teramonagi:20160816160300p:image

*1:まじめに調べてもいない

*2:より正確にはこの値の0.7倍程度か https://ja.wikipedia.org/wiki/%E5%8D%8A%E6%B8%9B%E6%9C%9F

2016-08-15

dplyrで時系列データの処理

| 22:01 | dplyrで時系列データの処理を含むブックマーク

で書いた、”dplyrでの時系列処理ってどうするんじゃい”というのの簡単な例。

下で示したようなやり方で大半は片付くのではないだろうか(たぶん)。

まず、ライブラリサンプルデータの作成。

lubridateパッケージが大事で、こいつで日付を操作しつつdata.frameのままdplyrで裁いてデータ処理をしていくという方針である。

library(xts)
library(lubridate)
library(dplyr)
# 適当なxtsデータ
x <- as.data.frame(sample_matrix)
x$index <- as.Date(rownames(sample_matrix))
> head(x)
               Open     High      Low    Close      index
2007-01-02 50.03978 50.11778 49.95041 50.11778 2007-01-02
2007-01-03 50.23050 50.42188 50.23050 50.39767 2007-01-03
2007-01-04 50.42096 50.42096 50.26414 50.33236 2007-01-04
2007-01-05 50.37347 50.37347 50.22103 50.33459 2007-01-05
2007-01-06 50.24433 50.24433 50.11121 50.18112 2007-01-06
2007-01-07 50.13211 50.21561 49.99185 49.99185 2007-01-07

では、このxというデータが月ごとに何件ずつデータをもっているのかをdplyrで集計してみると

> # 月次の集計
> x %>% group_by(month=month(x$index)) %>% summarize(n())
Source: local data frame [6 x 2]

  month   n()
  (dbl) (int)
1     1    30
2     2    28
3     3    31
4     4    30
5     5    31
6     6    30

となる。

このノリで全部いいんじゃないかな、xtsとかもういらない子なじゃないかなと思ってる。

dplyrを拡張してデータフレームっぽい俺俺クラスを動かす

| 20:34 | dplyrを拡張してデータフレームっぽい俺俺クラスを動かすを含むブックマーク

みんな大好きデータフレームのハンドリング用ライブラリdplyrを拡張して、時系列データのxts型も同じように操作できるtplyrパッケージを作ろうと思ったんだが、よくよく考えるとdplyrままでいいんじゃないかなって思って辞めた。

その際に試行錯誤した残骸を忘れないようにここにまとめる。

ここでは、dplyrのfilter関数に対してxts型が動作するように拡張することを考えよう。

そのために、以下のように3つの関数を用意する。

内容はコメントにあるとおりxtsとdata.frameの変換関数とxtsに対する実際のfilter処理だ。

# xtsをdata.frameへと変換する関数
xts_to_df <- function(x)
{
  data.frame(index=zoo::index(x), zoo::coredata(x))
}
# data.frameをxtsへと変換する関数
df_to_xts <- function(df)
{
  xts(dplyr::select(df, -index), order.by=df$index)
}
# xts型に対するfilter処理を記述(S3クラスのメソッドとする)
filter_.xts <- function(.data, ..., .dots) {
  dots <- lazyeval::all_dots(.dots, ...)
  df <- xts_to_df(.data)
  df_to_xts(dplyr::filter_(df, .dots = dots))
}

そして、ここがミソだが、ここで作成したfilter_.xts関数を無理やりdplyrのnamespaceに突っ込む。

こうすることで、dplyrパッケージの中の処理において、xts型のfilter処理に対してこの関数を呼んでくれるようになる*1

environment(filter_.xts) <- asNamespace("dplyr")

さて、実際に動くかxtsパッケージにあるサンプルデータで確かめてみよう。

> library(xts)
> data(sample_matrix)
> # 適当なxtsデータへと変換
> samplexts <- as.xts(sample_matrix)
> # 内容の確認
> head(samplexts)
               Open     High      Low    Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185
> # ちゃんと動いてるか確認
> dplyr::filter(samplexts, Close > 51)
               Open     High      Low    Close
2007-02-14 50.95283 51.04699 50.80317 51.04699
2007-02-15 51.06330 51.11401 50.94681 51.05185
2007-02-16 51.12879 51.12879 51.00613 51.02164
2007-02-17 50.97722 51.13653 50.95260 51.13653
2007-02-18 51.18414 51.32090 51.13713 51.15151
2007-02-19 51.29502 51.32342 51.13524 51.17899

・・・ちゃんとfilter処理できてる!やったぜ、親分

パッケージ作ってて、そのパッケージの中のある関数中で他のパッケージの関数呼ぶ場合、DESCRIPTIONのImports書いてるだけじゃだめぽなんですが、 これちゃんとNAMESPACEにImport書かないとダメなんでしたっけ…?

| 10:08 | パッケージ作ってて、そのパッケージの中のある関数中で他のパッケージの関数呼ぶ場合、DESCRIPTIONのImports書いてるだけじゃだめぽなんですが、 これちゃんとNAMESPACEにImport書かないとダメなんでしたっけ…?を含むブックマーク

完全備忘録

hoxo_n パッケージ作ってて、そのパッケージの中のある関数中で他のパッケージの関数呼ぶ場合、DESCRIPTIONのImports書いてるだけじゃだめぽなんですが、

これちゃんとNAMESPACEにImport書かないとダメなんでしたっけ…?

聖書にちゃんと書いてありました。。。

http://r-pkgs.had.co.nz/namespace.html#imports

```R functions

If you are using just a few functions from another package, my recommendation is to note the package name in the Imports: field of the DESCRIPTION file and call the function(s) explicitly using ::, e.g., pkg::fun().

If you are using functions repeatedly, you can avoid :: by importing the function with @importFrom pgk fun. This also has a small performance benefit, because :: adds approximately 5 µs to function evaluation time.

Alternatively, if you are repeatedly using many functions from another package, you can import all of them using @import package. This is the least recommended solution because it makes your code harder to read (you can't tell where a function is coming from), and if you @import many packages, it increases the chance of conflicting function names.

```

hoxo_m

```#' @importFrom package method

```

みたいに書いておくと roxygen がよしなにやってくれます

hoxo_n 親分は packagename::methodじゃなくてImportしちゃうマンかー

これ使用している関数の上に書いておかないと収集つかなくなりそうで、それならパッケージ名毎回書いてもいいかなって思いだしてました

あ、そうかいてるのね、おれもそうしよう

https://github.com/hoxo-m/bayesopt/blob/master/R/bayesopt.R

johndoe たまに書くと忘れているので、忘れていたら@hoxo_mか@hoxo_nのコードを読むと心に刻みました

hoxo_n おれもわすれているし、Pandasの勉強をしてモテたいです。

hoxo_m 昔は `pkg::fun()` マンだったけど CRAN 通すときにこの書き方知ってこっち派になりました。

hoxo_n Why are you usin `pkg::fun()` ? 事案だ (edited)

hoxo_m わかんないけど CRAN の中の人にこっちの書き方指定されたんですよねー。どっちでもいいのかもしれんけど。

hoxo_n よくわからないことに安易に迎合しない!それが俺たち(略

johndoe さきほどのhoxo_nの聖書からの引用に少し速くなるよ(マイクロ秒レベル)って書いてるのでそれくらいの利点と認識しました

hoxo_m 一つだけわかるのは `pkg::fun()` の書き方だと ` .onAttach()` が呼ばれないと思うんですよね。例えば rstan なんかは `.onAttach()` でいろいろ前処理やってるので `rstan::stan()` みたいな書き方だととうまく動かないと思います。 (edited)

johndoe みなさんパスをの文字列操作時は、paste()とかsprintf()じゃなくてfile.path()を使っていますか?

hoxo_m 昔は `file.path()` 使ってましたけど、最近はあんまりこだわってないです (edited)

*1:と思ったが、グローバル環境に定義してる分にはこれは不要で、パッケージ化してこれをやろうとすると必要になるんだった、失敬失敬。

2016-08-14

日経平均ボラティリティー・インデックス(日経平均VI)を高速フーリエ変換したが、なんもなかった

| 18:48 | 日経平均ボラティリティー・インデックス(日経平均VI)を高速フーリエ変換したが、なんもなかったを含むブックマーク

なんもなかった。が、なんもなかったなりの備忘録

データの取得

Quandle経由で取得。元のソースはコメントに記載。

PLOTを見る限り、ちゃんとした元の値っぽい。

library("Quandl")
# Sorce: https://www.quandl.com/data/NIKKEI/VLTL-Nikkei-Stock-Average-Volatility-Index
# Download Nikkei VI
vi <- Quandl("NIKKEI/VLTL")
# Plot
plot(vi$Date, vi$Close)

f:id:teramonagi:20160814185610p:image:w640

データをスケーリング(平均を0化のみ)し、FFTかけてPLOT。

データ自体はFFTの掟的にナイキスト周波数以降はぶった切ってる。

この図から見るに、低周波数=長周期が支配的っぽい。

# Scaling
x <- vi$Close - mean(vi$Close)
# FFT and convert into powerspectrum
vi_fft <- (abs(fft(x))^2)[1:(nrow(vi)/2)]
# Frequency range
f <- (seq(0, nrow(vi)-1, 1)/(nrow(vi)))[1:(nrow(vi)/2)]
#Plot frequency v.s. power spectrum(1)
plot(f, vi_fft)

f:id:teramonagi:20160814185412p:image:w640

↑が見づらいので、高周波数領域をやや削る。

#Plot frequency v.s. power spectrum(2)
plot(f[1:100], vi_fft[1:100])

f:id:teramonagi:20160814185411p:image:w640

周波数だとわかりにくいので、周波数→時間(周期)での表示に変換

500日強の周期を持つ波が強そうだが、これはノイズのlevel?&トレーディング的にそんなに待てないので、放置で。。。

#Plot frequency v.s. power spectrum
plot(1/f[1:100], vi_fft[1:100])

f:id:teramonagi:20160814185410p:image:w640

もうちょっと周期の短いところだけをPLOT。

う、うーん解釈に苦しむスペクトルだ。

#Plot frequency v.s. power spectrum
plot(1/f[5:100], vi_fft[5:100])

f:id:teramonagi:20160814185409p:image:w640

おしまい。

FFTのチェック

横軸の扱いがちょっと怪しいので確認。

周期123の正弦波のピークがちゃんとそこに立ってるかを確認(図示)。

# FFT check
T <- 123
t <- (0:500)
f <- (0:500)/500
wave = sin(2*pi*t/T)
plot(1/f[1:50], abs(fft(wave))[1:50])

f:id:teramonagi:20160814185413p:image

ついでに値も念のため確認。おおよそ対応する日にち(125日)にピークがあるので、あってるだろと。

> 1/f[1:50]
 [1]       Inf 500.00000 250.00000 166.66667 125.00000 100.00000  83.33333  71.42857  62.50000  55.55556  50.00000  45.45455
[13]  41.66667  38.46154  35.71429  33.33333  31.25000  29.41176  27.77778  26.31579  25.00000  23.80952  22.72727  21.73913
[25]  20.83333  20.00000  19.23077  18.51852  17.85714  17.24138  16.66667  16.12903  15.62500  15.15152  14.70588  14.28571
[37]  13.88889  13.51351  13.15789  12.82051  12.50000  12.19512  11.90476  11.62791  11.36364  11.11111  10.86957  10.63830
[49]  10.41667  10.20408
> abs(fft(wave))[1:50]
 [1]   1.8103671   2.9863890   6.1342361  14.6075087 246.2353106  21.4560098  11.1087158   7.7422349   6.0391541   4.9956202
[11]   4.2832646   3.7621472   3.3622089   3.0443135   2.7847859   2.5684082   2.3849216   2.2271388   2.0898629   1.9692353
[21]   1.8623257   1.7668661   1.6810707   1.6035134   1.5330400   1.4687057   1.4097291   1.3554581   1.3053438   1.2589205
[31]   1.2157905   1.1756122   1.1380900   1.1029669   1.0700187   1.0390482   1.0098819   0.9823661   0.9563642   0.9317544
[41]   0.9084279   0.8862869   0.8652434   0.8452180   0.8261388   0.8079406   0.7905640   0.7739551   0.7580643   0.7428466

参考

2016-08-03

foreachでprogressbarを表示しても進行状況は正しくないケースが多い

| 18:08 | foreachでprogressbarを表示しても進行状況は正しくないケースが多いを含むブックマーク

の辺でやられているforeach系パッケージでプログレスバーを出す方法だが、これは進捗を正しく表すものではないケースが多い*1

例えば、

library(foreach)
library(utils)
library(iterators)
library(doParallel)
cl <- makeCluster(4, type='SOCK')
registerDoParallel(cl)
f <- function(){
  pb <- txtProgressBar(min=1, max=n-1,style=3)
  count <- 0
  function(...) {
    count <<- count + length(list(...)) - 1
    setTxtProgressBar(pb,count)
    flush.console()
    c(...)
  }
}

n <- 100
# Run the loop in parallel
k <- foreach(i = seq(n), .final=c, .combine=f()) %dopar% {
  a <- rnorm(10^7)
  log2(i)
}

というコードを実行するとわかるが、全体の計算が終わった後、一気にプログレスバーが進むような挙動になり、

Combine(並列化したデータをまとめあげる処理)が全スレッド・全計算が終わってから走ってるような挙動に見える。


この辺は、下記のStackOverflowでも議論されており、

It won't work well with doParallel because doParallel only calls the combine function after all of the results have been returned, since it is implemented by calling the parallel clusterApplyLB function. This technique only with works well with backends that call the combine function on-the-fly, like doRedis, doMPI, doNWS, and (defunct?) doSMP.

という記述がある通り。

私もこれが正しいなと思い、やはりdoParallelがバックエンドだとコンバイン処理が最後にまとまってくるという解釈でよくプログレスバーを正しく出すのは難しいぞと*2

*1:並列化のバックエンド依存があるので、こういう物言いになっとる

*2:いう話を https://r-wakalang.slack.com/messages/r_chotowakaru/ でしていたのでそのまとめ