plyr なんて...（その8）

ある変数の値別に，別々のファイルにデータを書き出すにはどうするかという例

> f1 = function() d_ply(aq, .(Month), function(x) write.csv(x, paste0("airquality_", x$Month[1], ".csv")))
>
> f2 = function() invisible(sapply(split(aq, aq$Month), function(x) {
+         rownames(x) = NULL
+         write.csv(x, sprintf("airquality-%i.csv", x$Month[1]))
+         }))
> f1()
> f2()
> for (i in 5:9) {
+     d1 = read.csv(sprintf("airquality_%i.csv", i), header=TRUE)
+     d2 = read.csv(sprintf("airquality-%i.csv", i), header=TRUE)
+     cat(i, ":", identical(d1, d2), "\n")
+ }
5 : TRUE
6 : TRUE
7 : TRUE
8 : TRUE
9 : TRUE
>
> # そもそも入出力は他の処理と比べてオーダーが違う。他でどんなにうまいことやっても，出力時間にかき消される（replications=100）。引き分け。
> benchmark(f1(), f2(), columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), replications=100, order=NULL)
  test replications elapsed relative user.self sys.self
1 f1()          100   8.754    1.007     8.489    0.187
2 f2()          100   8.697    1.000     8.551    0.199

plyr なんて...（その7）

ある変数ごとに，全ての変数中に欠損値が幾つあるかを求めよという例

> f1 = function() ddply(aq, .(Month), colwise(nmissing))
> (ans = f1())
  Month Ozone Solar.R Wind Temp Day HighTemp
1     5   500     400    0    0   0        0
2     6  2100       0    0    0   0        0
3     7   500       0    0    0   0        0
4     8   500     300    0    0   0        0
5     9   100       0    0    0   0        0
>
> # 基本的には lapply
> f2 = function() {
+     a = as.data.frame(t(sapply(split(aq, aq$Month), sapply, nmissing)))
+     a$Month = as.integer(rownames(a))
+     a = a[c(5,1:4,6,7)]
+     rownames(a) = NULL
+     a
+ }
> identical(ans, f2())
[1] TRUE
>
> # by でもよいけど，遅い
> f3 = function() {
+     b = by(aq, aq$Month, colwise(nmissing))
+     a = cbind(as.integer(names(b)), do.call("rbind", b))[,-6]
+     colnames(a)[1] = "Month"
+     rownames(a) = NULL
+     a
+ }
> identical(ans, f3())
[1] TRUE
>
> # f2 は ddply に比べて 1.5 倍速
> benchmark(f1(), f2(), f3(), columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), replications=1000, order=NULL)
  test replications elapsed relative user.self sys.self
1 f1()         1000   6.547    1.099     6.411    0.201
2 f2()         1000   5.956    1.000     4.159    0.189
3 f3()         1000   7.379    1.239     7.217    0.209

plyr なんて...（その5）

2 つ以上の変数の組み合わせごとに，別の変数の平均値を求めるという例題

> aq = read.csv("aq", header=TRUE)
> aq = ddply(aq, .(Month), transform, HighTemp=ifelse(Temp-mean(Temp) > 0, 1, 0))
> f1 = function() ddply(aq1, .(Month, HighTemp), summarize, AveWind=mean(Wind))
> (ans = f1())
   Month HighTemp   AveWind
1      5        0 12.800000
2      5        1 10.518750
3      6        0 10.383333
4      6        1 10.091667
5      7        0  9.820000
6      7        1  8.118750
7      8        0 10.093750
8      8        1  7.406667
9      9        0 11.393750
10     9        1  8.792857
>
> f2 = function() data.frame(Month=rep(5:9, each=2), HighTemp=rep(0:1, 5), AveWind=c(by(aq$Wind, list(aq$HighTemp, aq$Month), mean)))
> # identical(ans, f2())
> all.equal(ans, f2())
[1] TRUE
>
> f3 = function() data.frame(Month=rep(5:9, each=2), HighTemp=rep(0:1, 5), AveWind=unname(sapply(split(aq$Wind, list(aq$HighTemp, aq$Month)), mean)))
> # identical(ans, f3())
> all.equal(ans, f3())
[1] TRUE
>
> # aggregate は列の順序が違ってくるので [c(2,1,3)] で入れ替えている（本来不要）
> f4 = function() {
+     a = aggregate(aq$Wind, list(aq$HighTemp, aq$Month), mean)[c(2,1,3)]
+     colnames(a) = c("Month", "HighTemp", "AveWind")
+     a
+ }
> identical(ans, f4())
[1] TRUE
>
> # f4 が遅いのは list のせいなので，sapply を二重にして回避
> f5 = function() {
+     a = split(aq, aq$Month)
+     b = c(sapply(a, function(d) sapply(split(d$Wind, d$HighTemp), mean)))
+     data.frame(Month=rep(5:9, each=2), HighTemp=rep(0:1, 5), AveWind=b, row.names=NULL)
+ }
> # identical(ans, f5())
> all.equal(ans, f5())
[1] TRUE
>
> # もう 1 つの回避法（sapply を for の代わりに使う）
> f6 = function() {
+     b = c(sapply(5:9, function(i) sapply(0:1, function(j) mean(aq[aq$Month==i & aq$HighTemp==j, "Wind"]))))
+     data.frame(Month=rep(5:9, each=2), HighTemp=rep(0:1, 5), AveWind=b, row.names=NULL)
+ }
> # identical(ans, f6())
> all.equal(ans, f6())
[1] TRUE
>
> # さらにもう 1 つの回避法（値の組合せを expand.grid で作る）
> f7 = function() {
+     l = expand.grid(0:1,5:9)
+     b = c(apply(l, 1, function(ij) mean(aq[aq$Month==ij[2] & aq$HighTemp==ij[1], "Wind"])))
+     data.frame(Month=l[,2], HighTemp=l[,1], AveWind=b, row.names=NULL)
+ }
> # identical(ans, f7())
> all.equal(ans, f7())
[1] TRUE
>
> # ddply は，たいして速くない。sapply を for 代わりに使うと 1.5 倍速
> # aggregate は遅い
> benchmark(f1(), f2(), f3(), f4(), f5(), f6(), f7(), columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), replications=1000, order=NULL)
  test replications elapsed relative user.self sys.self
1 f1()         1000   9.691    1.533     9.367    0.331
2 f2()         1000  15.140    2.395    14.905    0.256
3 f3()         1000  13.162    2.082    13.076    0.167
4 f4()         1000  43.494    6.881    41.944    0.580
5 f5()         1000  16.176    2.559    15.947    0.302
6 f6()         1000   6.321    1.000     5.998    0.362
7 f7()         1000   6.586    1.042     6.273    0.354

plyr なんて...（その6）

各変数ごとに欠損値が幾つあるかを求めよという例題

> nmissing = function(x) sum(is.na(x))
> # colwise は，do.call を使っているし，その本質は lapply だ（関数定義を見よ）
> f1 = function() colwise(nmissing)(aq)
> (ans = f1())
  Ozone Solar.R Wind Temp Month Day HighTemp
1  3700     700    0    0     0   0        0
>
> # 出力形式にこだわらないなら as.data.frame は不要
> f2 = function() as.data.frame(do.call("lapply", list(aq, nmissing)))
> identical(ans, f2())
[1] TRUE
>
> # 出力形式にこだわらないなら colwise よりは，速い
> f22 = function() unlist(do.call("lapply", list(aq, nmissing)))
> f22()
   Ozone  Solar.R     Wind     Temp    Month      Day HighTemp
    3700      700        0        0        0        0        0
>
> # 出力形式にこだわらないなら as.data.frame は不要
> f3 = function() as.data.frame(lapply(aq, nmissing))
> identical(ans, f3())
[1] TRUE
>
> # 出力形式にこだわらないなら colwise よりは，速い
> f32 = function() unlist(lapply(aq, nmissing))
> f32()
   Ozone  Solar.R     Wind     Temp    Month      Day HighTemp
    3700      700        0        0        0        0        0
>
> # colSums は データフレームには弱い
> # as.data.frame(t( )) で結果を無理矢理データフレームにしている
> f4 = function() as.data.frame(t(colSums(is.na(aq))))
> f4()
  Ozone Solar.R Wind Temp Month Day HighTemp
1  3700     700    0    0     0   0        0
>
> # 素のままの do.call("lapply",) や lapply の方が colwise より少しだけ速い
> benchmark(f1(), f2(), f22(), f3(), f32(), f4(), columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), replications=1000, order=NULL)
   test replications elapsed relative user.self sys.self
1  f1()         1000   0.710    1.703     0.654    0.068
2  f2()         1000   1.102    2.643     1.064    0.045
3 f22()         1000   0.426    1.022     0.380    0.049
4  f3()         1000   1.107    2.655     1.058    0.055
5 f32()         1000   0.417    1.000     0.369    0.049
6  f4()         1000   1.574    3.775     1.556    0.024

plyr なんて...（その4）

> # ifelse は使わない方がよい（Temp-mean(Temp) > 0 も Temp > mean(Temp) でよい）
> f1 = function() ddply(aq, .(Month), transform, HighTemp=ifelse(Temp-mean(Temp) > 0, 1, 0))
> ans = f1(); head(ans); tail(ans)
  Ozone Solar.R Wind Temp Month Day HighTemp
1    41     190  7.4   67     5   1        1
2    36     118  8.0   72     5   2        1
3    12     149 12.6   74     5   3        1
4    18     313 11.5   62     5   4        0
5    NA      NA 14.3   56     5   5        0
6    28      NA 14.9   66     5   6        1
      Ozone Solar.R Wind Temp Month Day HighTemp
15295    14      20 16.6   63     9  25        0
15296    30     193  6.9   70     9  26        0
15297    NA     145 13.2   77     9  27        1
15298    14     191 14.3   75     9  28        0
15299    18     131  8.0   76     9  29        0
15300    20     223 11.5   68     9  30        0
>
> # Reduce も場合によっては使い出がある。行名にこだわらなければ rownames(aq2) = NULL は不要。
> f2 = function() {
+     aq2 = Reduce(rbind, lapply(split(aq, aq$Month), function(x) {x$HighTemp = (x$Temp > mean(x$Temp))+0; x}))
+     rownames(aq2) = NULL
+     aq2
+ }
> identical(ans, f2())
[1] TRUE
>
> # ベクトル演算は強力だ
> f3 = function() {
+     aq3 = aq
+     m = c(by(aq$Temp, aq$Month, mean))
+     aq3$HighTemp = (aq$Temp > m[aq$Month-4])+0
+     aq3
+ }
> ans3 = f3()
> ans3 = Reduce(rbind, split(ans3, ans3$Month))
> rownames(ans3) = seq_len(nrow(ans3))
> identical(ans, ans3)
[1] TRUE
>
> # もっとも単純な f3 は，ddply, transform を使うものの，8 倍速だ
> benchmark(f1(), f2(), f3(), columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), replications=100, order=NULL)
  test replications elapsed relative user.self sys.self
1 f1()          100   2.986    8.967     2.924    0.092
2 f2()          100   1.742    5.231     1.655    0.100
3 f3()          100   0.333    1.000     0.314    0.022

plyr なんて...（その3）

ある変数が取る値ごとに別の変数の平均値を求める（ただし，変数には欠損値がありますよ）と言う例題

> f1 = function() ddply(aq, .(Month), summarize, AveOzone=mean(Ozone, na.rm=TRUE))
> (ans = f1())
  Month AveOzone
1     5 23.61538
2     6 29.44444
3     7 59.11538
4     8 59.96154
5     9 31.44828
>
> # sapply で使う関数に対するオプションは 第 3 引数以降に指定できる（無名関数を使う必要はない）
> f2 = function() {
+     a = split(aq$Ozone, aq$Month)
+     data.frame(Month=as.integer(names(a)), AveOzone=sapply(a, mean, na.rm=TRUE), row.names=NULL)
+ }
> identical(ans, f2())
[1] TRUE
>
> # この場合も同じく（基本は同じなので），ddply を使った f1 は，ほぼ 3 倍遅い
> benchmark(f1(), f2(), columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), replications=1000, order=NULL)
  test replications elapsed relative user.self sys.self
1 f1()         1000   5.500    3.519     5.306    0.206
2 f2()         1000   1.563    1.000     1.508    0.061

plyr なんて...（その2）

ある変数が取る値ごとに，別の変数の平均値と標準偏差を求めるという例題

> f1 = function() ddply(aq, .(Month), summarize, AveTemp=mean(Temp), SdTemp=sd(Temp))
> (ans = f1())
  Month  AveTemp   SdTemp
1     5 65.54839 6.744490
2     6 79.10000 6.488762
3     7 83.90323 4.246023
4     8 83.96774 6.479216
5     9 76.90000 8.216599
>
> # 最後の row.names=NULL は，行名を通し番号にするため
> f2 = function() {
+     a = t(sapply(split(aq$Temp, aq$Month), function(x) c(mean(x), sd(x))))
+     data.frame(Month=as.integer(rownames(a)), AveTemp=a[,1], SdTemp=a[,2], row.names=NULL)
+ }
> identical(ans, f2())
[1] TRUE
>
> # f2 と同じくらい速い
> f3 = function() {
+     a = split(aq$Temp, aq$Month)
+     data.frame(Month=as.integer(names(a)), AveTemp=sapply(a, mean), SdTemp=sapply(a, sd), row.names=NULL)
+ }
> identical(ans, f3())
[1] TRUE
>
> # ddply を使った f1 は，ほぼ 3 倍遅い
> benchmark(f1(), f2(), f3(), columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), replications=1000, order=NULL)
  test replications elapsed relative user.self sys.self
1 f1()         1000   5.824    3.274     5.634    0.225
2 f2()         1000   1.779    1.000     1.725    0.062
3 f3()         1000   1.793    1.008     1.708    0.098

plyr なんて...（その1）

plyr がもてはやされるが，いかがなものか。

「plyrパッケージ徹底入門」で紹介された用例について，検証してみよう。
http://www.slideshare.net/teramonagi/tokyo-r30-20130420

「データを分割し，分析し，結果を再結合して出力する」とはいうが，結果の出力書式は本来は意味のないものとは思う。
論文やレポートに収録するには，R のコンソール出力の丸々コピーでは済まないだろうから。

しかしまあ，一応できる限り同じ出力になるようにしてみよう。

そして，実際の処理時間を比較してみよう。

まずは，検証の環境についてまとめておこう。

> sessionInfo()
R version 3.2.2 (2015-08-14)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.1 (El Capitan)

> library(rbenchmark) # ver. 1.0.0
> library(plyr)       # ver. 1.8.3
> aq
> nrow(aq)
[1] 15300

まずは，ddply。

ある変数が取る値ごとに，別の変数の平均値を求めるという例題。

ベンチマーク・テストをするので，便宜的に引数のない関数とする。
> f1 = function() ddply(aq, "Month", summarize, AveTemp=mean(Temp))
> f2 = function() ddply(aq, .(Month), summarize, AveTemp=mean(Temp))
> f3 = function() {a=sapply(split(aq$Temp, aq$Month), mean);data.frame(Month=names(a), AveTemp=a, row.names=NULL)}

どれも，出力は同じになる
> f1()
  Month  AveTemp
1     5 65.54839
2     6 79.10000
3     7 83.90323
4     8 83.96774
5     9 76.90000
> f2()
  Month  AveTemp
1     5 65.54839
2     6 79.10000
3     7 83.90323
4     8 83.96774
5     9 76.90000
> f3()
  Month  AveTemp
1     5 65.54839
2     6 79.10000
3     7 83.90323
4     8 83.96774
5     9 76.90000

ベンチマーク・テスト
> benchmark(f1(), f2(), f3(), columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), replications=1000, order=NULL)
  test replications elapsed relative user.self sys.self
1 f1()         1000   3.998    3.483     3.704    0.346
2 f2()         1000   3.850    3.354     3.566    0.330
3 f3()         1000   1.148    1.000     1.075    0.087

考察

ddply のタイプ数はほぼ半分で済むが，処理時間は 3 倍以上かかる。

そうは言っても，1 回当たりの処理時間は 4 ミリ秒なんだから，速いの遅いの言っても意味がない。

時系列データだって言ってもねえ...(3)

最高気温について 2 日目以降に前日との差をとって，それと 2 日目以降の客数について回帰分析してみる。

データは，図のようになる。

両者の相関は 0.8001206 んんっ。前より高いし，ますます特定の相関係数を設定してデータを作ったのかと疑いの目でwww

単位根検定もパスした（いずれも P 値は 0.01）。

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  316.402      8.239    38.4  < 2e-16
x             14.181      2.579     5.5 3.91e-05

まあ，「昨日も暑かったけど，今日は更に暑いなあ」というのが「じゃあ，アイスクリームでも食べようか」という誘因になると考えて無理はなかろう。

はいはい。決定，決定。

暑い日は，アイスクリームが売れるんです。

終～～～了～～～

 

時系列データだって言ってもねえ...(2)

向後先生のアイスクリームデータ
http://kogolab.chillout.jp/elearn/icecream/chap4/sec1.html
は，土日を除く平日の 20 日分。気温と客数の相関係数は 0.8703519。

なんか，事前に相関係数=0.987 を指定して，二変数正規乱数を発生させてデータを作ったっぽいなあ。

で，PP.test（Phillips-Perron Test for Unit Roots）ってのをやってみると，

最高気温の P 値は 0.3427 単位根過程でないとはいえない，客数の P 値は 0.01 で単位根過程ではないと言う結果。

まあ，暑いさなかだと，前の日が暑ければ当日も翌日も暑いんだろうなあとは思う。それに比べ，客足は単位根過程とは思えない。

元データ（図の上側）を見ただけでは単位根過程かどうかは分かりづらいだろう。

ちなみに，20対のデータを並べ替えると，相関係数は変わらないが，単位根検定の結果は変わる。

並べ替えたデータが図の下側。両者で大違いということでもないように見えるし，サンプルサイズが 20 というのも少ない気がするが，検定結果は尊重すべきだろう。

向後先生がどうやってこのデータを作った（実際に採取した？）かは分からないが，データ対をもっと多くすればどちらに転がっていくだろうか。

作られたデータでなく実際のデータなら，最高気温は単位根過程になるんだろうか気象庁データで調べてみようか？

 

念のため。

本当は，共和分も調べなければならないのだけど，取りあえず。

時系列データだからダメなわけではない。単位根過程のデータかどうかをまずは調べる。

それらがいずれも大丈夫だとして，実質的にも両変数間に関連があるということに異論がなければ，時系列データに回帰分析したってかまわない...ん，だろ？

ブルートフォースで十分

「グー」「チョキ」「パー」の「じゃんけん」により多数決をする。
一番多くの人が出した手が勝つ。
「一度で」勝つ手が決まるような人数の組み合わせは何通りあるか。
例えば、4 人の場合は 12 通り。3 種の人数は以下の 4 通り× 3。
4 0 0
3 0 1
3 1 0
2 1 1

ブルートフォースで簡単に書けるが，出題元のコンピュータが遅いので，このままでは時間制限に引っかかる。

> f = function(n) {
+     count = 0
+     for (a in n:0) {
+         bc = n-a
+         for (b in bc:0) {
+             c = bc-b
+             if (a == max(c(a, b, c)) && a != max(c(b, c))) {
+                 count = count+1
+             }
+         }
+     }
+     3*count
+ }
> system.time(f(1234))
   ユーザ   システム       経過  
     1.301      0.009      1.299

時系列データだって言ってもねえ...

> https://twitter.com/KuboBook/status/660779858327269377/photo/1

「アイスクリームの販売数~気温」といった Y ~ X といった線形回帰・線形相関等が調べられていますが，X と Y が時系列データである場合，これは「見せかけの回帰」(spurious regression) とよばれるものになっています

まさかね。そこまで疑う？？

気温が高くなればアイスクリームが売れるというのが，見せかけの相関と？？

見せかけの相関を除いた真の相関関係はどうなっているのか，文句付けるだけで，何の分析もできないの？

時系列だから問題というなら，時系列のせいではない（時系列の縛りのない）データ，例えば，10日ごとのデータについて分析するならよいのかな？

いずれにしろ，気温が低い日よりは高い日の方がアイスクリームは売れそうだがなぁ？？

おお，向後先生だ
http://kogolab.chillout.jp/elearn/icecream/

http://www.seto.nanzan-u.ac.jp/st/gr-thesis/ms/2003/matsuda/00mm051.pdf

http://room.assign-navi.jp/column/3276/

http://xica-inc.com/blog/adelie-support-blog/2013/0828/

https://books.google.co.jp/books?isbn=4492045686

 

run length encoding

以下の文字列を出力するプログラムを書けと
bPP00$zzzuuuVVUUUQQ11000vvvWWWVVVRRR22211wwXXXWWWSSS333222xxYYYXXTTT444333yyyZZZYYUUU555444zzzaaZZZVV66555000bbaaaWW7766111cccbbbXXX8887722dddcccYYY999888333eeeddZZZAAA999444ffeeeaaaBBBAAA555ggfffbbCCBBB666hhhgggccDDCC77iiihhhdddEEEDD88jjjiieeeFFFEEE999kkkjjfffGGGFFFAAAllkkkggHHGGGBBBmmlllhhIIHHCCCnnnmmmiiiJJJIIDDooonnnjjjKKKJJJEEpppookkkLLLKKKFFFqqqpplllMMMLLLGGGrrqqqmmNNMMMHHHsssrrrnnOONNIItttsssttoooPPPOOJJuuuuupppQQQPPPKKKvvvvvvqqqRRRQQQLLLwwwwwrrrSSSRRRMMMxxxxxssTTSSSNNNyyyzzzyyytttUUUTTOOPP000zzuuuVVVUUUQQQ11100vvvWWWVVVRRR22111wwwXXXWWWSSS33222xxYYXXXTTT444333yyZZYYUU555444zzzaaaZZVVV66655000bbbaaaWWW77766111cccbbbXXX88777222dddcccYYY99988833eedddZZAAA99944ffeeaaBBBAA555gggffbbbCCCBB666hhhgggcccDDCCC777iiihhhdddEEDDD888jjjiiieeeFFFEEE99kkjjjffGGGFFFAAAlllkkggHHHGGBBBmmmlllhhhIIIHHCCCnnnmmmiiiJJIIIDDDooonnnjjjKKJJJEEppoookkkLLLKKKFFqqppllMMMLLLGGGrrrqqmmmNNNMMHHHsssrrrnnnOONNNIIItttssstttoooPPOOOJJJuuuuuupppQQQPPPKKvvvvqqRRRQQQLLwwwwrrSSSRRMMMxxxxxxsssTTTSSNNNyyyzzzyyytttUUTTTOOO

R には，"rle" Run Length Encoding という関数がある。しかし，それをそのまま使おうと思っても，短いプログラムにはならないので，あれこれこそくな手段を使う。

R による 389 文字解
g=gsub;w=c("zuVU",c(0:9,LETTERS,letters)[outer(c(26,1,0,57,32,31),0:23,"+")%%62+1],"toPOJupQPKvqRQLwrSRMxsTSNyzytUTO");cat(rep(unlist(strsplit(c("bP0$",w,"P0",w),"")),unlist(strsplit(g(7,323,g(8,33332,g(9,3883832282228,g(0,78832,"12213722338238023722238233920237222382380282223732253333633335333352233382232937883338273222277328383338232823293328338336333242332427637833")))),""))),sep="")

AWK による 472 文字解
{c="uVUQ10vWVR21wXWS32xYXT43yZYU54zaZV650baW761cbX872dcY983edZA94feaBA5gfbCB6hgcDC7ihdED8jieFE9kjfGFAlkgHGBmlhIHCnmiJIDonjKJEpokLKFqplMLGrqmNMHsrnONItstoPOJupQPKvqRQLwrSRMxsTSNyzytUTO";$0="12213722338238788322372223823392788322372223823878832282223732253333633335333352233382232937883338273222277328383338232823293328338336333242332427637833";gsub(9,3883832282228);gsub(8,33332);gsub(7,323);for(i=1;i<367;i++)for(j=0;j<substr($0,i,1);j++)printf substr("bP0$z"c"P0z"c,i,1)}