> 久保拓弥@KuboBook
> R の cor.test() で Kendall の tau を評価してる最中なんだけど，同じ cor.test() の "pearson" に比べて数百-数千倍の計算時間が必要とされてる模様…4-5 分を費やして計算終了．とうぜんながら無相関仮説が棄却されました，とゆー結果に．
> posted at 13:51:39
久保先生は検定なんかやらない人だと思っていたんだけどなぁ
Kendall の順位相関係数は，計算アルゴリズムからして計算時間がかかるのは当然といえば当然で
Spearman の順位相関係数でなくて Kendall の順位相関係数でなければならなかった理由はあったりするんだろ～かとか
http://t.co/ld3MR7r2yi
「for を捨てよ、foreach を書こう」だが，どういう場合に（何をどうやったら）効果的な時間短縮が得られるのか，実例を述べて欲しいなぁ
1:n までの平方根を求める
> n <- 10000
foreach を使う
> system.time(a <- foreach(i = 1:n, .combine=c) %do% sqrt(i))
ユーザ システム 経過
2.708 0.014 2.706
単なる for を使う
> system.time({b <- numeric(n);for(i in 1:n) b[i] <- sqrt(i)})
ユーザ システム 経過
0.019 0.001 0.043
ベクトル演算を使う
> system.time(c <- sqrt(1:n))
ユーザ システム 経過
0.000 0.000 0.001
同じ結果が得られることを確認しておこう
> all.equal(a, b)
[1] TRUE
> all.equal(a, c)
[1] TRUE
だれが foreach なんか使うか（ぺっ!）
時系列データの解析についてのコメントが流行っているのか...
「R: 時系列データ間の関係を状態空間モデルでみる」
「R: 時系列データ間の関係を状態空間モデルでみる（2）」 であるが...
状態空間モデルというのは，不勉強で知らなかったのだが，偏相関係数でも対応できるように思えたのでやってみた。
z を制御したときの x と y の間の偏相関係数を計算する関数
par.cor = function(x, y, z) {
xy = cor(x, y)
xz = cor(x, z)
yz = cor(y, z)
(xy-xz*yz)/sqrt((1-xz^2)*(1-yz^2))
}
第 1 部
本質的に無関係な x, y と，時間変数 z
> set.seed(1234)
>
> N.t = 60
> mu.x = numeric(N.t)
> x = numeric(N.t)
> mu.y = numeric(N.t)
> y = numeric(N.t)
> s1 = 2
> s2 = 1
>
> mu.x[1] = mu.y[1] = 10
>
> for (t in 1:(N.t - 1)) {
+ x[t] = rnorm(1, mu.x[t], s1)
+ mu.x[t + 1] = rnorm(1, mu.x[t] + 1, s2)
+ y[t] = rnorm(1, mu.y[t], s1)
+ mu.y[t + 1] = rnorm(1, mu.y[t] + 1, s2)
+ }
> x[N.t] = rnorm(1, mu.x[N.t - 1], s1)
> y[N.t] = rnorm(1, mu.y[N.t - 1], s1)
>
> z = 1:N.t # 時間変数
分析してみる。
> resid.x = lm(x ~ z)$residuals
> resid.y = lm(y ~ z)$residuals
> plot(resid.y ~ resid.x)
> text(-6, 6, sprintf("par.cor. = %.3f", par.cor(x, y, z)), pos=4)
時間を制御したときの x, y 間の偏相関係数は小さいことがわかる。
> summary(lm(y ~ x+z))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.5894 0.8172 5.616 6.06e-07
x 0.1275 0.1117 1.141 0.259
z 0.8144 0.1175 6.934 4.10e-09
重回帰分析でも，x は y を説明しているとは言いがたい。z によって説明される。
===============================================
第 2 部
本質的に関連のある x, y と，時間変数 z は前のまま（元記事の変数名 y を y2 に変更）
> # y2 は x から生成
> y2 = rnorm(N.t, x * 2, 8)
同じように分析してみる。
> resid.x = lm(x ~ z)$residuals
> resid.y2 = lm(y2 ~ z)$residuals
> plot(resid.y2 ~ resid.x)
> text(-6, 15, sprintf("par.cor. = %.3f", par.cor(x, y2, z)), pos=4)
時間を制御したときの x, y2 間の偏相関係数は大きい（x と y2 には本質的な相関関係がある）。
> summary(lm(y2 ~ x + z))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.09809 2.24666 1.379 0.173
x 1.90342 0.30705 6.199 6.77e-08
z 0.04705 0.32293 0.146 0.885
重回帰分析でも，x は y2 を有意に説明している。
それに反して，時間変数 z は y2 を説明していない（x が説明しているので，それで十分）。
中澤さんの 「誤解を生むグラフ」に示されているデータ
> cor(x)
YEAR CUCI BEEFCC BEEFSP
YEAR 1.0000000 0.9402893 0.7666001 0.7501338
CUCI 0.9402893 1.0000000 0.5671291 0.5509674
BEEFCC 0.7666001 0.5671291 1.0000000 0.9971032
BEEFSP 0.7501338 0.5509674 0.9971032 1.0000000
CUCI と BEEFCC, BEEFSP には，ほどほどの正の相関（0.567, 0.551）
しかし，いずれも YEAR とは，高度の正相関
YEAR を制御した（YEAR の影響を除いた）偏相関係数を求めてみると，
> # partial correlation r(CUCI, BEEFCC . YEAY)
> (cor(x$CUCI, x$BEEFCC) - cor(x$CUCI, x$YEAR) * cor(x$BEEFCC, x$YEAR)) /
+ sqrt((1-cor(x$CUCI, x$YEAR)^2)*(1-cor(x$BEEFCC, x$YEAR)^2))
[1] -0.7032116
> # partial correlation r(CUCI, BEEFSP . YEAY)
> (cor(x$CUCI, x$BEEFSP) - cor(x$CUCI, x$YEAR) * cor(x$BEEFSP, x$YEAR)) /
+ sqrt((1-cor(x$CUCI, x$YEAR)^2)*(1-cor(x$BEEFSP, x$YEAR)^2))
[1] -0.6858511
いずれも，やや強い負の相関になった。
YEAR で説明できない部分の散布図を描いておく。左上の 1 点を除いても，負の相関だ。
cuci.year = lm(CUCI ~ YEAR, data=x)$residuals
beefcc.year = lm(BEEFCC ~ YEAR, data=x)$residuals
plot(cuci.year ~ beefcc.year)
text(-1, 3000, sprintf("r = %.3f (partial cor. coef.)",
cor(cuci.year, beefcc.year)), pos=4)
cuci.year = lm(CUCI ~ YEAR, data=x)$residuals
beefsp.year = lm(BEEFSP ~ YEAR, data=x)$residuals
plot(cuci.year ~ beefsp.year)
text(-0.5, 3000, sprintf("r = %.3f (partial cor. coef.)",
cor(cuci.year, beefsp.year)), pos=4)
昔の英文ワープロ（monospace font で，単語間に空白を挿入して右端を揃える）
昔，FORTRAN で実機を使わずにプログラムしたことがあったなぁ～
func = function(s, w) {
func2 = function(s) {
n = length(s)
len = sum(sapply(s, nchar))
spaces = w - len
left.part.spaces = spaces%/%(n - 1)
left.part.words = n - 1 - spaces%%(n - 1)
if (left.part.words == n - 1) {
left.part.words = left.part.words + 1
}
sp = paste(rep(" ", left.part.spaces+1), collapse = "")
s2 = paste(s[1:left.part.words], collapse = substr(sp, 2, nchar(sp)))
if (n > left.part.words) {
s3 = paste(s[(left.part.words + 1):n], collapse = sp)
s2 = paste(s2, s3, sep = substr(sp, 2, nchar(sp)))
}
s2
}
s = unlist(strsplit(s, " "))
repeat {
len = 0
n.words = length(s)
for (i in 1:n.words) {
len = len + nchar(s[i]) + 1
if (len == w + 1) {
cat(sprintf("%s\n", paste(s[1:i], collapse = " ")))
s = s[-(1:i)]
break
} else if (len > w + 1) {
cat(sprintf("%s\n", func2(s[1:(i - 1)])))
s = s[-(1:(i - 1))]
break
}
}
if (i == n.words) {
cat(sprintf("%s\n", paste(s, collapse = " ")))
break
}
}
}
###### テスト
s = "One child, one teacher, one book and one pen can change the world. Education is the only solution. Education first."
for (w in 20:30) {
func(s, w)
}
