餡子付゛録゛

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

プレビュー画面で特定情報を出力しないための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関数が遅いんじゃないかと言う疑惑が・・・

時系列データの季節調整をしてみよう

時系列データには、季節バイアスが入りがちです。年中行事はもちろんのこと、天候の変化も周期的に発生します。ゴールデンウィークのある5月よりも6月の方が行楽客が少ないとしても、行楽客が減少し出したとは言えないでしょう。そもそも月ごとに日数も異なりますし。月次データを見るときは、季節バイアスの影響を考慮する必要があります。これを数字の処理で行なうのが、季節調整です。実際にRで試してみて、どのぐらい調整できるか見てみましょう。単純移動平均とloessアルゴリズムを用います。

1. 精度を確認する手順

実データでは季節バイアスがどれぐらい入っているのか真実は誰も知らないので、季節バイアスが入った擬似データを作成します。二次方程式から作った非季節変動値に、乱数から作った真の季節バイアスを加えて、観測データを作ります。この観測データを処理して計算された季節バイアスを取り出し、真の季節バイアスと就き合わせてみましょう。

2. 擬似データを作成する

1996年1月から60ヶ月の月次データとして、12ヶ月周期の擬似データを作成します。

# 結果が同じになるように乱数シードを固定
set.seed(20150119)

# 12ヶ月5年分で観測数60
n <- 12*5

# 季節バイアスを乱数から生成
s_bias <- runif(12, min=0, max=4)
# 平均0に正規化しておく
s_bias <- s_bias - mean(s_bias)

# 季節バイアスを連結して、5年分の季節変化をつくる
# tsで時系列データ型にしておく
s_chg <- ts(rep(s_bias, 5), start=c(1996, 1), frequency=12)

# 二次方程式から非季節変動値を作る
x <- 1:n
e <- 0 # 誤差はゼロとしておく。入れるときは rnorm(n) などn個のベクトルを指定。
nsv <- ts(0.9 + 0.003*(x-2*n/3)^2 + 0.1*x, start=c(1996, 1), frequency=12)

# 非季節変動値と季節バイアスを足して、原数値を作成する
y <- ts(nsv + s_chg + e, start=c(1996, 1), frequency=12)

3. 単純移動平均

12ヶ月の算術平均をとる単純移動平均は、単純な計算の割りには効果的に季節バイアスを除去します。最初の11ヶ月はデータが無い欠損値になること、非季節変動値に変化があってから効果が見えるのが遅れことが欠点になりますが、目視するには悪く無いです。

# 移動平均
# 最初11ヶ月分の値はでない。
ma <- ts(numeric(n-11), start=c(1996, 12), frequency=12) # グラフを描くために、時系列データ型にしておく
for(i in 12:n){
  ma[i-11] <- mean(y[(i-11):i])
}

4. loessアルゴリズム

Rの組み込み関数stl()で使えるloessアルゴリズムは、局所重み付け回帰関数を用いた高機能な季節調整方法です。ここでは詳しい手順は考えずに、Rに計算させましょう。なお、引数で与える変数が時系列データ型でないと動いてくれないので、注意してください。

stl.y <- stl(y, s.window="periodic")

5. グラフを描いて比較してみる

どれぐらいの精度で季節変動を補正できるのか、グラフを描いて比較してみましょう。この例では、かなりの精度で一致します。移動平均も遅れて動くわけですが、十分機能することが分かります。

plot(nsv + e, ylim=c(0, 10), main="季節調整方法の精度", lwd=2, xlab="", ylab="")
lines(y, lty=2, lwd=1)
lines(y - stl.y$time.series[,1], lty=3, col="red", lwd=2)
lines(ma, lty=4, col="blue", lwd=2)

legend("topleft", c("真の非季節変動値", "原数値", "loessアルゴリズム", "12ヶ月移動平均"), col=c("black", "black", "red", "blue"), lty=c(1, 2, 3, 4), lwd=c(2, 1, 2, 2), bg='white', box.col = 'black', bty="n")

ただし、これは非季節変動値が滑らかな二次関数だからで、ゼロにしている誤差が大きくなると当てはまりが悪くなります。季節変動の標準偏差と同じ大きさの標準偏差を持つ誤差項を入れると、トレンドの変化を掴むのに問題があるほどでは無い*1ですが、真の値と計算値の乖離は目で見て分かる程度になります。


6. X-12-ARIMA

最近はX-13-ARIMAが出てきたようですが、官公庁などでは米国商務省が開発した移動平均型季節調整法であるX-12-ARIMAを使っています。これはアプリケーションが配布されていて、Rからも扱うためのx12パッケージもあります。
Windows機の場合、C:/WinX12にX-12-ARIMAを置いたら、Rでx12パッケージをインストールしたあと、以下のように使います。

library("x12")
x12path("C:/WinX12/x12a/x12a.exe")
x12.y <- x12(y)
# 季節調整値を表示
slot(x12.y, "d11")

オプションや同時計算される系列が多いので扱いが難しいのですが、季節バイアスの変化にも対応できるし精度は高いです。何十年分の統計を抱えている官公庁以外、使い道が限られそうですが。

*1:実際、回帰分析を行なっても、真の値でy=0.003*x^2 - 0.14*x、loessアルゴリズムによる推定値でy=0.002998*x^2 - 0.1402*xとなり、ほとんど差が無い。