裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

どれだけ大きいデータなんだろう?

2015年07月24日 | ブログラミング

> 久保拓弥@KuboBook

> R の cor.test() で Kendall の tau を評価してる最中なんだけど,同じ cor.test() の "pearson" に比べて数百-数千倍の計算時間が必要とされてる模様…4-5 分を費やして計算終了.とうぜんながら無相関仮説が棄却されました,とゆー結果に.

> posted at 13:51:39

 

久保先生は検定なんかやらない人だと思っていたんだけどなぁ

Kendall の順位相関係数は,計算アルゴリズムからして計算時間がかかるのは当然といえば当然で

Spearman の順位相関係数でなくて Kendall の順位相関係数でなければならなかった理由はあったりするんだろ~かとか

コメント

私も相当な新しもの好きではあるんだが...

2015年07月22日 | ブログラミング

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 なんか使うか(ぺっ!)

コメント

2変数が本質的に相関のある場合とない場合の時系列データの解析

2015年07月09日 | 統計学

時系列データの解析についてのコメントが流行っているのか...

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 が説明しているので,それで十分)。

 

コメント

時間を調整した偏相関係数をみると...

2015年07月07日 | 統計学

中澤さんの 「誤解を生むグラフ」に示されているデータ

> 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)


コメント

両端揃え

2015年07月02日 | ブログラミング

昔の英文ワープロ(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)
}

コメント