ダメ出し：R で学ぶデータ・プログラミング入門 ―RStudioを活用する― その1

まだ，ご本家の Web サイトに正誤表などは掲載されていないようなので(2012/12/05現在）

● 6 ページ 　下から 3 行目
表示されたファイルですが　→ 表示されたウインドウですが
　開かれたのはエディタ・ウィンドウでしょう。それを保存したらファイルです。

● 17 ページ　上から 6 行目
ホィールキー　→　ホイール
　あれはホィールと呼ぶのでは？

● 32 ページ　9 行目より
この二つの引数の後ろに NULL とあるのは，この引数は省略しても構わないことを意味しています
　そうではなく，デフォルトの値が NULL になるということであって，引数の値が NULL であるときの動作（意味）はプログラムに書かれているわけです

● 34 ページ　上から 3 行目
パッケージに定義　→　パッケージで定義

● 34 ページ　コラムの 3 行目より
関数の名前の先頭に半角のクエスチョンマーク(?)を付け，関数から丸括弧を取り去った
→ 関数の名前の先頭にクエスチョンマーク(?)をつけ，
　丸括弧は関数名には含まれない。
　ちなみに，本書では関数を示すときに「seq() 関数」のように記述しているが冗長である。「seq 関数」か「seq()」でよいだろう。

● 36 ページ　コラムの 3 行目
replicate　→　replace

● 41 ページ　上から 2 行目
TRUE　→　FALSE

● 42 ページ　下から 2 行目より
左端にあるのは行番号で，これは出力の際に加えられる情報であり，データそのものには含まれていません。
　これは不正確。左端にあるのは行名。出力の際に加えられるものではなく，データそのものに含まれている。そのままの状態で dput(sleep) で書き出したものと，rownames(sleep) <- letters[1:10]; dput(sleep) によって書き出したものとを比較してみれば分かる。

● 47 ページ　上から 3 行目
ヨーロッパ人の髪と目の対応　→　アメリカ人学生の髪と目の色の対応
　survey of students at the University of Delaware　デラウエアはアメリカ合衆国ニューアーク州デラウエア

● 51 ページ　下から 9 行目
if() 関数　→　if 文
　R では if も関数ではあるが，if を 関数として使うときは x <- 3; "if"(x %% 2, print("odd"), print("even")) のような使い方になる。すなわち，例に挙げているのは if 文である。

● 56 ページ　上から 8 行目
for() 関数　→　for 文
　if の場合と同じ理由。つまり，for 関数と for 文は，書き方が違う。

● 57 ページ　上から 1 行目
while() 関数　→　while 文
　if の場合と同じ理由。

● 58 ページ　上から 8 行目
画面を改行する　→　画面上の出力行を改行する
　画面は改行できない。

● 61 ページ　上から 9 行目
行なわれ　→　行われ

関数名と，関数名に続く ( の間に空白を置くのは，気持ち悪い（for, if, while は例外...関数じゃないから）

) と { の間に空白を置かないのは，気持ち悪い

本日は取りあえずこれまで。

ダメ出し：どちらが正しいかを示すには

不偏分散と1/(n-1) では

単純に，どちらがもっともらしいか示すだけで十分では？

極端な場合を見せれば十分かな？

> n <- 2
> loop <- 10000
> u <- v <- numeric(n)
> for (i in 1:loop) {
+ x <- rnorm(n)
+ u[i] <- var(x)
+ v[i] <- var(x)*(n-1)/n
+ }
> layout(matrix(1:2, 2))
> hist(u)
> hist(v)
> layout(1)
> mean(u)
[1] 1.002552
> mean(v)
[1] 0.5012759

もうちょっとあり得そうな例

> n <- 20
> loop <- 10000
> u <- v <- numeric(n)
> for (i in 1:loop) {
+ x <- rnorm(n)
+ u[i] <- var(x)
+ v[i] <- var(x)*(n-1)/n
+ }
> layout(matrix(1:2, 2))
> hist(u)
> hist(v)
> layout(1)
> mean(u)
[1] 1.004137
> mean(v)
[1] 0.9539302

ダメ出し： カイ二乗統計量の元になる分割表データ

分散の偏り で

前にも指摘したことがあるが，カイ二乗統計量を計算する元の分割表は，カウントデータの集計表。

使われたデータのような実数値ではない。

> tmp.m
                 [,1]             [,2]
[1,] 351.515921600163 350.484078399837
[2,]  75.484078399837 321.515921600163

しかも，期待値 m.e が間違い。転置しないと正しい答にならない。

> addmargins(tmp.m)
                 [,1]             [,2] [,3]
[1,] 351.515921600163 350.484078399837  702
[2,]  75.484078399837 321.515921600163  397
[3,] 427.000000000000 672.000000000000 1099

> addmargins(m.e)
                 [,1]             [,2] [,3]
[1,] 272.751592356688 154.248407643312  427
[2,] 429.248407643312 242.751592356688  672
[3,] 702.000000000000 397.000000000000 1099

期待値を求めるプログラムも汚い。

> m <- matrix(v,2,2)
> a <- apply(m,1,sum)
> b <- apply(m,2,sum)
> m.e <- matrix(c(a[1]*b[1],a[1]*b[2],a[2]*b[1],a[2]*b[2]),2,2)/sum(v)

こんな風にしてはいかが?

> m
     [,1] [,2]
[1,]  276  426
[2,]  151  246

> (m.e <- outer(rowSums(m), colSums(m)) / sum(m)) # 簡単でしょ?
                 [,1]             [,2]
[1,] 272.751592356688 429.248407643312
[2,] 154.248407643312 242.751592356688

ダメ出し：データに NA が含まれるときには注意

Regression Analysis on R　kosugitti

標準化偏回帰係数を求めるためにもとのデータを変数ごとに標準化したものを使って重回帰分析をするという手順を示している。

reg2 <- lm(kokugo ~ sansuu + eigo, data = sample)
summary(reg2)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  22.0649     6.7036    3.29   0.0014 **
sansuu       -0.0789     0.0872   -0.91   0.3675    
eigo          0.8036     0.0340   23.64   <2e-16 ***

sample.z <- scale(sample[4:10])
sample.z <- data.frame(sample.z)
reg3 <- lm(kokugo ~ sansuu + eigo, data = sample.z)
summary(reg3)

しかし，出てきた結果がおかしい。

             Estimate Std. Error t value Pr(>|t|)    
 (Intercept)  0.00409    0.03882    0.11     0.92    
 sansuu      -0.03618    0.03996   -0.91     0.37    
 eigo         0.92324    0.03905   23.64   <2e-16 ***

どこがおかしいか分からない人もいるとは思うが，変数を標準化して重回帰分析をすると，(Intercept) は 0 にならなければならないのだ。

なぜこうなったかという理由は，元のデータ sample に欠損値があるためである。
第 1 には，scale 関数は，変数ごとに NA を自動的に除いて，平均値=0，標準偏差=1 になるように標準化する。

> (d <- data.frame(x=c(1, 3, 4, NA, 7), y=c(3, NA, 7, 3, 2)))
   x  y
1  1  3
2  3 NA
3  4  7
4 NA  3
5  7  2
> (d2 <- data.frame(scale(d)))
     x          y
1 -1.1 -0.3382407
2 -0.3         NA
3  0.1  1.4657098
4   NA -0.3382407
5  1.3 -0.7892283
> apply(d2, 2, mean, na.rm=TRUE)
            x             y
-6.938894e-18  0.000000e+00
> apply(d2, 2, sd, na.rm=TRUE)
x y
1 1

第 2 には，lm 関数は分析に関与する変数（独立変数と従属変数）のどれかが NA のものは分析に使われない。
従って，lm(y ~ x, data=d2) には，2 行目と 4 行目のデータは使われない。
2 行目と 4 行目のデータを除いて d2 の平均値と標準偏差を見るとどうなるか。

> apply(d2[-c(2, 4),], 2, mean, na.rm=TRUE)
        x         y
0.1000000 0.1127469
> apply(d2[-c(2, 4),], 2, sd, na.rm=TRUE)
       x        y
1.200000 1.193201

分析に使われるデータは「標準化されたものではない」。

正しくやるためにはどうするか。
まず，実際に分析に使用する変数だけを含むデータフレームを作り，na.omit なりを使って欠損値を含まないデータフレームにする。その後，scale で標準化する。そのデータフレームを lm で分析する。

> (d3 <- data.frame(scale(na.omit(d2))))
              x          y
1 -1.000000e+00 -0.3779645
3  1.156482e-17  1.1338934
5  1.000000e+00 -0.7559289
> apply(d3, 2, mean, na.rm=TRUE)
            x             y
-7.016595e-17 -1.848565e-17
> apply(d3, 2, sd, na.rm=TRUE)
x y
1 1
> summary(lm(y ~ x, data=d3))

Call:
lm(formula = y ~ x, data = d3)

Residuals:
      1       3       5
-0.5669  1.1339 -0.5669

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.205e-17  8.018e-01   0.000    1.000
x           -1.890e-01  9.820e-01  -0.192    0.879

Residual standard error: 1.389 on 1 degrees of freedom
Multiple R-squared: 0.03571,    Adjusted R-squared: -0.9286
F-statistic: 0.03704 on 1 and 1 DF,  p-value: 0.879

さらには，普通は標準化偏回帰係数はこんな風にしては求めない（2 回も重回帰分析をする必要はない）。

ダメ出し：なぜ oneway.test を使わないのだろうか?

kosugitti

> result.aov <- aov(height ~ class, data = sample)
> summary(result.aov)
            Df Sum Sq Mean Sq F value Pr(>F)
class        2    116   58.21   0.706  0.496
Residuals   97   8000   82.47               

何故こんなに回りくどいことをやるのだろうか?

しかも，これは「各群の分散が等しい」ことを仮定している。

t.test を使うなら，oneway.test を使うべきであろう。これは，「各群の分散が等しいことを仮定しない」がデフォルトになっている（t.test と平仄が合っている）。

> oneway.test(height ~ class, data=sample)

    One-way analysis of means (not assuming equal variances)

data:  height and class
F = 0.6541, num df = 2.000, denom df = 64.509, p-value =
0.5233

> oneway.test(height ~ class, data=sample, var.equal=TRUE)

    One-way analysis of means

data:  height and class
F = 0.7058, num df = 2, denom df = 97, p-value = 0.4962

ダメ出し：par(new=TRUE) は使わない

# P 君の最初の位置と、A さん B さんそれぞれの位置
P <- c(0, 0)
A <- c(0, 10)
B <- c(6, 8)

library(gplots)
niter <- 1000
a <- 1/2
b <- 1/2
pos <- matrix(NA, nr=niter, nc=2)
pos[1, ] <- a * (A - P)
for(i in 2:niter){
    if(odd(i)){
        pos[i, ] <- (pos[i - 1,] - P) + a * (A - pos[i - 1, ])
    }else{
        pos[i, ] <- (pos[i - 1,] - P) + b * (B - pos[i - 1, ])
    }
}
pos <- rbind(P, pos)

tri <- rbind(P, A, B)
plot(tri, xlim=range(tri[, 1]), ylim=range(tri[, 2]), xlab="", ylab="")
# par(new=TRUE)
# plot(pos, type="l", xlim=range(tri[,1]), ylim=range(tri[,2]), xlab="", ylab="")
lines(pos)
abline(0, pos[niter - 1, 2] / pos[niter - 1, 1], lty=2)
abline(0, pos[niter, 2] / pos[niter,1], lty=2)
polygon(tri[, 1], tri[, 2])

関数言語の特性を生かしたプログラミング

2進64桁を16進16桁に変換 にて

# 64桁の01の数値列
set.seed(1)
a <- sample(0:1,64,replace =TRUE)
# 4桁の01の数値列を0,1,2,...fに対応づける
b <- as.matrix(expand.grid(rep(list(0:1),4)))
b <- b[,4:1]
d <- c(as.character(0:9),letters[1:6])
# 2進64桁を16進16桁に直す
chage64to16<- function(a){
    b <- as.matrix(expand.grid(rep(list(0:1),4)))
    b <- b[,4:1]
    d <- c(as.character(0:9),letters[1:6])
    n <- length(a)%/%4
    ret <- rep(0,n)
    for(i in 1:n){
        tmp <- apply((t(b) - a[(1+4*(i-1)):(4*i)]==0),2,prod)
        ret[i] <- d[which(tmp==1)]
    }
    ret
}
a
chage64to16(a)
を
bin2hex <- function(a) {
    as.hexmode(apply(matrix(a, 4), 2, function(byte) sum(c(8, 4, 2, 1)*byte)))
}
のように書く。
実行速度もずいぶん速くなるけど，元のプログラムの方が，「2進数・16進数の変換を学ぶ」という目的にはあっているのか?

分かりやすくプログラミングしよう

kosugitti

「国語の試験の成績を A, B, C, D にする」ということですけど ifelse は面倒。cut も面倒。findInteval がよい。ついでに，途中結果を代入したりしないで，直接 factor 関数に引き渡す。

何回も，同じ変数名（しかも長い）を入力するのは面倒（間違っても直せばよいだけだけど）。あとでプログラムを読み直すときに長い名前はわずらわしい。longlongname と longlongname2 なんてあると，間違えるのが必至。

というようなことで，以下のように。

だいじなこと。「先生！わたし，８０点なのに B になります」

分かりやすくプログラミングするとミスも減る。念のためにテストデータで正しい答えになることを確認。

> sample <- data.frame(kokugo = c(0, 1, 58:62, 68:72, 78:82, 98:100))
> sample$kokugo.class <- ifelse(sample$kokugo > 80, 1, ifelse(sample$kokugo > 70, 2, ifelse(sample$kokugo > 60, 3, 4)))
> sample$kokugo.class <- factor(sample$kokugo.class, labels = c("A", "B", "C", "D"))
> sample$kokugo.class2 <- factor(cut(sample$kokugo, breaks = c(-1, 59, 69, 79, 101)), labels = LETTERS[4:1])
> sample$kokugo.class3 <- factor(findInterval(sample$kokugo, c(60, 70, 80)), labels = LETTERS[4:1])
> sample
   kokugo kokugo.class kokugo.class2 kokugo.class3
1       0            D             D             D
2       1            D             D             D
3      58            D             D             D
4      59            D             D             D
5      60            D             C             C
6      61            C             C             C
7      62            C             C             C
8      68            C             C             C
9      69            C             C             C
10     70            C             B             B
11     71            B             B             B
12     72            B             B             B
13     78            B             B             B
14     79            B             B             B
15     80            B             A             A
16     81            A             A             A
17     82            A             A             A
18     98            A             A             A
19     99            A             A             A
20    100            A             A             A

 

ダメ出し：オンラインヘルプ以前の問題

kosugitti

χ二乗検定をするにはchisq.test関数

chisq.test(data) 

## ## Pearson's Chi-squared test 
## ## data: data 
## X-squared = 16.67, df = 2, p-value = 0.0002404 

イェーツの補正を行わない場合

chisq.test(data, correct = FALSE) 

## 
## Pearson's Chi-squared test 
## 
## data: data 
## X-squared = 16.67, df = 2, p-value = 0.0002404 

２×２分割表以外では correct 引数は無意味

プログラミングの練習ならよいけれど

oddcount <- function(x){ # 奇数の数をカウントする
    k <- 0　
    for(n in x){
        if (n%%2 == 1) k<- k+1　# %%は剰余
    }
    return(k)
}
x <- sample(500, 1000, replace=TRUE)
oddcount(x)

以下で十分

length(x[x%%2])

sum(x%%2) で十分と書こうと思ったら既にコメントが付いていた..（呵々大笑）