『Rの基礎とプログラミング技法』メモ書き

石田基広先生による訳の『Rの基礎とプログラミング技法』を読んだメモ書きです。
あくまで自分のためのメモという類のものなので、見づらかったらすみません。

第1章 Rとは何か
「はじめにSがあった」
S言語の仕様書籍は、その表紙の色で呼ばれることが多い。
最初のS(1984)…『茶本』
バージョン2(1988)…『青本』
バージョン3(1992)…『白本』
バージョン4(1998)…『緑本』
第2章 基礎の習得
オブジェクト名は文字で始まり、間には数字や(演算子を除く)記号を含んでよい。
ただしこれを満たさない名前のオブジェクトをつくることも可能である。
その場合、バッククォーテーション(`)を用いて囲むことでアクセスする。
> (=A=) <- 100
エラー:   予想外の '=' です  ( "(=" の)
> `(=A=)` <- 100
> `(=A=)`
[1] 100
> "(=A=)" # これだと単に(=A=)というcharacterを打っただけ
[1] "(=A=)"
> "こう書いてるのと同じこと"
[1] "こう書いてるのと同じこと"
NAを対象とした計算結果は再びNA
ただし結果が一意に定まる場合にはこの限りではない。
> NA + 1
[1] NA
> NA && TRUE
[1] NA        # NA如何で結果が変わりうるため、欠損値が返る
> NA && FALSE
[1] FALSE     # NAがなんであれFALSE
ベクトルを初期化する際には、空の角括弧を用いる。
> (x <- 1:3)
[1] 1 2 3
> x[] <- 0 # xの各要素に0を代入
> x
[1] 0 0 0
> x <- 1 # xを1というスカラで上書き
> x
[1] 1
Rにおける行列や配列は、ベクトルにdimという属性を持たせることで実装されている。
よって数学では「ベクトルが配列の」特殊な場合だが、
Rでは「配列がベクトルの」特殊な場合(すなわちなくてもよいdimという属性を持っている場合)である。
> x <- 1:6
> (x <- 1:6)
[1] 1 2 3 4 5 6
> dim(x) <- c(2, 3)
> x
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
with()はデータフレームを受け取り、そこから一時的に構築された環境内でのデータハンドリングを可能にする。
> # irisのデータにおいて、がく片の幅が2.5未満で、花弁の幅が1.0より大きいサンプルの種
> # わざわざiris$Speciesのように書かなくてOK!
> with(iris, Species[Sepal.Width < 2.5 & Petal.Width > 1.0])
[1] versicolor versicolor versicolor versicolor virginica
Levels: setosa versicolor virginica
同様にsubset()は、データフレームからデータを部分取出しするのに適す。
すなわちsubset()は、第2引数の論理式を一時環境内で評価し、合致した行からなるデータフレームを返す。
> subset(iris, Sepal.Width < 2.5 & Petal.Width > 1.0)
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
54           5.5         2.3          4.0         1.3 versicolor
69           6.2         2.2          4.5         1.5 versicolor
81           5.5         2.4          3.8         1.1 versicolor
88           6.3         2.3          4.4         1.3 versicolor
120          6.0         2.2          5.0         1.5  virginica
要因ごとにデータフレームを分割するにはsplit()が便利
ただし上記のwith()subset()と違い一時環境は構築されないので、指定の仕方に注意
> split(iris, Species)
 以下にエラー split.default(seq_len(nrow(x)), f, drop = drop, ...) :
   オブジェクト "Species" は存在しません
> split(iris, iris$Species) # これでOK
$setosa
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1           5.1         3.5          1.4         0.2  setosa
2           4.9         3.0          1.4         0.2  setosa
				:
				:
$versicolor
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
51           7.0         3.2          4.7         1.4 versicolor
52           6.4         3.2          4.5         1.5 versicolor
				:
				:
$virginica
    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
101          6.3         3.3          6.0         2.5 virginica
102          5.8         2.7          5.1         1.9 virginica
				:
				:(結果一部のみ表記)
if...else構文はその機能の定義上ベクトル化されていない。
よって論理値ベクトルを与えると、その最初の値だけを分岐判断に用いる。
一方ifelse()関数はベクトル化されており、各判定での真・偽の場合の返り値もベクトルで指定しておける。
> if(c(TRUE, FALSE, TRUE)) { print("これはダメ") }
[1] "これはダメ"
Warning message:
In if (c(TRUE, FALSE, TRUE)) { :
   条件が長さが2以上なので,最初の一つだけが使われます 
> ifelse(c(TRUE, FALSE, TRUE), c("1回目T","2回目T","3回目T"), c("1回目F","2回目F","3回目F"))
[1] "1回目T" "2回目F" "3回目T"
switch()による条件判定も可能である。
このとき条件の指定は「ケースの~番目という数字」「ケースの名前」のどちらも可
存在しないケースを指定すると、無名のケースとして与えたデフォルト値が使われる。
> switch(2, a=101, b=102, c=103, 200)
[1] 102
> switch("c", a=101, b=102, c=103, 200)
[1] 103
> switch("z", a=101, b=102, c=103, 200)
[1] 200 # デフォルト値
ベクトルxの各要素への操作をforループを用いて行なうとき、繰り返しは1:length(x)ではなくseq(along=x)で指定されるべき。
なぜならxが長さ0のベクトルのとき、1:length(x)c(1, 0)を生成してしまい、ループが回ってしまうからである。
> x <- c(1,3,5,7)
> for(i in seq(along=x)) { print("cycle") }
[1] "cycle"
[1] "cycle"
[1] "cycle"
[1] "cycle"
> for(i in 1:length(x)) { print("cycle") } # この場合は同じ
[1] "cycle"
[1] "cycle"
[1] "cycle"
[1] "cycle"
> x <- integer(0)
> # xに要素がないのだから、ループは回らないで欲しい
> for (i in seq(along=x)) { print("cycle") } # printは一度も実行されない(ループに入っていない)
> for (i in 1:length(x)) { print("cycle") } # iがc(1, 0)の間、2回ループが回っている
[1] "cycle"
[1] "cycle"
第3章 データの入力と出力
dump()関数はオブジェクト名を受け取り、そのオブジェクトを構成するためのスクリプトを生成する。
> dump("women", file="women.txt")
> cat(readLines("women.txt"), sep="\n") # women.txtの内容の確認
`women` <-
structure(list(height = c(58, 59, 60, 61, 62, 63, 64, 65, 66, 
67, 68, 69, 70, 71, 72), weight = c(115, 117, 120, 123, 126, 
129, 132, 135, 139, 142, 146, 150, 154, 159, 164)), .Names = c("height", 
"weight"), row.names = c(NA, -15L), class = "data.frame")
> # 実行すればdump()に示したオブジェクトが生成できるようなコードが保存されている
第4章 R言語の詳細
ピリオド仮引数"..."は実行時に与えられた実引数のうち、他の仮引数に割り当てられなかったものをすべてうける。
関数内においてこれらの仮引数は"..."によってさらに他の下位関数に渡すことが可能
このような仕様は、上位関数から多くの設定を引き継ぐ必要のあるグラフィック関連の関数で多用される。
> func <- function(x, ...) { matrix(x, ...) }
> func(1:6, ncol=2, byrow=TRUE)
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6
Rにおけるワークスペースは.GlobalEnvという環境であり、サーチパスの一番前に位置する。
関数を呼び出すとこの.GlobalEnvの前に、関数独自の環境が生成される。
この仕組みにより、関数内において生成・操作されたオブジェクトが同名のワークスペースのオブジェクトに干渉することはない。
ただしRは対象がみつかるまでサーチパスを順に検索するため、関数内で生成していない名前のオブジェクトが、意図しないかたちで他の環境から呼ばれてしまうことがある。
> search() #現在のサーチパスの確認
[1] ".GlobalEnv"        "package:stats"     "package:graphics"  "package:grDevices" "package:utils"
[6] "package:datasets"  "package:methods"   "Autoloads"         "package:base"     
> x <- 5 # .GlobalEnvに生成される
> scope <- function(do.it) {
+     if (do.it) x <- 3
+     inner <- function() print(x)
+     inner()
+ }
> scope(TRUE) # scope()関数内でx <- 3が実行されているため、inner()実行時に直近の環境であるscope()の環境から呼ばれる
[1] 3
> scope(FALSE) # xはinner()環境にもscope()環境にもみつからず、.GlobalEnv環境においてみつかったxが使われる
[1] 5
名前空間をもつパッケージにおいては、明示的にエクスポートされていない限り、同じ名前空間内の関数からしかアクセスされない。
これにより名前空間内では自由にオブジェクトを命名できる。
ある名前空間のオブジェクトにアクセスすることを明示的に示したい場合は"::"演算子を用いる。
また名前空間からエクスポートされていない関数に関しても、getFromNamespace()や":::"演算子によりアクセス可能
一般にオブジェクトの定義やその値の取得は、"<-"による代入およびオブジェクト名の表記によってなされる。
これらを多数自動的に行なうには、assign()get()を用いる。
ただし言うまでもなく、このようなオブジェクトをそれぞれつくるよりもリストとして構成するほうが管理がラク
> rm(list=ls())
> for (i in 1:3) assign(paste("letter", i, sep=""), LETTERS[i]) # LETTERSはRのビルトイン定数
> ls()
[1] "i"       "letter1" "letter2" "letter3"
> letter2
[1] "B"
ある関数Aから関数Bを呼び出したとき、関数B内でエラーが生じると、関数Aも直ちに終了されてしまう。
try()関数を用いれば、エラー終了を抑制し、下位関数でエラーが生じたときの対応を上位関数で変化させられる。
try()は、与えた式がエラーが生じなかった場合には単にその評価結果を返す。
エラーが生じた場合、try-errorクラスのオブジェクトを返す。
> B <- function(x) { return(x+1) }            # 入力値に1を足して返すだけ
> A <- function(x) {
+     y <- B(x)^2                             # B()の返り値を二乗し
+     cat("関数Bのあとに実行されています\n")  # コンソールに文字を出力して
+     return(y)                               # おしまい
+ }
> A(3)
関数Bのあとに実行されています
[1] 16                                        # (3+1)^2
> A("ムリぽ")                                 # B()でエラー終了するとそれ以下(cat())も実行されない
以下にエラー x + 1 :  二項演算子の引数が数値ではありません
>
> A <- function(x) {
+     y <- try(B(x)^2)                             # B()の結果に二乗しようとしてみて
+     cat("関数Bのあとに実行されています\n")       # コンソールに文字を出力し
+     if (class(y)=="try-error") y <- "エラー発生"  # もしB()が失敗していたらyを書き換え
+     return(y)                                    # 返す
+ }
> A("ムリぽ")
Error in x + 1 :  二項演算子の引数が数値ではありません 
関数Bのあとに実行されています                      # B()は失敗してもcat()は実行されている
[1] "エラー発生"                                   # 値も返っている(A()は強制終了されていない)
第5章 効率的なプログラミング
「読者はベクトルを処理単位としたプログラミングを心得てほしい!」
「読者はベクトルを処理単位としたプログラミングを心得てほしい!」
大切なことだから、二度言いました(・ω・)
実装済みの関数の多くは、CやFortranなどによって高速なライブラリ化されている。
わざわざ作り直す意味はない。
また必要があれば、自作の関数を上記のような他言語に移植することも可能
そのためにRには、コンパイルされたライブラリを作成し、リンクする手段が用意されている。
プログラミングは
(1) 汎用性:可能な限り汎用的に
(2) 利便性:多くのコメントや、場合によってはドキュンメンテーションを含み
(3) 可読性:簡潔すぎず、冗長すぎず、また適度にタブや半角スペースでみやすくし
(4) 効率性:速度に関して最適に
に注意して書かれるべきである。
これはたとえば
(1) 当座は必要ないパラメータも引数化→上位関数から変更したくなったときにすぐ指定可能
(2) ほとんどの場合コメントは、存在しないか、あっても簡潔すぎて意味が分からない
(3) カンマのあとには半角スペースを入れ、異なる操作は別の行にする。
(4) どんな場合でも(たとえその時点では不要でも)ベクトル単位のプログラミングを心掛ける
などのことを念頭においておくことを意味している。
ループは避けるに越したことはないが、すべてのプログラムがループなしで実装できるわけではない。
たいせつなのは、ループにおいて何をすると時間的・空間的計算量で損をするかを把握すること。
ループの結果を格納する変数は、長さが分かっているなら先に初期化しておく。
毎回の処理ごとに要素を足していくようなプログラムは、メモリ領域の確保と解放を繰り返すので、効率が悪い。
長さが分からない場合でも、十分な長さの変数を確保しておき、最後に使わなかった分を削った方が、計算量的に良い場合さえある。
> a <- NULL
> system.time( for(i in 1:100000) a <- c(a, i) )
   ユーザ   システム       経過  
     86.25       0.06      89.91 
> 
> a <- numeric(100000)
> system.time(for (i in 1:100000) a[i] <- i)
   ユーザ   システム       経過  
      1.15       0.00       1.18 
引数の適切性のチェックはループに入る前に、計算結果のチェックはループ終了後に、いずれもベクトル単位で行なう。
ループ内で1回の計算ごとに毎回行なうのはばかばかしい。
同じ理由で、すべての結果に対して行なう計算もループ終了後に、ベクトル単位で行なう。
たとえば割合の結果を百分率にしたいとき、"*100"はループの終了後にできた計算結果に対して行ない、毎回のループのなかで行なったりしない。
(こう言われると「当たり前じゃん」と思うが、意外にやってしまうので注意)
apply()系関数をうまく利用して、可読性・効率性をますことができる。
apply(X, MARGIN, FUN, ...):行列・配列の次元ごとのベクトル処理
lapply(X, FUN, ...):リストの要素ごとでの処理
sapply(X, FUN, ...)lapply()と同じだが、結果をリストよりシンプルな形(ベクトルや行列)で返す
tapply(X, INDEX, FUN, ...):INDEXで指定されたfactorに基づいて、グループごとでの処理
mapply(FUN, ...):多数の引数をとる関数でのベクトル処理
これらの関数群はCライブラリによって最適化されており、速い。
ただし「どんな場合でもループよりapply()系関数が速くて良い」というのは、迷信
計算対象にもよるが、ループも適切に組めばこれらとほとんど変わらないことも。
時間を測定する関数を利用して、コードの速度を計測する。
system.time(評価式)は、コードの実行にかかった時間を表示する。
Rprof(); 評価式; Rprof(NULL)は、ワーキングディレクトリにRprof.outというレポートファイルを作成する。
これはsummaryRprof()によって参照でき、コード内の使用時間の内訳をみることができる。
Rで使用可能なメモリ量はmemory.limit()で確認できる。
またmemory.limit(MB)によって、R起動後にメガバイト単位で使用可能メモリを変更することも可能。
(本ではmemory.limit()による「表示」はバイト単位で「指定」はメガバイト単位とあるが、どちらもメガバイト単位になった模様)
これらの設定は実行時のオプション「 --max-mem-size」でも変更可なので、必要があればショートカットで設定しておく。
Rをしばらく使用していると、削除済みオブジェクトに割り当てられていたメモリなど、不要なメモリ領域がたまって不安定になる。
gc()によってそのようなメモリを解放できる。
また直前の結果を保存している隠しオブジェクト.Last.valueが大きい場合、gc()を2回連続で行なうとよい。
一度目のgc()によって.Last.valuegc()の返り値で上書きされ、二度目のgc()で先の大きな.Last.valueが使用していた領域が解放されるからである。

ちなみにわたしは、gc()を2回ではなぜかちゃんと解放されないので、invisible(replicate(20, gc()))みたいな使い方をしている。
有効性は不明だが、とりあえずほっておくと落ちてしまうプログラムが通るようになった。
> gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 105974  2.9     350000  9.4   350000  9.4
Vcells  74630  0.6     786432  6.0   353908  2.8
> x <- rnorm(10000000)
> rm(x)
> gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 106125  2.9     350000  9.4   350000  9.4
Vcells  74970  0.6    8658815 66.1 10075926 76.9 # 66Mb!!
> gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 106132  2.9     350000  9.4   350000  9.4
Vcells  74993  0.6    6927051 52.9 10075926 76.9 # まあちょっと減ったが
> 
> 
> invisible(replicate(20, gc()))
> gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 106150  2.9     350000  9.4   350000  9.4
Vcells  75212  0.6     786432  6.0 10075926 76.9 # やっと起動直後レベル
第6章 オブジェクト指向プログラミング