餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

Rの並列処理(doParallel/foreach)で繰り返し処理を中断する方法

確認してみたのですが、試行錯誤してみた結果、結論は「ありません」でした*1。仕方が無いので、OpenMP 4.0以前の代替措置*2と同じ事をしてみたいと思います。速度的には改善になるはずです。

library("doParallel")

cl <- makeCluster(detectCores())
registerDoParallel(cl)

f <- TRUE
# 処理を if(f){ ... } で囲み、f==TRUEの間だけ計算するようにする
r <- foreach(i=1:100, .combine='c') %dopar% if(f){
# √計算するだけの処理
  v <- sqrt(i)
  if(v > 6){
# 目標達成をしたら f にFALSEを入れる
    f <- FALSE
  }
  v
}

# 平行処理終了
stopCluster(cl)

# 結果を確認する。f <- FALSE をコメントアウトして比較すると違いが分かるはず。
r

1から100までのルート計算をしていくループですが、ループは止められないけれども、6より大きい結果が得られたら計算自体はスキップすると言う代替です。ループ内の処理が重い場合は、高速化になります。

%dopar% if(f) を %:% when(f) %dopar% に書き換えたくなるのですが、平行処理スコープから親環境の変数は書き換え不可なので、狙った通りには動きません。f <- FALSE を assign("f", FALSE, envir=.GlobalEnv) に変えてみたりしたのですが、狙ったとおりには動きません。どうも平行処理スレッドごとに .GlobalEnv もコピーして保持しているようで、シェアードしないみたいですね。

*1: 初出と記事内容が大きく変わっていましたが、大きく勘違いしていたので訂正しています。動作検証時に %dopar% を %do% と書いてしまうポカをしていました。

*2: OpenMP並列forループの中断処理 - yohhoyの日記

プレビュー画面で特定情報を出力しないためのBloggerテンプレートの書き方

Bloggerで外部JavaScriptを読み込ませるとプレビュー画面だけ妙に時間がかかる謎現象が発生するので、その回避のために、プレビュー画面だけそこの部分の出力をしないことにしました。Bloggerの記事のプレビューを爆速にするカスタマイズ | クラウド番外地を参照して、以下のようにします。ちょっと変えてありますが、not (A and B) は、(not A) or (not B) なのでほぼ同じです。

<b:if cond='data:blog.pageType != &quot;item&quot; or data:blog.url != data:blog.homepageUrl'>

<script src='...' type='text/javascript'/>

</b:if>

何分もかかっていたのが、あっさり表示されるようになりました。とても謎です。

JavaでOpenMPもどき

気づくとマルチコア化が一般化したため、数値演算をするときにカジュアルにマルチスレッド対応にしたい状況も随分と増えてきました。計算時間が数分の一になる程度ですが、何時間もかかる計算もあります。GoやFortran 2008のように並列化容易なような言語もありますが、C言語など枯れた言語を使い続けたいと言うニーズもあります。この用途で人気なのはOpenMPです。

OpenMPコメントアウトを利用したプリプロセッサとして実現しています。元ソースコード中のコメントにOpenMPへの命令を書いておき、OpenMPにそれに従って元ソースコードをSMP対応ソースコードに修正してから、コンパイルして実行ファイルを得るわけです。コンパイラが対応していると、コンパイルオプションを指定するだけで利用できます。

浮動小数点演算では最速のケースも多々あるJavaは数値演算でも使われるケースが多いでしょうが、OpenMPは公式にはサポートしていません。しかし、OpenMP互換プリプロセッサomp4jと言うのがあるので試しに使ってみました。チュートリアルの通り手軽に使える反面、エラーメッセージから問題究明に慣れが必要そうでした。

使うには、まずソースコードのforループなどの前に、コメントアウトした omp で始まる命令を書きます。単に命令を追加すれば良いだけではなく、平行処理するループ中ではbreakできないなど縛りがあるので、注意してください。

// omp parallel for
for (int nProductivity = 0; nProductivity < nGridProductivity; nProductivity++) {

次に、omp4j-1.2.jarを使ってコンパイル

java -jar omp4j-1.2.jar RBC_Java.java

コンパイルができたら、実行。

java RBC_Java

ベルマン方程式によって簡単なマクロ経済モデルを解いたところ、倍以上、速くなりました。なお、言語が違うわけで当然ですが、公式と挙動が完全に一致するわけでは無いようです。例えばomp atomicでローカル変数の更新が上手くいかないので、volatileをつけたメンバー変数を更新させるなど工夫が要ると思います。

Rでテキストファイルを読み込んで分解

RでCSVファイルを上手く読み込めないときは、read.tableにあわせてテキスト・ファイルを書き換えた方が手軽な気がしますが、それが都合で出来ないときは、テキスト・ファイルとしてCSVファイルを読み込んでから、カンマで分離することもできます。

con <- file("example.csv", "r", blocking = FALSE)
lines <- readLines(con)
close(con)
lst <- strsplit(lines, ",")
row <- 1
col <- 1
sprintf("%d行%d列 %s", row, col, lst[[row]][col])
# as.numeric()で数値型に直すなど、適時型変換が必要

なお、strsplitは正規表現が使えます。
特定の列のデータをベクトルにしたい場合は以下のようにすると手軽です。

sapply(lst, "[[", 1) # 1列目を取り出す

お気軽にAR(1)の構造転換点を区間推定してみる

学術的な意義は無いのですが、時系列データの構造転換点を求めたいときがあります。記述的に転換点を仮定して推定することが多いわけですが、主観的になりやすく言い合いになることがあるからです。統計学者に殺されそうな荒業ですが、尤度を使って構造転換点を探し、さらに区間推定してみましょう。
時系列データであれば何でも良いのですが、AR(1)を考えます。

1. モデル

 y_t = C + \beta y_{t-1} + \gamma (t \gt T) D + \epsilon

 y_tt期の観測値、Cは切片項、Tは推定する構造転換点、t \gt TはtがTより大のときに1、他は0をとる構造転換ダミー、\beta\gammaは推定する係数、\epsilonは誤差項です。
なお、AR(1)にダミーを入れていいのかは謎です。

2. データセット

モデルにそって、データセットを作成します。

set.seed(20150214)
n <- 100 # 観測数
y <- numeric(n) # 観測値
C <- 3 # 切片項
gamma <- 4 # 構造転換ダミー
T <- round(n/2) # 構造転換点
beta <- 0.8 # β
epsilon <- rnorm(n) # ε
y[1] <- 15 + epsilon[1]
for(i in 2:n){
  y[i] <- C + beta * y[i-1] + epsilon[i] + (i>T)*gamma
}

# 生成データを確認する
plot(y, type="l", main="データセット", xlab="t", ylab="y")

3. 構造転換点が分かっている場合

構造転換点Tが固定で分かっているのであれば、簡単に推定できます。

d <- (2:n>50)*1 # tが50以下は0、51以上は1になる
r <- lm(y[-1] ~ y[-n] + d)
summary(r)

しかし、実際にはT=50は分かりません。40〜60のどれを仮定しても、それらしい結果になります。

4. 対数尤度を最大にする点を探す

Tを変えて説明変数の数は増えたりしないので、対数尤度を最大化する点を探しましょう。対数尤度関数を書いてニュートン・ラフソン法などで探したくなりますが、構造転換点は不連続なため上手く推定できません。ループして力技で最大値を求める方が確実です。

logLik <- numeric(n) # 対数尤度の保存用ベクトル
logLik[1] <- logLik[n] <- NA
ml_scp <- 0 # 尤度が最大の点
r <- lm(y[-1] ~ y[-n]) # 尤度が最大の推定結果
for(scp in 2:(n-1)){
  d <- (2:n>scp)*1
  tmp <- lm(y[-1] ~ y[-n] + d)
  logLik[scp] <- logLik(tmp) # 対数尤度は蓄えておく
  if(logLik(r) < logLik(tmp)){
    r <- tmp
    ml_scp <- scp
  }
}
summary(r)
sprintf("最尤推定された構造転換点: %d", ml_scp)

5. 構造転換点の区間推定を行なう

構造転換点が分かれば十分なときも多いわけですが、信頼性を疑われるので区間推定が行なえるように標準誤差を求めたいです。しかし、対数尤度の集合は離散データなので、このままでは最尤推定値の標準誤差は求まりません。

そこでスプラインで補間して、連続な対数尤度関数モドキを作ってしまいましょう。

t <- 2:(n-1)
sp <- smooth.spline(t, logLik[t])
f <- function(p){
  predict(sp, p)$y
}

プロットするとそこそこの精度で近似できていることが分かります。全体としては丸くなるので、区間推定の範囲は広くなりそうですが、狭くなるよりは良いでしょう。

plot(logLik, type="p", main="対数尤度", xlab="t", ylab="y")
lines(f(1:n))

作った対数尤度関数モドキから、構造転換点を再推定します。

r_ml <- nlm(function(p){
  -f(p)
}, ml_scp, hessian=TRUE)
SEs <- sqrt(diag(solve(r_ml$hessian))) # 標準誤差

これで構造転換点の標準誤差が出ました。区間推定もしてみましょう。

interval <- function(beta,se,range,nof){
  a <- 1 - range/2
  sprintf("%.3f(%d%%信頼区間%.3f〜%.3f)",beta,as.integer((1-range)*100),beta-se*qt(a, nof),beta+se*qt(a, nof))
}
interval(r_ml$estimate[1], SEs[1], 0.05, summary(r)$df[2])

厳密な方法ではありませんが、tは49から51ぐらいを見とけば良いとなります。

特定月の分散が大きいときの大数の法則

別のブログで使った大数の法則を示すシミュレーションのソースコードです。

# 乱数から分析データを作る
set.seed(20150209)
x.sd <- 16.96534 # 想定標準偏差
n <- 24 # 分析年数
x <- as.numeric(n*12) # 分析データ
for(i in seq(1, (n-1)*12, 12)){
  x[i:(i+11)] <- c(rnorm(3, sd=x.sd), rnorm(1, sd=x.sd*2), rnorm(8, sd=x.sd))
}

# 分析データの状態を確認
plot(x, type="l", main="集計前データ", xlab="", ylab="")

# n年分を合計してプロットする
x.sum <- numeric(12)
for(i in 1:(n-1)){
x.sum <- x.sum + abs(x[(i*12-11):(i*12)])
}
max <- (as.integer(max(x.sum)/n/5)+1)*5
B <- barplot(x.sum/n, names.arg=sprintf("%d月", 1:12), ylim = c(0, max), main=sprintf("シミュレーション(%d年分)", n))

# 偏差値を計算してみる
y.mean <- mean(x.sum/n)
y.sd <- sd(x.sum/n)
sprintf("4月の偏差値: %.2f", ((x.sum/n)[4]-y.mean)/y.sd*10+50)

Rで各集団における大きい方から2番目の値を調べる

ツイッターで見かけた御題なのですが、(学校の課題などで)ありそうなので、置いておきます。

# お試しデータを作る
set.seed(20150208)
n <- 30
df <- data.frame(type=c("A", "B", "C")[round(runif(n, min=0.5, max=3.5))], value=60 - ((1:n)-10)^2)

# タイプごとに最大値を求める
tapply(df$value, df$type, max)

# 2番目の値を求める
tapply(df$value, df$type, function(x){
  # 降順ソートをして2番目を返す
  # 最大値が複数ありえる場合は、後述の「最大値を除外する: O(n)が2回」のコードか、追記したコードを使ってください
  sort(x, decreasing=TRUE)[2]
})

奥村晴彦氏から計算量がO(n log n)でオセーヨってツッコミが来たので、O(n)になるように改良してみましたが、300万件でもuserで2割ぐらいしか時間を削れませんでした*1小選挙区の数は300しかないので(謎)、遅いコードで許してください。
以下は時間計測に使ったコードです。

set.seed(20150208)

# 300万件でテスト
n <- 3000000

# データは散らばるようにしておく
df <- data.frame(type=c("A", "B", "C")[round(runif(n, min=0.5, max=3.5))], value=runif(n, min=0, max=n))

# ソートする遅い版: O(n log n)
gc();gc();system.time({
  tapply(df$value, df$type, function(x){
   sort(x, decreasing=TRUE)[2]
  })
})

# 最大値を除外する: O(n)が2回
gc();gc();system.time({
  tapply(df$value, df$type, function(x){
    max( x[x!=max(x)] )
  })
})

# コメントで示唆されたもの: O(n)が2回
gc();gc();system.time({
  tapply(df$value, df$type, function(x){
    x[which(x==max(x))] <- NA;
    max(x, na.rm=TRUE)
  })
})

追記(2018年5月19日)

これでは同値で1位があったときに2番目に大きい数字が取れないと言う指摘があったので、(上述の最大値を除外する: O(n)が2回のコードでも良いのですが)max2関数を定義して正しく処理できるようにします。

# データセットを作る
set.seed(20180519)
df2 <- data.frame(type=c("A", "B", "C")[round(runif(n, min=0.5, max=3.5))], value=round(runif(n)*10))

# データセットのタイプ別の値をソートして並べ替える
tapply(df2$value, df2$type, sort)

#
# 2番目もしくはranking番目に大きい数字を戻す関数
# 該当がなければNA
#
max2 <- function(x, ranking=2){
  y <- sort(x, decreasing=TRUE)
  m <- ranking
  p <- 1
  while(1<m && p<length(y)){
    if(y[p]>y[p+1]){
      m <- m - 1
    }
    p <- p + 1
  }
  if(1<m){
    return(NA)
  }
  y[p]
}

# 2番目を求める
tapply(df2$value, df2$type, max2)

# 3番目を求める
tapply(df2$value, df2$type, function(x){
  max2(x, ranking=3)
})

*1:3000万件でも3割ぐらいしか削れないので、Rのsort関数が速く、max関数が遅いんじゃないかと言う疑惑が・・・