> 実験は二手に分けて行い、一方のグループだけＦＢの利用を禁止。１週間後に「生活にどの程度満足しているか」を１０点満点で聞いたところ、利用継続グループは実験後７・７５点と大きな変化はなかったが、禁止グループは７・５６点が８・１２点にアップした。 実験後、禁止グループは８８％が「幸福だ」、８４％が「人生を楽しんでいる」とし、継続グループの各８１％、７５％を上回った。禁止グループで「悲しい」と答えたのは２２％、「孤独だ」は１６％で、継続グループは３４％、２５％だった。
アップしたのではなく，元に戻ったんだろうね
たぶん，前のページを読んだ人は，「なぜ fread を使わないのだ」と，怒りに震えていることだろうwww
> library(data.table)
> system.time(dat1 <- read.csv("aq"))
ユーザ システム 経過
44.938 1.164 45.952
に対して
> system.time(dat2 <- fread("aq"))
Read 15300000 rows and 6 (of 6) columns from 0.266 GB file in 00:00:05
f 6) columns from 0.266 GB file in 00:00:05
ユーザ システム 経過
4.785 0.131 4.856
なので，10倍速い。AWK なんていうものも，及びが付かない（当たり前だ。AWK はインタプリタでしかも数値も内部では文字として保持しているのだから）。
key を設定しておくといろいろ便利なこともあるが，しておかなくても十分高速だ。
> system.time(setkey(dat2, Month))
ユーザ システム 経過
0.538 0.025 0.559
「データテーブル[i, j, その他いろいろな引数]」という書式で，かなりのデータ処理が出来る。
しかも，「by reference」でやるので，余計なメモリを使わない。by reference（参照渡し） は，プログラム言語をやった人は知っている人が多い。関数に引数をわたすとき，call by value と call by reference の区別がある。普通の R 関数は call by value（値渡し） で引数を受けるので，関数内部で値を変更しても，呼び出した側の値は変わらない。大元を変えるためには，関数の戻り値を代入してやらないといけない。そんなこんなで，メモリも（コピーを含む）処理時間もコストがかかる。call by reference は，メモリ番地を引き渡すので，ある番地を書き換えたら，大元が書き換わる（というか，大元しかないのだが）。
なんやかやいったけど，結局何がいいたいかといえば，「以下のようにすれば，高速処理が出来る」。
> system.time(print(dat2[, list(m=mean(Temp), s=sd(Temp)), by=Month]))
Month m s
1: 5 65.54839 6.743403
2: 6 79.10000 6.487682
3: 7 83.90323 4.245338
4: 8 83.96774 6.478172
5: 9 76.90000 8.215231
ユーザ システム 経過
0.281 0.015 0.284
また，dplyr というやつを使えば，（好き嫌いはあるが）以下のようにして，ある程度高速処理ができる。
> library(dplyr)
> system.time(dat2 %>% group_by(Month) %>% summarize(m=mean(Temp), s=sd(Temp)) %>% print)
Source: local data table [5 x 3]
Month m s
(int) (dbl) (dbl)
1 5 65.54839 6.743403
2 6 79.10000 6.487682
3 7 83.90323 4.245338
4 8 83.96774 6.478172
5 9 76.90000 8.215231
ユーザ システム 経過
0.397 0.046 0.431
しかし，そんなことしなくてもねぇ！
data.frame で split，sapply を使って，Month ごとの Temp の平均値と標準偏差を計算してみても，結構速い。2，3倍遅いと言っても，これだけの計算を1秒未満でやれているのだから，無問題。
> system.time(print(sapply(split(dat1$Temp, dat1$Month), function(x) c(m=mean(x), s=sd(x)))))
5 6 7 8 9
m 65.548387 79.100000 83.903226 83.967742 76.900000
s 6.743403 6.487682 4.245338 6.478172 8.215231
ユーザ システム 経過
0.929 0.129 1.041
お勧めは，入力だけは fread でさっさと済ませて，あとは慣れ親しんだdata.frame で処理する。
data.table は setDF で data.frame に，by reference で変換できる（as.data.frame を使う例を示している人が多いけど，as.data.frame では，by value なのでね）。
> class(dat2)
[1] "data.table" "data.frame"
> setDF(dat2)
> class(dat2)
[1] "data.frame"
data.frame に読み込むには，data.table=FALSE オプションをつければよい（これも，言及している人が少ない。オンラインヘルプを隅から隅まで読めば，目から鱗がポロポロ落ちる...ということもあるものだ）。
> system.time(dat3 <- fread("aq", data.table=FALSE))
Read 15300000 rows and 6 (of 6) columns from 0.266 GB file in 00:00:06
f 6) columns from 0.266 GB file in 00:00:06
ユーザ システム 経過
4.903 0.153 4.950
> class(dat3)
[1] "data.frame"
逆に data.frame を data.table にするには，setDT を使う。
> class(dat1)
[1] "data.frame"
> setDT(dat1)
> class(dat1)
[1] "data.table" "data.frame"
airquality を 100000 倍した CSV データ（15300000行）の Month ごとの Temp の平均値を求める
R でやるなら，
> system.time({a <-read.csv("aq"); print(by(a[,4], a[,5], mean))})
a[, 5]: 5
[1] 65.54839
------------------------------------------------------------------------------------------------------------------------
a[, 5]: 6
[1] 79.1
------------------------------------------------------------------------------------------------------------------------
a[, 5]: 7
[1] 83.90323
------------------------------------------------------------------------------------------------------------------------
a[, 5]: 8
[1] 83.96774
------------------------------------------------------------------------------------------------------------------------
a[, 5]: 9
[1] 76.9
ユーザ システム 経過
47.388 1.475 48.647
> nrow(a)
[1] 15300000
と，47 秒もかかる（入力時間がほとんど）
> system.time({print(by(a[,4], a[,5], mean))})
a[, 5]: 5
[1] 65.54839
------------------------------------------------------------------------------------------------------------------------
a[, 5]: 6
[1] 79.1
------------------------------------------------------------------------------------------------------------------------
a[, 5]: 7
[1] 83.90323
------------------------------------------------------------------------------------------------------------------------
a[, 5]: 8
[1] 83.96774
------------------------------------------------------------------------------------------------------------------------
a[, 5]: 9
[1] 76.9
ユーザ システム 経過
1.583 0.361 1.923
つまり，47 秒の 97% が入力に必要な時間。
AWK で書いたら，
$ cat a.awk
BEGIN { FS = "," }
NR > 1 { x[$5] += $4; n[$5]++ }
END {
for (i in x) print i, x[i]/n[i]
}
のように，簡単なスクリプトで，
$ time gawk -f a.awk aq
5 65.5484
6 79.1
7 83.9032
8 83.9677
9 76.9
17.274u 0.141s 0:17.43 99.8% 0+0k 0+0io 3pf+0w
と，ほぼ 3 倍速である。
【コードゴルフ】シンプル・ライフゲーム
シンプルなライフゲームです。できるだけ短いコードを書いてください。
締切:11月19日（木）AM10:00
【ライフゲーム】
ライフゲームでは、格子状のフィールドの各マス目（セル）に対して、「生」と「死」の2つの状態を初期値として与え、次のルールによって世代変化します。
各セルについて、
「生」の状態の場合
・周囲の8セルのうち、「生」の状態のセルが2つまたは3つ存在する場合は、「生」の状態を継続する
・周囲の8セルのうち、「生」の状態のセルが1つ以下または4つ以上の場合は、「死」の状態になる
「死」の状態の場合
・周囲の8セルのうち、「生」の状態のセルが丁度3つ存在する場合は、「生」の状態になる
☆端のセルの扱い
フィールドの上下左右端のセルはループします。
つまり、左端のセルの左は右端のセル、上端のセルの上は下端のセルと考えてください。
[問題]
フィールドの初期状態と世代数Nを標準入力から受け取り、初期状態からN世代後の状態を出力してください。
なるべく短いプログラムを書けということなので，プログラム言語によってハンデがあるだろう。
AWK で書くと，以下のようなプログラム（可能な限り空白類は削除する）
function m(i,j){return(j==1?a[i,W]:a[i,j-1])+a[i,j]+(j==W?a[i,1]:a[i,j+1])}
function k(i,j){return(i==1?m(H,j):m(i-1,j))+m(i,j)+(i==H?m(1,j):m(i+1,j))}
BEGIN{getline N;getline H;getline W;for(i=1;i<=H;i++){getline s;for(j=1;j<=W;j++)a[i,j]=substr(s,j,1)=="*"}
for(l=1;l<=N;l++){for(i=1;i<=H;i++)for(j=1;j<=W;j++)b[i,j]=a[i,j]
for(i=1;i<=H;i++){for(j=1;j<=W;j++){c=k(i,j)-a[i,j];if(!a[i,j]&&c==3)b[i,j]=1
else if(a[i,j]&&(c<2||c>3))b[i,j]=0}}for(i=1;i<=H;i++)for(j=1;j<=W;j++)a[i,j]=b[i,j]}
for(i=1;i<=H;i++){for(j=1;j<=W;j++)printf a[i,j]==1?"*":".";print""}}
R で書くと以下のようになる（空白などは残す）
con = file("stdin", "r")
N = as.numeric(readLines(con, 1))
H = as.numeric(readLines(con, 1))
W = as.numeric(readLines(con, 1))
a = readLines(con)
a = matrix(unlist(strsplit(a, ""))=="*", ncol=nchar(a[1]), byrow=TRUE)
m = function(i, j) {
ifelse(j == 1, a[i, W], a[i, j-1])+a[i, j]+ifelse(j == W, a[i, 1], a[i, j+1])
}
k = function(i, j) {
ifelse(i == 1, m(H, j), m(i-1, j))+m(i, j)+ifelse(i == H, m(1, j), m(i+1, j))
}
for (l in seq_len(N)) {
b = a
for (i in 1:H) {
for (j in 1:W) {
c = k(i, j)-a[i, j]
if (!a[i, j] && c == 3) {
b[i, j] = 1
} else if (a[i, j] && (c < 2 || c > 3)) {
b[i, j] = 0
}
}
}
a = b
}
for (i in 1:H) {
cat(paste(gsub("0", "\\.", gsub("1", "\\*", paste(a[i, ], collapse=""))), "\n", sep=""))
}
締切：11月13日（金）AM10:00まで
「セカンド・イクエイション」問題
ちなみに，二次方程式は an equation of the second degree。短くいえば，quadratic equation。
second equation は The second equation is a polynomial equation or a rational equation. のように使う（らしい。知らんけど）
3つの自然数 a, b, c に対し、これらを係数に持つ次の二次方程式を考えます。
a x2 ＋ b x ＋ c ＝ 0 … (※)
例えば、(a, b, c) = (1, 6, 4) のとき、方程式(※)は、x2 ＋ 6x ＋ 4 = 0 となります。
この方程式は実数解を持ちます。x = －3＋√5 と x = －3－√5 の二つです。
一方で、(a, b, c) = (1, 6, 10) のとき、方程式(※)は実数解を持ちません。
b ≦ 3 の範囲で、方程式(※)が少なくとも一つの実数解を持つような自然数の組 (a, b, c) は、全部で 4 組あります。
(a, b, c) = (1, 2, 1), (1, 3, 1), (1, 3, 2), (2, 3, 1) です。
自然数 m に対し、b ≦ m の範囲で方程式(※)が少なくとも一つの実数解を持つような自然数の組 (a, b, c) の個数を F(m) と定義します。
例えば F(3) = 4、F(4) = 12、F(10) = 287、F(100) = 620174 です。
標準入力から、自然数 m（1 ≦ m ≦ 3000）が与えられます。
標準出力に F(m) の値を出力するプログラムを書いてください。
短いよ。驚くな。先頭 2 行と終わりの 1 行は入出力のため。残り 4 行でプログラム。
ある数列の和になるのだ。公式があるwww
con=file("stdin", "r")
b=as.integer(readLines(con, 1))
n=0
for (i in (2:b)^2 %/% 4) {
n = n+2*sum(i %/% 1:as.integer(sqrt(i)))-as.integer(sqrt(i))^2
}
cat(n)
解けた。しかし，締切が 12月04日（金）AM10:00 なので，解答公開はその後で...
「ツイン・トライアングル」問題
数学の問題をプログラミングで解こう！三角形と外接円に関する問題です。
2つの自然数 a, b（a ≦ b）に対し、直角をはさむ2辺の長さが a と b である直角三角形を考えます。
この直角三角形の外接円を描きます。
さらに、その円に外接し、もとの直角三角形に相似となる直角三角形を描きます。
この直角三角形の直角をはさむ2辺の長さを a’, b’（a’ ≦ b’）とおきます。
例えば、(a, b) = (6, 8) のとき、(a’, b’) = (15, 20) となることが確かめられます。
整数 L に対し、a ≦ b ≦ L の範囲で、a’ と b’ がいずれも整数となるような整数の組 (a, b) の個数を F(L) と定義します。
例えば F(10) = 1、F(50) = 7、F(1000) = 173 となることが確かめられます。
標準入力から、自然数 L（1 ≦ L ≦ 105）が与えられます。
標準出力に F(L) の値を出力するプログラムを書いてください。
他の統計解析ソフトにはない，画期的で，すごく，すごく，便利な機能www
このようなデータの場合には，折れ線グラフを選ばないと悲惨な結果になる（特に今までのバージョンでは）
使いたいなら使えば...お好きにどうぞ。
それを使ったからといって，スペシャリストのお墨付きが付くわけでもない...結果次第だねえ。
基礎ができていないと，宝の持ち腐れかも...
基礎があれば，そんなモノ使わなくても，立派に生きていける。
ここに書いたもののうち，最速なものを関数化すれば，plyr と同じインターフェースでより高速な関数パッケージを書けるということになるのかな？？
私は，全くやる気がないが。
今更なんだけど，ある変数が取る値ごとのサンプルサイズ
> f1 = function() count(aq, .(Month))
> (ans = f1())
Month freq
1 5 3100
2 6 3000
3 7 3100
4 8 3100
5 9 3000
>
> f2 = function() {
+ a = ddply(aq, "Month", nrow)
+ colnames(a)[2] = "freq"
+ a
+ }
> identical(ans, f2())
[1] TRUE
>
> f3 = function() {
+ a = sapply(split(aq, aq$Month), nrow)
+ data.frame(Month=as.integer(names(a)), freq=a, row.names=NULL)
+ }
> identical(ans, f3())
[1] TRUE
>
> f4 = function() {
+ a = table(aq$Month)
+ data.frame(Month=as.integer(names(a)), freq=c(a), row.names=NULL)
+ }
> identical(ans, f4())
[1] TRUE
>
> f5 = function() {
+ b = 5:9
+ a = sapply(b, function(i) sum(aq$Month==i))
+ data.frame(Month=b, freq=a)
+ }
> identical(ans, f5())
[1] TRUE
>
> # sapply を for 代わりに使った f5 は，count に比べて 2 倍速
> benchmark(f1(), f2(), f3(), f4(), f5(), 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 1.853 2.230 1.686 0.183
2 f2() 1000 4.336 5.218 4.039 0.317
3 f3() 1000 3.919 4.716 3.632 0.334
4 f4() 1000 4.861 5.850 4.698 0.184
5 f5() 1000 0.831 1.000 0.769 0.066
複数キー（一部に逆順）のソート
> f1 = function() arrange(aq, Ozone, desc(Solar.R))
>
> # 複数個指定したり，その一部が逆順の場合も，同じように対応する
> f2 = function() {
+ d = aq[order(aq$Ozone, -aq$Solar.R),]
+ rownames(d) = NULL
+ d
+ }
> identical(f1(), f2())
[1] TRUE
>
> # ほぼ互角
> 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 9.144 1.015 8.970 0.200
2 f2() 1000 9.010 1.000 8.867 0.188
逆順ソート
> f1 = function() arrange(aq, desc(Ozone))
>
> # 逆順は decreasing=TRUE
> f2 = function() {
+ d = aq[order(aq$Ozone, decreasing=TRUE),]
+ rownames(d) = NULL
+ d
+ }
>
> # または，数値変数なら符号付け替え（desc がやっていること）
> f3 = function() {
+ e = aq[order(-aq$Ozone),]
+ rownames(e) = NULL
+ e
+ }
> a = f1()
> identical(a, f2())
[1] TRUE
> identical(a, f3())
[1] TRUE
>
> # order を使う方が若干速い
> 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.340 1.011 5.174 0.201
2 f2() 1000 5.282 1.000 5.103 0.223
3 f3() 1000 5.308 1.005 5.079 0.261
単なるソートなんだけど？
> f1 = function() arrange(aq, Ozone)
> a = f1(); head(a); tail(a)
Ozone Solar.R Wind Temp Month Day HighTemp
1 1 8 9.7 59 5 21 0
2 1 8 9.7 59 5 21 0
3 1 8 9.7 59 5 21 0
4 1 8 9.7 59 5 21 0
5 1 8 9.7 59 5 21 0
6 1 8 9.7 59 5 21 0
Ozone Solar.R Wind Temp Month Day HighTemp
15295 NA 145 13.2 77 9 27 1
15296 NA 145 13.2 77 9 27 1
15297 NA 145 13.2 77 9 27 1
15298 NA 145 13.2 77 9 27 1
15299 NA 145 13.2 77 9 27 1
15300 NA 145 13.2 77 9 27 1
>
> # order を使う
> # rownames(d) = NULL は，行名を無視するため（arrange も，最後にやっている）
> f2 = function() {
+ d = aq[order(aq$Ozone),]
+ rownames(d) = NULL
+ d
+ }
> identical(f1(), f2())
[1] TRUE
>
> # 両者互角（引き分け）
> 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.322 1.000 5.193 0.160
2 f2() 1000 5.341 1.004 5.160 0.214
複数のパラメータを持つ分析を記録するという例
> # mdply は，基本的には mapply
> f1 = function() mdply(data.frame(mean=1:3, sd=1:3), rnorm, n=2)
> set.seed(1); (ans = f1())
mean sd V1 V2
1 1 1 0.3735462 1.1836433
2 2 2 0.3287428 5.1905616
3 3 3 3.9885233 0.5385948
>
> # data.frame は本来不要
> f2 = function() {
+ mean=1:3
+ sd=1:3
+ d = data.frame(mean, sd, t(mapply(function(x, y) rnorm(2, x, y), mean, sd)))
+ colnames(d)[3:4] = paste("V", 1:2, sep="")
+ d
+ }
> set.seed(1); identical(ans, f2())
[1] TRUE
>
> # mapply は mdply の 5 倍速
> 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 2.118 4.555 2.122 0.012
2 f2() 1000 0.465 1.000 0.463 0.003
ある分析を何回か繰り返し，その結果を保存しておくという例
> f = function() {lm(Ozone ~ Temp, aq[sample(nrow(aq), 100, replace=TRUE),])}
> f1 = function() rlply(100, f)
> set.seed(1); ans1 = f1(); head(ans1, 2)
[[1]]
Call:
lm(formula = Ozone ~ Temp, data = aq[sample(nrow(aq), 100, replace = TRUE),
])
Coefficients:
(Intercept) Temp
-153.799 2.559
[[2]]
Call:
lm(formula = Ozone ~ Temp, data = aq[sample(nrow(aq), 100, replace = TRUE),
])
Coefficients:
(Intercept) Temp
-147.962 2.418
>
> f2 = function() lapply(1:100, function(i) f())
> set.seed(1); ans2 = f2(); head(ans2, 2)
[[1]]
Call:
lm(formula = Ozone ~ Temp, data = aq[sample(nrow(aq), 100, replace = TRUE),
])
Coefficients:
(Intercept) Temp
-153.799 2.559
[[2]]
Call:
lm(formula = Ozone ~ Temp, data = aq[sample(nrow(aq), 100, replace = TRUE),
])
Coefficients:
(Intercept) Temp
-147.962 2.418
>
> # lm の結果が全部がほしいわけではないというならば以下のように
> # rlply の定義の主要部は replicate とほとんど同じ
> f3 = function() replicate(100, f()$coefficients[2])
> set.seed(1); ans3 = f3(); head(ans3, 2)
Temp Temp
2.559170 2.418056
>
> # 結果的にも，引き分け（replicate が若干速い？）
> 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 15.485 1.063 15.488 0.120
2 f2() 100 15.295 1.050 15.298 0.090
3 f3() 100 14.572 1.000 14.582 0.078
