裏 RjpWiki

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

フェイスブック利用者

2015年11月29日 | 統計学

フェイスブックやめると満足度アップ

> 実験は二手に分けて行い、一方のグループだけFBの利用を禁止。1週間後に「生活にどの程度満足しているか」を10点満点で聞いたところ、利用継続グループは実験後7・75点と大きな変化はなかったが、禁止グループは7・56点が8・12点にアップした。 実験後、禁止グループは88%が「幸福だ」、84%が「人生を楽しんでいる」とし、継続グループの各81%、75%を上回った。禁止グループで「悲しい」と答えたのは22%、「孤独だ」は16%で、継続グループは34%、25%だった。

アップしたのではなく,元に戻ったんだろうね

コメント

fread はお勧め,それ以外は...

2015年11月26日 | ブログラミング

たぶん,前のページを読んだ人は,「なぜ 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"

コメント

ちょっとしたことをやるなら R より AWK

2015年11月25日 | ブログラミング

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 倍速である。

コメント

シンプル・ライフゲーム

2015年11月19日 | ブログラミング

【コードゴルフ】シンプル・ライフゲーム

シンプルなライフゲームです。できるだけ短いコードを書いてください。

締切: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=""))
}

コメント

実数解を持つ二次方程式の係数

2015年11月17日 | ブログラミング

締切: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)

コメント

「ツイン・トライアングル」問題

2015年11月16日 | ブログラミング

解けた。しかし,締切が 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

2015年11月15日 | ブログラミング

他の統計解析ソフトにはない,画期的で,すごく,すごく,便利な機能www

 

コメント

Excel の散布図と折れ線ブラフ

2015年11月15日 | ブログラミング

このようなデータの場合には,折れ線グラフを選ばないと悲惨な結果になる(特に今までのバージョンでは)

コメント

plyr 結論

2015年11月12日 | ブログラミング

使いたいなら使えば...お好きにどうぞ。

それを使ったからといって,スペシャリストのお墨付きが付くわけでもない...結果次第だねえ。

基礎ができていないと,宝の持ち腐れかも...

基礎があれば,そんなモノ使わなくても,立派に生きていける。

ここに書いたもののうち,最速なものを関数化すれば,plyr と同じインターフェースでより高速な関数パッケージを書けるということになるのかな??

私は,全くやる気がないが。

コメント

plyr なんて...(その14)

2015年11月12日 | ブログラミング

今更なんだけど,ある変数が取る値ごとのサンプルサイズ

> 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

コメント

plyr なんて...(その13)

2015年11月12日 | ブログラミング

複数キー(一部に逆順)のソート

> 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

コメント

plyr なんて...(その12)

2015年11月12日 | ブログラミング

逆順ソート

> 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

コメント

plyr なんて...(その11)

2015年11月12日 | ブログラミング

単なるソートなんだけど?

> 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

コメント

plyr なんて...(その10)

2015年11月12日 | ブログラミング

複数のパラメータを持つ分析を記録するという例

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

コメント

plyr なんて...(その9)

2015年11月12日 | ブログラミング

ある分析を何回か繰り返し,その結果を保存しておくという例

> 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

コメント