![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
COLOR(red){SIZE(25){R プログラミング Tips 大全}} #contents R でプログラミングする際役に立つと思われる Tips を脈絡無く集めます。なお、ある 程度溜ったら適当にカテゴリー化 しますので宜しく。また、お寄せ頂いた Tips は 参照に適当と思われる場合は、合体・ 修正を勝手に行う場合がありますので予めご承知ください。また、各自ご協力宜しく。 *本項について [#o6b41d5e] #contents ~ R でプログラミングする際役に立つと思われる Tips を脈絡無く集めます。なお、ある程度溜ったら適当にカテゴリー化しますので宜しく。 また、お寄せ頂いた Tips は 参照に適当と思われる場合は、合体・修正を勝手に行う場合がありますので予めご承知ください。 各自ご協力宜しく。 *Rの関数定義の基本 [#ve5e8148] [[Rの関数定義の基本]] 参照 *Rの基本データ構造：ベクトル・行列・配列・リスト [#v99667c8] [[ベクトル・行列・配列・リストTips大全]] を参照。 *R 言語の実行制御フロー [#k311bf23] R は多くの計算機言語と同じような Algol 風制御命令のセットをもつが、より柔軟である。 実行文 expr は単純実行文でも、(波括弧で括った)複合実行文(同一行に並べるにはセミコロンで区切る)でもよい。 ** 繰り返し for [#he751a77] 書式 (ループ範囲 range の各要素 arg に対して expr を実行 ) for(arg in range) expr 注意：for ループは一般に実行速度を遅くするボトルネックになりやすい。またコードが長くなり勝ちである。apply 関数ファミリの使用や、特にベクトル・行列・配列の成分ごとのループは専用高速関数が用意されているのでその使用を考える。 ループ範囲にベクトルを取る(基本) > x = 1:4 > for (i in x) cat(i,"\n") 1 2 3 4 ループ範囲に文字ベクトルを取る > x = c("a", "b", "1", "2") > for (i in x) cat(i,"\n") a b 1 2 ループ範囲に行列を取る(実際はベクトルとしてアクセス) > x = matrix(1:4, c(2,2)) > for (i in x) cat(i, "\n") 1 2 3 4 ループ範囲に配列を取る(実際はベクトルとしてアクセス) > x=array(1:8, c(2,2,2)) > for (i in x) cat(i, "\n") 1 2 3 4 5 6 7 8 ループ範囲にリストを取る > x = list("a", 1:4, 3) > for (i in x) cat(i, "\n") a 1 2 3 4 3 ループ範囲は最初に与えられたもので実行される。実行文中でループ変数を変更しても影響をうけない。 > for (i in 1:4) {i <- i+2; cat(i, "\n")} 3 4 5 6 # i=1,2,3,4 で4回実行されている > x=1:4 > for (i in x) {x[i] <- i+10; cat(i,"\n")} 1 2 3 4 # 実行文中での範囲変数 x の変更はループ範囲を変更しない > x [1] 11 12 13 14 # x 自身は変更される ## next はループを一回パスする > for (i in 1:3) + for (j in 1:3) {if (j==2) next; cat(i, j, "\n")}} 1 1 1 3 2 1 2 3 3 1 3 3 ## すべての x 中の i, j, k (i, j, k は異なる) に関する三重和 for (i in x) { for (j in x[-i]) { for (k in x[-c(i,j)]) { ..................... }}} ## または for (i in x) { for (j in x) { if (i==j) next for (k in x) { if (k==i || k==j) next ......................................... }}} (重複の無い)数値ベクトルから長さ 3 の組み合わせをループ > temp <- function(x){ + x <- sort(x) # ソートしておかないとまずい + for (i in x) + for (j in x[x>i]) + for (k in x[x>j]) + cat(i,j,k,"\n")} > temp(1:4) 1 2 3 1 2 4 1 3 4 2 3 4 > temp(runif(4)) 0.151133 0.223455 0.264883 0.151133 0.223455 0.993983 0.151133 0.264883 0.993983 0.223455 0.264883 0.993983 for を使わない方法，重複があってもよいし，ソートも不要 > (x <- runif(4)); t(combn(4, 3, function(i) x[i])) [1] 0.2769620 0.2577154 0.7625690 0.1760667 [,1] [,2] [,3] [1,] 0.2769620 0.2577154 0.7625690 [2,] 0.2769620 0.2577154 0.1760667 [3,] 0.2769620 0.7625690 0.1760667 [4,] 0.2577154 0.7625690 0.1760667 ** 条件実行 if, if else [#b139ac47] # condがTRUEのとき、exprを実行する。 if (cond) expr # condがFALSEのとき、alt.exprを実行する。 if (cond) cons.expr else alt.expr # if else 構文は入れ子にできる。 if (cond1) expr1 else if (cond2) expr2 else expr3 # if else 構文は値を返す。if (x < 0) y <- NA else y <- sqrt(x) と同じ。 x <- 1; y <- if (x < 0) NA else sqrt(x) 注意：ifelse(cond, expr1, expr2) 関数は if(cond) expr1 else expr2 構文を簡略化したもの 簡略化したものではない!! > x <- -3:3 > (a <- ifelse(x < 0, NA, sqrt(x))) [1] NA NA NA 0.000000 1.000000 1.414214 1.732051 警告メッセージ： In sqrt(x) : 計算結果が NaN になりました > (b <- if(x < 0) NA else sqrt(x)) [1] NA 警告メッセージ： In if (x < 0) NA else sqrt(x) : 条件が長さが2以上なので，最初の一つだけが使われます ** 条件が満たされる限り実行 while [#sda71355] 書式 (条件、実行文に関する注意は if 文と同様) while(cond) expr 注意：条件 cond は必ず一度は評価される > i=0 > while( (i<-i+1) <= 5 ) cat(i, "\n") 1 2 3 4 5 > i # 実行後の i の値 [1] 6 > i = j = 0 > while( (i <- i+1) < 1 ) {j <- j+1; cat(i, "\n")} > i [1] 1 # cond は一回実行され i の値は 1 になっている > j # expr は実行されていない [1] 0 ** 単純繰り返し repeat [#oc258f0d] 書式 (R 自身を中断しない限り、repeat 文を中断する唯一の手段は expr 中の break 文) repeat expr ** 実行中断 break [#b1ddf9a2] 書式 (プログラムの実行を中断し、これを含む制御文の次の実行文に移る) break > i=0 > repeat {if ((i <- i+1) == 5) {cat(i,"\n"); break}} # i==5 になったら中断 5 指定秒数プログラムをアイドリングする関数(実はこのため用の関数 Sys.sleep がある -> その他の項参照) > Sleep <- function (n) { startingtime <- unclass(Sys.time()) repeat { if (unclass(Sys.time())-startingtime >= n) break } } > system.time(Sleep(20)) [1] 18.66 0.46 19.12 0.00 0.00 # 19.12 秒中断 > system.time(Sleep(10)) [1] 9.07 0.24 9.31 0.00 0.00 # 9.31 秒中断 ** 次の実行サイクルへ飛ぶ next [#o3a5495d] 書式 (制御文の現在の実行を中断し、次の実行サイクルに移る) next > for (i in 1:10) {if (i%%2==0) next else cat(i,"\n")} # 偶数ならスキップ 1 3 5 7 9 ** 多重選択 switch [#md2f0835] 書式 (EXPR は整数もしくは文字列(を与える式)) switch(EXPR, ...) 任意個数の引数 ... は選択肢リスト もし EXPR が整数値ならば該当番号の選択肢が評価され、返される もし EXPR が文字列を返すなら、それが選択肢である名前付き項目とマッチングされ、該当する項目の値が返される。 もしマッチした項目の値が欠損していれば、その次の項目の値が返される。 もしマッチする項目が一切なければ、switch にそれ以外の引数があればそれが、さもなければ NULL が返される。 > switch("cc", a=1, cc=, d=2) 2 # 項目 cc= にマッチするが値が与えられて無いので、その次の値のある項目の値 2 が返される > switch(3, 11,12,13) [1] 13 # 3番目の数値が返される > centre <- function(x, type) { switch(type, mean = mean(x), median = median(x), trimmed = mean(x, trim = .1)) } > x <- rcauchy(10) > centre(x, "mean") # 文字列 "mean" にマッチした項目の値 mean(x) が返される [1] -1.163613 > ccc <- c("b","QQ","a","A","bb") > for(ch in ccc) cat(ch,":",switch(EXPR = ch, a=1, b=2:3), "\n") b : 2 3 # ch="b" は項目 b=2:3 にマッチ QQ : # ch="QQ" はマッチする項目が無いので NULL が返された a : 1 A : bb : b : 2 3 > for(ch in ccc) cat(ch, ":", switch(EXPR = ch, a=, A=1, b=2:3, "Otherwise: last"), "\n") b : 2 3 QQ : Otherwise: last # マッチする名前付き項目がないので最後の文字列が返された a : 1 A : 1 bb : Otherwise: last ## EXPR が数値の時は「その他」項目は許されない > for(i in c(-1:3,9)) print(switch(i, 1,2,3,4)) NULL # i=-1 に該当する項目は無いので NULL が返された NULL [1] 1 [1] 2 [1] 3 NULL * 関数引数関係 [#l6e65392] ** missing() 関数(仮)引数が存在するかどうかチェック [#p6311080] R 関数の仮引数の数は場合により変わり得るので、ある引数が実際に存在するかどうかチェックする必要が起こる 書式 (x は関数の仮引数) missing(x) 注意：missing(x) は x が関数中で変更されると信頼できなくなる(例えば x <- match.arg() の後では必ず FALSE になる。 myplot <- function(x, y) { if(missing(y)) { y <- x # もし仮引数 y が与えられなければ y <- x とする x <- 1:length(y) } plot(x,y) } ** 関数引数のマッチング match.arg [#ld898533] 関数の引数に与えられた値が候補のどれに一致するかどうかチェックする 書式 (choice 引数は候補リストのベクトル) match.arg(arg, choices) 部分的マッチングが行なわれるが、返り値は非省略形で与えられる。該当する項目がなければエラーメッセージがでる center <- function(x, type = c("mean", "median", "trimmed")) { type <- match.arg(type) switch(type, # match.type の結果により switch mean = mean(x), median = median(x), trimmed = mean(x, trim = .1)) } x <- rcauchy(10) > center(x, "t") # 部分的マッチングで type="trimmed" とされた [1] 0.478473 > center(x, "med") # 部分的マッチングで type="median" とされた [1] 0.4619571 > center(x, "m") # エラー(m で始まる候補が二つある) Error in match.arg(type) : ARG should be one of mean, median, trimmed ** 関数の形式引数を得る、設定する [#v10bfc9f] > f <- function(x) a+b # 関数の定義 > fomrals(f) # 関数 f の形式的引数は？ $x # x > formals(f) <- alist(a=, b=3) # 形式的引数 x を a ,b (b は既定値 3 付き) > f # f の定義を確認 function (a, b = 3) a + b > f(2) [1] 5 > f(2, 5) [1] 7 > f <- function() x > formals(f) <- alist(x = , y = 2, ... = ) # 「その他大勢」引数 ... の設定 > f function (x, y = 2, ...) x * 関数本体の操作 [#r3170f3d] ** 関数本体を表示、設定する body() (2004.1.31) [#je825414] > foo <- function() x # 関数を定義 > body(foo) # 関数本体を表示 x > body(foo) <- expression(x^2) # 関数本体を再定義 > foo # 確かに変わっている function () x^2 fd <- deriv(~ 100*(x2 - x1*x1)^2 + (1-x1)^2, c("x1", "x2")) # 関数の数式一階偏微分 fdd <- function(x1, x2) {} # 空の関数定義 body(fdd) <- fd # 関数の本体を fd とする * 論理判断 [#j94dd634] ** 基本論理演算子 [#h33732c0] 注意：これらはベクトル化されている # 一致・大小判断は数値ベクトルのみならず R のオブジェクトにもそれぞれの意味に応じた定義がされていることがある #文字列なら辞書式順序での一致・大小判断の意味になる x == y # 一致 x != y # 不一致 x > y # x が大きい x >= y # x が大きいか等しい x < y x<= y x || y # 論理和 (x が TRUE か y が TRUE なら TRUE) x && y # 論理積 (x, y が ともに TRUE なら TRUE) !x # x の否定 x & y # ベクトルの論理積 x | y # ベクトルの論理和 > x <- (1:4)>2 > x [1] FALSE FALSE TRUE TRUE > !x [1] TRUE TRUE FALSE FALSE > "abc" > "abd" # FALSE > "abc" > "aba" # TRUE > "abc" > "abca" # FALSE > x <- rnorm(10) > x > 0 [1] TRUE TRUE TRUE FALSE TRUE TRUE FALSE FALSE TRUE TRUE > x[x > 0] [1] 2.0583593 0.5014932 0.8126365 0.1641549 0.6721926 0.3478276 1.6111553 注意：if 文等の検査で ==, != を使うことは避け、identical 文を使うべきである。if 文では単一の TRUE, FALSE が要求される~ 注意：また数値誤差の観点から実数の正確な比較は無意味になることがある > x1 <- 0.5 - 0.3 > x2 <- 0.3 - 0.1 > x3 <- 0.4 - 0.2 > x1 == x2 # ほとんどの機械で FALSE になる > x1 == x3 # (私の機械では) TRUE > identical(all.equal(x1, x2), TRUE) # いつも TRUE > a<-c(TRUE,FALSE,TRUE,FALSE) > b<-c(TRUE,TRUE,FALSE,FALSE) > a & b [1] TRUE FALSE FALSE FALSE > a | b [1] TRUE TRUE TRUE FALSE ** 安全な完全一致判断。 identical() [#w275f4d4] 二つの R オブジェクトが完全に一致するかどうかを検査する安全で信頼できる方法。 注意：通常の一致検査 x==y, x != y はしばしば困難をもたらす。x,y が異なったタイプである場合や、長さが異なる場合等。また返り値がベクトルになる可能性がある。 # 書式 identical(x, y) # x, y は二つの R オブジェクト > if(identical(all.equal(x,y),TRUE) # R お勧めの二つのオブジェクトの一致検査 > if (x == y) # 好ましくない一致検査 > identical(1, 1.) ## R では 1 と 1. はともに倍精度実数で表現 [1] TRUE > identical(1, as.integer(1)) ## as.integer(1) は整数型なので FALSE になる [1] FALSE > x <- 1.0; y <- 0.99999999999 > identical(all.equal(x, y), TRUE) # 安全な一致比較法 [1] TRUE > all.equal(options(), .Options) [1] "Modes: list, pairlist" # 一致しないので食い違いが表示される ** すべて該当するか？all() [#wd0ceaf1] すべてが TRUE なら TRUE 書式 all(x) x は論理値ベクトルもしくはそれを与える式 range(x <- sort(round(rnorm(10) - 1.2,1))) if(all(x < 0)) cat("all x values are negative\n") ** ほとんど該当するか？ all.equal() [#g6c2b6f7] 二つの R オブジェクトが「殆んど等しいか」どうかをチェックする。「殆んど等しい」の意味はオブジェクトにより異なる。数値ベクトルでは誤差の範囲無いで一致するかどうかを判断する。 書式 (all.equal.numeric は数値オブジェクトに対する all.equal のメソッド) all.equal(target, current, ...) all.equal.numeric(target, current, tolerance= .Machine$double.eps ^ 0.5, scale=NULL, ...) tolerance は誤差範囲の指定パラメータ # R で階乗 n! を計算する基礎は gamma(n+1) として計算することであるが、 # gamma(n+1) の値は内部的に実数計算されるので、真の意味で n! に一致しないことを注意する > all.equal(gamma(2:14), cumprod(1:13)) # 既定の誤差範囲無いで TRUE [1] TRUE > all (gamma(2:14) == cumprod(1:13)) # 真の一致は無いので FALSE [1] FALSE > all.equal(gamma(2:14), cumprod(1:13), tol=0) # 食い違いを見る [1] "Mean relative difference: 5.198216e-15" > all.equal(options(), .Options) [1] "Modes: list, pairlist" ** 属性が一致するか？ attr.all.equal() (2004.08.19) [#v1f377c7] ## 以下を比較せよ > all.equal(1, c(a=1)) # c(a=1) は名前属性付きの数値 1 [1] TRUE > identical(1, c(a=1)) [1] FALSE > attr.all.equal(1, c(a=1)) [1] "names for current but not for target" > 1 == c(a=1) a TRUE ** ひとつでも該当するか？ any() [#o266d461] どれかが TRUE なら TRUE 書式 any(x) x は論理値ベクトルもしくはそれを与える式 range(x <- sort(round(rnorm(10) - 1.2,1))) if(any(x < 0)) cat("x contains negative values\n") ** すべてが TRUE でなければ stop する。 stopifnot() [#jd417fa1] # 書式 stopifnot(A, B, ...) # A, B, ... は真偽値を持つ任意個数のオブジェクト # A, B, ... のどれかが FALSE なら最初の FALSE となった項目を表示して stop する > stopifnot(1 == 1, all.equal(pi, 3.14159265), 1 < 2) # すべてが TRUE なので stop しない > m <- matrix(c(1,3,3,1), 2,2) > stopifnot(m == t(m), diag(m) == rep(1,2)) # すべてが TRUE なので stop しない > stopifnot(1==1, 2< 2) Error: 2 < 2 is not TRUE # stop した > stopifnot(1:4==c(1,2,3,5)) # stop した Error: 1:4 == c(1, 2, 3, 5) is not TRUE > stopifnot(1:4==c(1,2,3,4)) # stop 無し *通し番号付きのオブジェクト名をプログラム中で生成 assign() [#f3455df4] 文字ベクトルや数値を用いてあるオブジェクト名をプログラム中で生成し、値を付値する には > assign(paste("file", 1, "max", sep=""), 1) > ls() [1] "file1max" 変数 var1,var2,...,var10 を生成し、それぞれ値 1,2,...,10 を代入 for (i in 1:10) assign(paste("var", i, sep=""), i) file[1],...,file[10] はデータファイル名の文字列であるとして、次の 関数はそれを読み込み、順に f1,f2,...,f10 という名前のオブジェトに付値 temp <- function() { for (i in 1:10) { assign(paste("f",i,sep=""), read.table(file[i],skip=1)) } } 順に通し番号付きの名前 d1, d2, d3 の変数を作り、それに値 x1-x2, y1-y2, z1-z2 の値を代入 (r-help 記事より、2004.09.23) x1 <- rnorm(10); x2 <- rnorm(10) y1 <- rnorm(10); y2 <- rnorm(10) z1 <- rnorm(10); z2 <- rnorm(10) names. <- c("x", "y", "z") # 文字ベクトル for(i in names.){ # (順に)変数 x1, y1, z1 の値を res1 に代入 res1 <- get(paste(i,"1",sep="")) # (順に)変数 x2, y2, z2 の値を res2 に代入 res2 <- get(paste(i,"2",sep="")) # (順に)変数 d1, d2, d3 を作り、それに値 res1-res2 を代入 assign(paste("d",i,sep=""), res1-res2) } *日付・時刻 date(), Sys.time() [#k8599892] > date() [1] "Fri Mar 21 15:25:05 2003" > Sys.time() [1] "2003-03-21 15:25:12 JST" > unclass(Sys.time()) [1] 1048227920 # この数値は1970年元旦(Unix 元年)からの経過 *条件分岐 ifelse() [#o4d24cd8] 論理値ベクトル test に対し ifelse(test, x, y) は test[i]==TRUE なら x[i], test[i]==FLASE なら y[i] を並べた test と同じ長さのベクトルを返す。関数の値を変数値に応じて変えるために便利。 > x <- (-1):3 > ifelse(x>0, log(x), NA) # x <=0 なら NA を返 [1] NA NA 0.0000000 0.6931472 1.0986123 Warning message: # 途中で log(x) を計算するので NaN が発生、警告がでる NaNs produced in: log(x) > log(ifelse(x>0,x,NA)) # x=-1,0 では 先ず x <- NA となり、ついで log(NA)=NA が返される [1] NA NA 0.0000000 0.6931472 1.0986123 # 警告無し > ifelse(x>=0, x, -x) # abs(x) と同じ [1] 1 0 1 2 3 * apply 関数ファミリ apply, mapply, lapply, sapply, tapply [#mbbf5f3c] apply 関数ファミリーは一つの関数を複数のオブジェクトに適用して得られた結果をベクトル、行列、リストの形で一括して返す。for ループを陽に使ったり、複数オブジェクトを個別に与える必要がなく、コードが簡潔になる。また関数中で使用すれば、 R 関数プログランミングの基本精神「意味がある時は(ない時も)常にベクトル化せよ」を簡単に実現できる。 ** 配列のあるマージンに関数を適用 apply() [#qa7b69a9] -書式 COLOR(red){apply(X, MARGIN, FUN, ...)} > x <- matrix(1:12, ncol=3) > x [,1] [,2] [,3] [1,] 1 5 9 [2,] 2 6 10 [3,] 3 7 11 [4,] 4 8 12 > apply(x,1,max) # 各行の最大値を求める [1] 9 10 11 12 > apply(x,2,max) # 各列の最大値を求める [1] 4 8 12 # 例 (行列の各列毎に正規化) # apply(x, 2, mean), apply(x, 2, sd) で行列 x の列毎の平均、標準偏差からなるベクトルをつくる。 # scale(x, a, b) は行列 x の第 i 列から a[i] を引き、b[i] で割った結果を表す x と同じサイズの行列を計算する。 # つまり xx[i,j] = ( x[i,j] - a[j] ) / b[j] > x <- matrix(runif(100), c(20,5)) > xx <- scale(x, apply(x, 2, mean), apply(x, 2, sd)) # チェック > apply(xx, 2, mean) # 数値誤差範囲内で各列の平均零 [1] -1.554312e-16 5.551115e-17 -1.915135e-16 1.387779e-16 -1.609823e-16 > apply(xx, 2, sd) # 標準偏差はすべて零 [1] 1 1 1 1 1 # scale 関数は列毎の操作専用なので、行毎に正規化したければまず転置しておいてから、 # 最後にまた転置する。 > xx <- t( scale(t(x), apply(t(x), 2, mean), apply(t(x), 2, sd)) ) ** mapply 関数。結果 FUN(x[1],y[1]), FUN(x[2],y[2]),...,FUN(x[n],y[n]) を一括して返す [#hd19cc61] 関数 FUN を任意(個数)引数 ... 中の値を順に引数にとって得られるベクトルのリストを返す -書式 COLOR(red){mapply(FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE, USE.NAMES = TRUE)}。例えば mapply(FUN, 1:4, 1:2) は FUN(1,1), FUN(2,2), FUN(3,1), FUN(4,2) を並べたリスト(足りない引数はリサイクル規則が適用される) # ベクトル rep(1,4), rep(2,3), rep(3,2), rep(4,1) を並べたリストを返す > mapply(rep, 1:4, 4:1) [[1]] [1] 1 1 1 1 [[2]] [1] 2 2 2 [[3]] [1] 3 3 [[4]] [1] 4 # rep(4,1), rep(3,2), rep(2,3), rep(1,4) を並べたリストを返す > mapply(rep, times = 1:4, x = 4:1) [[1]] [1] 4 [[2]] [1] 3 3 [[3]] [1] 2 2 2 [[4]] [1] 1 1 1 1 # rep(42, 1), rep(42, 2), rep(42,3), rep(42,4) を並べたリストを返す > mapply(rep, 42, 1:4) [[1]] [1] 42 [[2]] [1] 42 42 [[3]] [1] 42 42 42 [[4]] [1] 42 42 42 42 # mapply(rep, 1:4, 42) はすべての行ベクトルが 1:4 であるサイズ 42x4 の行列になる > x <- matrix(1:10, byrow = TRUE, nrow = 2) > x [,1] [,2] [,3] [,4] [,5] [1,] 1 2 3 4 5 [2,] 6 7 8 9 10 > mapply(function(x1, x2) x1+x2, x[1, ], x[2, ]) # 列和 [1] 7 9 11 13 15 > mapply(function(x1, x2) x1-x2, x[1, ], x[2, ]) # 列差 [1] -5 -5 -5 -5 -5 ** lapply, sapply 関数。リストもしくはベクトルの各成分に関数を適用した結果を一括して返す [#oc814bc1] -書式 COLOR(red){lapply(X, FUN, ...)}, COLOR(red){sapply(X, FUN, ..., simplify = TRUE, USE.NAMES = TRUE)}。リストもしくはベクトル変数 x の各成分に順に関数 FUN を適用して得られた結果をリストで返す。sapply は結果をベクトル、行列の形で返す。 # リストの三つの成分の平均値からなるリストを計算(論理値は整数 1,0 と見做される) > x <- list(a = 1:10, beta = exp(-3:3), logic = c(TRUE, FALSE, FALSE, TRUE)) > lapply(x, mean) $a [1] 5.5 $beta [1] 4.535125 $logic [1] 0.5 > lapply(x, quantile, probs = 1:3/4) $a 25% 50% 75% 3.25 5.50 7.75 $beta 25% 50% 75% 0.2516074 1.0000000 5.0536690 $logic 25% 50% 75% 0.0 0.5 1.0 # x の中央値を計算(結果は行列の形になる) > sapply(x, quantile) a beta logic 0% 1.00 0.04978707 0.0 25% 3.25 0.25160736 0.0 50% 5.50 1.00000000 0.5 75% 7.75 5.05366896 1.0 100% 10.00 20.08553692 1.0 # 数列からなるリストの各成分の５数要約からなる行列を計算 > str(i39 <- sapply(3:9, seq)) List of 7 $ : int [1:3] 1 2 3 $ : int [1:4] 1 2 3 4 $ : int [1:5] 1 2 3 4 5 $ : int [1:6] 1 2 3 4 5 6 $ : int [1:7] 1 2 3 4 5 6 7 $ : int [1:8] 1 2 3 4 5 6 7 8 $ : int [1:9] 1 2 3 4 5 6 7 8 9 > sapply(i39, fivenum) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1.0 1.0 1 1.0 1.0 1.0 1 [2,] 1.5 1.5 2 2.0 2.5 2.5 3 [3,] 2.0 2.5 3 3.5 4.0 4.5 5 [4,] 2.5 3.5 4 5.0 5.5 6.5 7 [5,] 3.0 4.0 5 6.0 7.0 8.0 9 ** tapply 関数。 因子で定義されるグループ毎に関数を適用した結果を一括して返す [#l4e16427] -書式 COLOR(red){tapply(X, INDEX, FUN = NULL, ..., simplify = TRUE)}。X は普通ベクトルで INDEX は X の要素をグループに分ける因子の組み合わせ。各グループに関数 FUN を適用した結果をベクトルもしくはリストで返す > n <- 17 > fac <- factor(rep(1:3, len = n), levels = 1:5) > fac [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 # 17個のデータに与えた因子値 Levels: 1 2 3 4 5 # 因子水準は5とする > table(fac) # table 関数で因子水準毎に度数を集計 fac 1 2 3 4 5 # 数列 1:17 中の因子値が 1,2,3,4,5 であるものの度数はそれぞれ 6.6.5.0.0 6 6 5 0 0 > tapply(1:n, fac, sum) # 結果を(ラベル付き)ベクトルで返す 1 2 3 4 5 # 因子値が 1,2,3 であるグループの総和はそれぞれ 51,57,45 51 57 45 NA NA # 因子値が 4,5 のグループは空なので値 NA が返される > tapply(1:n, fac, sum, simplify = FALSE) # 結果をリストで返す $"1" [1] 51 $"2" [1] 57 $"3" [1] 45 $"4" NULL $"5" NULL > tapply(1:n, fac, prod) # 積の場合 1 2 3 4 5 58240 209440 29160 NA NA ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// * ベクトル、行列、配列の(一般化)外積 outer、クロネッカー積 kronecker [#a92a9177] [[配列の外積、クロネッカー積]] を参照 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////// * R の基本(数値)関数 [#i2ec7fa1] R における他の関数と同様にここで紹介する関数はベクトル化されており、引数は数値・整数ベクトルが可能で、結果も対応するベクトルになる。 ////////////////////// ** 基本数値演算子 [#df9e9ec2] |x + y | | |x - y | | |x * y | | |x / y | | |x ^ y | 1^y, y^0 は常に 1 になる。x,y が 正負の inf の時も可能なら合理的な値が返る| |x %% y| x mod y| |x %/% y | いわゆる整数商で値は整数。関係 x = (x%/%y)*y + x %%y が成立 | 注意：x％％y は適当な整数 n で 0<= x - ny < y となる時の値 x - ny を表す。その時の整数 n が x%/%y になる。 > x [1] -1 0 1 2 3 4 5 6 7 8 9 10 11 12 > x%%5 [1] 4 0 1 2 3 4 0 1 2 3 4 0 1 # すべて 0,1,2,3,4 の範囲の値 > x%/%5 [1] -1 0 0 0 0 0 1 1 1 1 1 2 2 2 > (x%/%5)*5+(x%%5) [1] -1 0 1 2 3 4 5 6 7 8 9 10 11 12 # x を復元 > x=1+0.5i; y= 3+2i > c(x+y, x-y, x*y, x/y, x^y) # 複素数も問題無し [1] 4.000000000000000000+2.50000000000000000i [2] -2.000000000000000000-1.50000000000000000i [3] 2.000000000000000000+3.50000000000000000i [4] 0.307692307692307709-0.03846153846153846i [5] -0.023927552111312745+0.55238103053669563i > c(x%%y, x%/%y) # これはさすがに駄目 Error: unimplemented complex operation ///////////////////////////// ** 初等数値関数 [#n4cd6033] ***三角関数ファミリー(引数の単位はラディアンでベクトル化されている) [#c21a1a0e] |cos(x) | | |sin(x) | | |tan(x) | | |acos(x) | arccos | |asin(x) | arcsin | |atan(x) | arctan | |atan2(y, x) | x 軸と、原点と座標 (x,y) を結ぶベクトルのなす角 | | | x、y>0 なら atan2(y,x) = atan(y/x)| *** hyperbolic 関数ファミリー [#od3b35c3] |cosh(x) | (exp(x) + exp(-x))/2 | |sinh(x) |(exp(x) - exp(-x))/2 | |tanh(x) | sinh(x)/cosh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x)) | |acosh(x) | cosh(x) の逆関数 | |asinh(x) |sinh(x) の逆関数 | |atanh(x) |tanh(x) の逆関数 | *** 対数、指数関数ファミリー [#y59b1d45] 注意：底 base は実数。引数 x は実数もしくは複素数 | log(x, base = exp(1)) | 底 base (既定値 exp(1)) の対数 | | logb(x, base = exp(1)) | log(x, base=exp(1)) と同じに見えるが？| | log10(x) | 常用対数 log(x, base=10) | | log2(x) | 底 2 の対数 | | exp(x) | 指数関数 | | expm1(x) | x の絶対値 が1より(かなり)小さい時、exp(x)-1 をより正確に計算| | log1p(x) | x の絶対値が1より(かなり)小さい時、log(1+x) をより正確に計算| | | x が -1 に近い時は逆により不正確 | > log(-1) [1] NaN Warning message: NaNs produced in: log(x) ## log も exp も所詮近似計算！精度が問題になることも？ > options(digits=22) # 表示桁数最大化 > log(0.001+1) [1] 0.0009995003330834232 > log1p(0.001) ## |x| が小さければ log1p の方が相対的に正確(13桁目で食い違い) [1] 0.000999500333083533 > exp(0.01)-1 [1] 0.010050167084167949 > expm1(0.01) ## |x| が小さければ expm1 の方が相対的に正確(14桁目で食い違い) [1] 0.010050167084168058 ## 複素数もOK > log(1+2i) [1] 0.80471895621705+1.10714871779409041i > exp(1+2i) ***組合せ論的関数 (2004.09.23) [#zee970e2] -COLOR(red){choose(n, k)} n 個から k 個を取り出す場合の数 -COLOR(red){lchoose(n, k)} choose(n,k) の自然対数 -COLOR(red){factorial(x)} x の階乗 x! (R 1.9 より登場。実際は単に gamma(x+1) を計算するだけで、整数でない x も可能) -COLOR(red){lfactorial(x)} factorial(x) の自然対数 (R 1.9 より登場。実際は単に lgamma(x+1) を計算するだけで、整数でない x も可能) ***番外 [#h328fd71] -COLOR(red){sign(x)} 実数 x の符号。x < 0, = 0, x > 0 に応じて sign(x) = -1,0,1 -COLOR(red){abs(x)} 実数もしくは複素数 x の絶対値 -COLOR(red){sqrt(x)} 整数の平方根、もしくは複素数の平方根(主値) > sqrt(-1) # 負実数の平方根は存在しないので NaN が返される [1] NaN Warning message: NaNs produced in: sqrt(-1) > sqrt(-1+0i) # しかし複素数と考えれば難なく計算 [1] 0+1i ***数値ベクトルに対する逐次処理関数 (2004.09.23) [#ye526a03] -COLOR(red){cumsum(x)} 累積和ベクトルを返す -COLOR(red){cumprod(x)} 逐次積ベクトルを返す -COLOR(red){cummax(x)} 逐次最大値ベクトルを返す -COLOR(red){cummin(x)} 逐次最小値ベクトルを返す > (x <- sample(1:10)) [1] 6 7 1 9 5 2 3 10 8 4 > cumsum(x) [1] 6 13 14 23 28 30 33 43 51 55 > cumprod(x) [1] 6 42 42 378 1890 3780 11340 113400 907200 [10] 3628800 > cummin(x) [1] 6 6 1 1 1 1 1 1 1 1 > cummax(x) [1] 6 7 7 9 9 9 9 10 10 10 //////////////////////////// ** 超越関数 [#wbca0b51] ***ガンマ関数ファミリー [#v983f4a7] -COLOR(red){beta(a, b)}。ベータ関数 B(a,b) = (gamma(a)gamma(b))/(gamma(a+b))。引数 a,b は零および負整数を除く数値(のベクトル)。 -COLOR(red){lbeta(a, b)}。ベータ関数の自然対数。直接計算しており、log(beta(a,b)) として計算しているわけではない。 -COLOR(red){gamma(x)}。ガンマ関数。引数 x は零および負整数を除く数値(のベクトル)。正の整数 n に対する階乗 n! は gamma(n+1) として計算するのが基本。 -COLOR(red){lgamma(x)}。ガンマ関数の自然対数。直接計算しており、log(gamma(x)) として計算しているわけではない。 -COLOR(red){digamma(x)}。lgamma の一階微分。 -COLOR(red){trigamma(x)}。lgamma の二階微分。 -COLOR(red){tetragamma(x)}。lgamma の三階微分。 -COLOR(red){pentagamma(x)}。lgamma の四階微分。 -COLOR(red){choose(n, k)}。二項係数(n 個から k 個を選ぶ場合の数)。gamma(n+1)/(gamma(k+1)gamma(n-k+1)) として計算しているので n, k が整数でなくても値は得られるが、結果は整数化されるようである(そもそも意味はなさそうであるが)。 -COLOR(red){lchoose(n, k)}。二項係数 choose(n,k) の自然対数。直接計算しており、log(choose(n,k)) として計算しているわけではない。 *** ベッセル関数ファミリー [#rddb5fc0] 引数 x は非負値、次数 nu は実数(負の値の場合を含む)~ もし `expon.scaled = TRUE' ならば桁溢れ、桁落ちを避けるため exp(-x) I(x;nu) もしくは exp(x) K(x;nu) を計算~ 参考: Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions.~ Dover, New York; Chapter 9: Bessel Functions of Integer Order. -COLOR(red){besselI(x, nu, expon.scaled = FALSE)} Modified Bessel function of first kind. -COLOR(red){besselK(x, nu, expon.scaled = FALSE)} Modified Bessel function of second kind. -COLOR(red){besselJ(x, nu)} Bessel function of first kind. -COLOR(red){besselY(x, nu)} Bessel function of second order. -COLOR(red){gammaCody(x)} ガンマ関数(gamma 関数とはコードが異なり、計算が少し早いが gamma よりは精度が落ちる) #ref(bessel1.png) #ref(bessel2.png) #ref(bessel3.png) #ref(bessel4.png) #ref(bessel5.png) #ref(bessel6.png) #ref(bessel7.png) //////////////////////////// ** 丸め関数 [#q23fcdce] 注意：R における数値丸め処理の基本は round 関数であるが、これは IEEE 規約に基づき、いわゆる「四捨五入」とは微妙に異なる点が有るので注意が必要。詳しくは [[JIS,ISO式四捨五入]] を参照。 -COLOR(red){ceiling(x)} 「天井関数」。x は数値ベクトル、結果は整数ベクトル。 x 未満でない最小の整数 > y=ceiling(1.3) > y [1] 2 > is.integer(y) # y は整数を表す数値であり、整数保管モードではないことを注意 [1] FALSE > is.numeric(y) [1] TRUE > ceiling(1.0) [1] 1 > ceiling(-1.0) [1] -1 > ceiling(-1.3) [1] -1 -COLOR(red){ floor(x)} 「床関数」。x は数値ベクトル、結果は整数(をあらわす数値ベクトル)。x 以上でない最大の整数 (いわゆるガウス記号 [x] が計算する整数) > floor(c(1.3, 1.0, 0, -1.1)) [1] 1 1 0 -2 -COLOR(red){ round(x, digits = 0)} 丸め関数。digits で指定した(小数点以下)桁で IEEE 式丸めを行なう。既定では digits = 0、つまり小数以下を丸める。IEEE 方式は「一番近い偶数に丸める」が基本で、場合によると「5捨」になることもある。 > round(c(1.32, 1.5, -0.2, -1.2, -1.5)) [1] 1 2 0 -1 -2 > round(23.5) [1] 24 > round(23.52, 1) # 小数点以下第１桁で丸める(i.e. digits =1) [1] 23.5 > round(23.52, -1) # 小数点以上第１桁で丸める [1] 20 > round(c(23.5, 24.5)) # round(24.5) は 25 でないことを注意() [1] 24 24 > round(1.3-1.2i) # 複素数も受け付ける(実・虚部をそれぞれ丸める?) [1] 1-1i -COLOR(red){ signif(x, digits = 6)} digits で指定された有効桁数(既定値 6 桁)に丸める > x2 <- pi * 100^(-1:3) > round(x2, 3) [1] 0.031 3.142 314.159 31415.927 3141592.654 > signif(x2, 3) # 有効数字3桁になる [1] 3.14e-02 3.14e+00 3.14e+02 3.14e+04 3.14e+06 -COLOR(red){ trunc(x)} x を「0 へ向かって」整数化する > trunc(c(1.2, 3, 3.1, -0.2, -0.3, -0.31, -1.3)) [1] 1 3 3 0 0 0 -1 > trunc(c(0.2, -0.2)) # 「一番近い整数」に変換 [1] 0 0 > trunc((-10):10/3) [1] -3 -3 -2 -2 -2 -1 -1 -1 0 0 0 0 0 1 1 1 2 2 2 3 3 -COLOR(red){ zapsmall(x, digits= getOption("digits"))} round(x, digits = dr) の丸め位置桁 dr を「0 に近い値」がゼロになるように定めて、丸める(今一意味が分かりにくいが?) > x2 <- pi * 100^(-1:3) > x2/1000 [1] 3.141593e-05 3.141593e-03 3.141593e-01 3.141593e+01 3.141593e+03 > zapsmall(x2/1000, digits = 4) [1] 0.0 0.0 0.3 31.4 3141.6 > round(x2/1000, digits = 4) [1] 0.0000 0.0031 0.3142 31.4159 3141.5927 > zapsmall(exp((0+1i) * 0:4 * pi/2)) # 複素数も可 [1] 1+0i 0+1i -1+0i 0-1i 1+0i *集合演算 [#s702e58b] x, y 等は(同一モードの)ベクトルで、重複は無いとする。もし重複があっても、union, intersect, setdiff は重複の無い結果を返す。 | union(x, y) | 和集合| | intersect(x, y) |積集合、共通部分| | setdiff(x, y) |差集合 | | setequal(x, y) |集合として等しいか？| | is.element(el, set) | el は set に含まれるか？| | %in% | el %in% set は is.element(el, set) の簡易形| | unique | 重複のある集合から重複元を除いて一意化する| > (x <- c(sort(sample(1:20, 9)), NA)) [1] 1 2 5 10 12 15 17 19 20 NA > (y <- c(sort(sample(3:23, 7)), NA)) [1] 4 5 6 9 13 14 16 NA > union(x, y) # 和集合 [1] 1 2 5 10 12 15 17 19 20 NA 4 6 9 13 14 16 > intersect(x, y) # 共通集合 [1] 5 NA > setdiff(x, y) # 差集合(x の要素の内、y に属さないものの全体) [1] 1 2 10 12 15 17 19 20 > setdiff(y, x) [1] 4 6 9 13 14 16 > setequal(x, y) # 集合としての同等性検査 [1] FALSE > setequal(1:3, 3:1) [1] TRUE # 集合としては等しい！ > is.element(x, y) # x の各要素毎に y に所属するかどうか検査 [1] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE > y %in% x [1] FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE > x <- unique(c(1:5, 3:7, 8:2)) [1] 1 2 3 4 5 6 7 8 # 重複要素を一意化 * ソート、オーダー、ランク関数 sort(), order(), rank() [#f3514e96] [[ソート、オーダー、ランク]] を参照 * 基本統計処理関数群 [#jbd74732] max(..., na.rm=FALSE) min(..., na.rm=FALSE) pmax(..., na.rm=FALSE) pmin(..., na.rm=FALSE) range(..., na.rm = FALSE) range.default(..., na.rm = FALSE, finite = FALSE) mean(x, ...) mean.default(x, trim = 0, na.rm = FALSE, ...) weighted.mean -COLOR(red){ave(x, ..., FUN = mean)} 因子レベル毎に x の部分集合を平均。~ 引数：x は数値ベクトル。... は x と同じ長さのグループ化変数(群)。普通因子。FUN は各因子レベルの組合せに適用される関数(既定で平均値)。~ 返り値：x と同じ長さのベクトル y。 もし ... が二つのグルーピング変数 g1,g2 なら、y[i] = FUN( {x[k]; g1[k]=g1[i] & g2[k]=g2[i]} )。 > ave(1:5) # グルーピング変数無し。x[1],...,x[5] 毎に mean(1:5) が返される [1] 3 3 3 3 3 > ave(1:3, factor(c(1,1,2))) # y[1]=y[2]=(x[1]+x[2])/2, y[3]=x[3] [1] 1.5 1.5 3.0 > ave(1:3, factor(c(1,2,2))) # y[1]=x[1], y[2]=y[3]=(x[2]+x[3])/2 [1] 1.0 2.5 2.5 > ave(1:3, factor(c(1,2,2)), FUN=max) # y[1]=max(x[1]), y[2]=y[3]=max(x[2],x[3]) [1] 1 3 3 > data(warpbreaks); attach(warpbreaks) > warpbreaks breaks wool tension 1 26 A L 2 30 A L 3 54 A L --- 途中省略 --- 53 16 B H 54 28 B H > wool [1] A A A A A A A A A A A A A A A A A A A A A A A A A A A B B B B B B B B B B B [39] B B B B B B B B B B B B B B B B Levels: A B > ave(breaks, wool) # wool の種類 A, B 毎の平均 [1] 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 [9] 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 [17] 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 [25] 31.03704 31.03704 31.03704 25.25926 25.25926 25.25926 25.25926 25.25926 [33] 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 [41] 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 [49] 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 > unique(ave(breaks, wool)) [1] 31.03704 25.25926 # A, B に対する二種類の平均 median(x, na.rm=FALSE) quantile(x, probs = seq(0, 1, 0.25), na.rm = FALSE, names = TRUE, ...) fivenum(x, na.rm = TRUE) IQR mad cumsum colSUms var(x, y = NULL, na.rm = FALSE, use) cor(x, y = NULL, use = "all.obs") cov(x, y = NULL, use = "all.obs") sd(x, na.rm = FALSE) cov.mt diff(x, ...) diff.default(x, lag=1, differences=1, ...) *R のオブジェクト・命令を表す文字列を評価実行 [#q06776d6] ** 構文 eval(parse(text=...)) [#hce83178] String2Eval を R のオブジェクト・命令を表す文字列とすると eval(parse(text = String2Eval)) はそれを解釈実行する > eval(parse(text="ls()")) > eval(parse(text = "x[3] <- 5")) > hoge <- function(a, b) { x1 <- paste("temp <- function(x) {",sep="") # 関数 temp を表現する文字列の一部 x2 <- paste("a*x + b} ",sep="") # 関数 temp を表現する文字列の一部 eval(parse(text=paste(x1,x2,sep=""))) # 関数 temp(x) を定義 temp(9) } # temp(9) を即座に実行する > hoge(1,2) # 内部で関数 temp(x) が定義され、temp(9) を計算する [1] 11 # 以下でも同じこと > hoge <- function(a, b) { x1 <- paste("temp <- function(x) {",sep="") x2 <- paste(a, "*x + ", b, "} ",sep="") eval(parse(text=paste(x1,x2,sep=""))) temp(9) } # 匿名関数を利用する例 hoge <- function(a, b) { x1 <- paste("(function(x) { ", sep="") # 関数 temp を表現する文字列の一部 x2 <- paste("a*x + b} )(9)", sep="") # 関数 temp を表現する文字列の一部 y <- eval(parse(text=paste(x1,x2,sep=""))) # 関数 temp(x) を定義 y } hoge(1,2) [1] 11 # 内部で定義された関数は hoge 終了時には消える。hoge を呼び出した環境中に定義するには次のようにする # (temp(x) は hoge 中でも実行可能) > rm(temp) > hoge <- function(a, b) { x1 <- paste("temp <- function(x) {",sep="") x2 <- paste(a, "*x + ", b, "} ",sep="") eval.parent(parse(text=paste(x1,x2,sep=""), n=1)) } > hoge(1,2) > temp # 関数 temp(x) が定義されて（残って）いる function(x) {1*x + 2} # 応用問題。通し番号つきの関数をプログラム中で作成 > hoge <- function(a, b) { # a,b は共通パラメータ for (i in 1:3) { code <- paste("temp", i, " <- function(x) ", a, "*x^(", i, ")+ ", b, sep="") eval.parent(parse(text=code), n=1) } } > hoge(1,2) > temp1 function(x) 1*x^(1)+ 2 > temp2 function(x) 1*x^(2)+ 2 > temp3 function(x) 1*x^(3)+ 2 **文字列を与えて、それを名前に持つオブジェクトの値を得る, get() 関数 (2004.1.31) [#q2dd65c5] > xyz <- 1:5 > get("xyz") [1] 1 2 3 4 5 > eval(parse(text="xyz")) # これでも良い [1] 1 2 3 4 5 > foo <- function(x) x^2 > get("foo") function(x) x^2 > eval(parse(text="foo")) # これでも良い function(x) x^2 * 関数から複数の値を返す (リスト返り値) [#e0fd8647] 一つの関数から複数の値を返すにはベクトル、配列等のオブジェクトを返り値とすれば良い。return 関数に複数の引数を与えると、それらは自動的にリストとして返される。リストの各成分には元の変数名が名前タグとして自動的に付加される。~ COLOR(red){注意：この機能は将来は廃止の予定、R 1.8 からは警告が出る。}) 代わりに返り値として明示的にリストを与える。 > test <- function(x){y=log(x);z=sin(x);return(x,y,z)} > test(1:3) $x [1] 1 2 3 $y [1] 0.0000000 0.6931472 1.0986123 $z [1] 0.8414710 0.9092974 0.1411200 名前タグを自前で指定するには明示的に名前付きリストを返り値にする > test <- function(x){y=log(x);z=sin(x);return(list(value=x,log=y,sin=z))} > test(1:3) $value [1] 1 2 3 $log [1] 0.0000000 0.6931472 1.0986123 $sin [1] 0.8414710 0.9092974 0.1411200 * 関数内部での関数定義 [#qee54dda] R の関数内部には別の関数定義を置き、その場で実行できる "test" <- function() { "hello" <- # 関数内部での関数定義 function() { x <- 1:3 y <- 23:34 z <- cbind(x,y) return(x,y,z) } c <- hello() # その実行 return(c$z,c$x, c$y) } *永続付値 <<- の秘密 [#xec7f997] 永続付値演算子は現在の環境の親R環境中のオブジェクトに対する操作であり、特に関数中で用いるときは注意がいる。永続付値演算子を関数中で用いるとき陥りやす間違いは次の例でみられる: > v <- 1 > test <- function () {v <- 2; v <<- 3; cat("v=",v,"\n")} > test() v= 2 > v [1] 3 この例では関数内部での変数 v と永続付値で生成された親環境(つまり関数 test を 実行した環境)中の変数 v という、同じ名前だが実体が異なる二つの変数が登場する。関数中では単なる変数 v は局所的な変数を指し、関数実行後は残らない。一方変数v は関数実行後も残る #もう少し微妙な例: > test <- function () {v <- 1:3; v[2] <<- 0; cat("v=",v,"\n")} > test() v= 1 2 3 # 局所的な v の値 > v [1] 1 0 3 # 大局的な v の値(いつの間にかベクトルになっていることに注意) 教訓: 永続付値は便利だが使用には注意がいる。使うときはできるだけ個性的な名前 を用いる。決して親環境(場合によるとそのまた親環境等)、もしくは関数中の局所的変数と紛らわしい名前は使わない! 特に R の大局的環境(起動時の環境)は COLOR(red){.GlobalEnv} という名前が与えられている。大局的環境中のオブジェクトは任意の関数から直ちに参照でき、特に副関数に複雑な引数を(参照渡し、コピー無し)引き渡すとき便利である。大局的環境中の変数の操作には明示的に COLOR(red){assign} 関数を用いるのが安全な方法。 assign(x, 2, env=.GlobalEnv) # 大局的環境中の変数 x に値 2 を代入 多くの場合これは x <<- 2 と同等の効果を持つが、環境が何重にも入れ子になっている場合は致命的なエラーを(気づかぬままに)侵す可能性がある。 ** 一次的作業環境の使用 new.env(TRUE,NULL) [#ua70a97c] 永続付値の項で述べた同じ名前でも環境(簡単にいえばメモリースペース)が異なれば違 た実体を持つという点は次の積極的に仮の環境を使う例で一層はっきりする。多少面倒だが多数の関共有する変数はこうした一次的作業環境に置いておくと混乱が無いかも知れない。またこれは R でオブジェクトのコピー無しで他から参照する(参照渡し)を実現する手段でもある。 > rm(v) # 大局的環境に変数 v があれば抹消 > tempenv <- new.env(TRUE,NULL) # 新しい環境を作り tempenv と名付ける > assign("v",3,env=tempenv) # 環境 tempenv 中に変数 v (注意：文字列で与える)を作り、値を 3 とする > v Error: Object "v" not found # 大局的環境には v という変数は無いのでエラー > ls(envir=tempenv) # 環境 tempenv 中の変数一覧 [1] "v" # v という変数がある! > get("v",env=tempenv) # 環境 tempenv 中の変数 v の値を得る [1] 3 > rm(tempenv) # 環境 tempenv は大局的環境 .GlobalEnv 中のオブジェクト、rm 関数で抹消できる > ls(envir=tempenv) Error in ls(envir = tempenv) : Object "tempenv" not found * コマンドラインから文字列を入力 readline() [#mc47eaf5] ** キー入力があるまでプログラム、プロットの実行を停止 readline() [#mf0ad4cb] キー入力があるまでプログラム、プロットの実行を停止。指定メッセージを出力し 何かをキー入力したら以下を実行 > temp <- function(x) { cat("sum of x is -->", sum(x), "\n") readline("Can I proceed? ") cat("prod of x is -->", prod(x),"\n") } > temp(1:4) sum of x is --> 10 Can I proceed? #ここで何かをキー入力すると次を実行 prod of x is --> 24 ** 次のプロットを出力する前に確認のプロンプトを提示する par(ask=TRUE) [#zcc03118] これは example 関数によるデモを実行する際、ひき続くグラフィックス表示を 一つずつ眺めるためにも有効である。par(ask=FALSE) を実行すれば元に戻る ** コマンドラインからオブジェクト名を入力 (2004.1.31) [#e96115aa] > a [1] 0.557755694 0.756125389 0.062838319 0.116221907 0.287878570 0.490844365 [7] 0.988744111 0.681430858 0.681793307 0.008137233 > y <- readline() a > y # y は文字列でそのままではオブジェクト名ではない [1] "a" > eval(parse(text=y)) # 文字列 y をオブジェクト名 a に変換し評価した [1] 0.557755694 0.756125389 0.062838319 0.116221907 0.287878570 0.490844365 [7] 0.988744111 0.681430858 0.681793307 0.008137233 > get(y) # 同じこと。get 関数は文字列を与え、対応するオブジェクトの値を返す [1] 0.557755694 0.756125389 0.062838319 0.116221907 0.287878570 0.490844365 [7] 0.988744111 0.681430858 0.681793307 0.008137233 **数値をコマンドラインから入力する [#m6c45052] readline 関数による入力は文字列となるので、数値を入力するにはひと工夫がいる。 > x <- readline() 1.2 # 1.2 を入力 > x [1] "1.2" # 文字列になっている > xx <- as.numeric(x) # 数値に強制変換 > xx [1] 1.2 もし複数の数値をいっぺんに入力したければ次のようにします。つまり readline 関数で読み込んだ複数数値を表す文字列(半角空白一つを間にはさむ)を、strsplit 関数で文字列リストに変換、unlist 関数で文字列ベクトルに変え、最後に as.numeric で数値に強制変換します。 > x <- as.numeric(unlist(strsplit(readline(), split=" "))) 1.0 2 -3.9 > x [1] 1.0 2.0 -3.9 もし入力に複数の半角空白が分離記号として入っていれば正規表現 [[:space:]]+ (半角空白の一個以上の繰返し)を split オプションに指定。 > x <- readline() 2 3 1.14 > x [1] "2 3 1.14" # 入力先頭の半角空白は無視される > xx <- as.numeric(unlist(strsplit(x,split="[[:space:]]+"))) > xx [1] 2.00 3.00 1.14 # 整数 2,3 が欲しければすこし厄介(後から as.integer 関数を使う？) * 自前の関数を保存し、次回に使えるようにする save(), source() [#ha49b5ab] 関数の定義をテキストファイル(例えば "myfun.R")に書き source("myfun.R") 関数で 読み込むR の初期設定ファイル .Rprofile に行 source("myfun.R") を入れておけばいつも自動に読み込まれる) R セッション中に関数定義を save() 関数で Rdata ファイルに保存しておく。関数 attach() で定義ファイルを読み込むことができる > myfun <- function(x) cat(x,"\n") > save(myfun, file="myfun.Rdata") .......... 次の R セッション............ > attach("myfun.Rdata") * 数値ベクトルの対話的入力 readline() [#pfede9cb] readline 関数で指定プロンプトを提示し、対話的に文字列を入力(任意の分離記号を指定 できる)結果 Str は文字列 "10,20,30"。strsplit 関数で文字列リストに変換し、unlist 関数で文字列ベクトル化 as.numeric 関数で数値化する > Str <- readline("Enter x: ") Enter x: 10,20,30 # プロンプトと、入力数値 > x <- as.numeric(unlist(strsplit(Str, ","))) # Str は文字ベクトル > x [1] 10 20 30 # x は数値ベクトル c(10,20,30) *ファイル、ディレクトリ操作 [#d892419d] **現在のディレクトリを R 中から変更 getwd(),setwd() [#u33077be] R からファイルを操作するには完全なディレクトリを込めて指定するのが一番安全。 現在のディレクトリを R から変更することもできる。但し R から(Unix の)OS命令を実行する基本である system("命令名")は必ずしも希望の結果を得られないこともあるらしい。 > getwd() # 現在のディレクトリを表 [1] "/home/hoge" > system("cd ..") # Unix 命令 "cd .." を実行 > getwd() [1] "/home/hoge" # しかし実際は変更はされて > getwd() [1] "/home/hoge/" > setwd("..") > getwd() [1] "/home" # 変更されている > setwd("~") > getwd() [1] "/home/hoge" **ファイル操作関数あれこれ [#ea9d9da6] 以下で COLOR(red){...}, COLOR(red){file1},COLOR(red){file2}, COLOR(red){to}, COLOR(red){from} はファイル名を表す文字列(ベクトル)。COLOR(red){path} は単一のパス名を含む文字ベクトル。COLOR(red){overwrite} は論理値でファイルを上書きするかどうかを指示する。 -file.create(...) -file.exists(...) -file.remove(...) -file.rename(from, to) -file.append(file1, file2) -file.copy(from, to, overwrite = FALSE) -file.symlink(from, to) -dir.create(path) The '...' arguments are concatenated to form one character string: you can specify the files separately or as one vector. All of these functions expand path names: see 'path.expand'. 'file.create' creates files with the given names if they do not already exist and truncates them if they do. 'file.exists' returns a logical vector indicating whether the files named by its argument exist. 'file.remove' attempts to remove the files named in its argument. 'file.rename' attempts to rename a single file. 'file.append' attempts to append the files named by its second argument to those named by its first. The R subscript recycling rule is used to align names given in vectors of different lengths. 'file.copy' works in a similar way to 'file.append' but with the arguments in the natural order for copying. Copying to existing destination files is skipped unless 'overwrite = TRUE'. The 'to' argument can specify a single existing directory. 'file.symlink' makes symbolic links on those Unix-like platforms which support them. The 'to' argument can specify a single existing directory. 'dir.create' creates the last element of the path. Value: 'dir.create' and 'file.rename' return a logical, true for success. The remaining functions return a logical vector indicating which operation succeeded for each of the files attempted. Author(s): Ross Ihaka, Brian Ripley See Also: 'file.info', 'file.access', 'file.path', 'file.show', 'list.files', 'unlink', 'basename', 'path.expand'. cat("file A\n", file="A") cat("file B\n", file="B") file.append("A", "B") file.create("A") file.append("A", rep("B", 10)) if(interactive()) file.show("A") file.copy("A", "C") dir.create("tmp") file.copy(c("A", "B"), "tmp") list.files("tmp") setwd("tmp") file.remove("B") file.symlink(file.path("..", c("A", "B")), ".") setwd("..") unlink("tmp", recursive=TRUE) file.remove("A", "B", "C") * コード中の長い行の一括注釈化 if(0){...} [#m6a9e835] R コード中の # はそれ以降行末までをコメントとみなし無視する(日本語もOK) コード中の各行は先頭に # をつければ注釈できるが、長い行をひとかたまりに 注釈化するには次の構文で囲む。つまり if 文の条件が常に FALSE だから括弧内はすべ て無視される if(0) {(長い行........)} 応用：foo(x, 1) ならデバッグ用コードを実行、foo(x) ならしない foo <- function(x, debug=0) { # 既定ではデバッグ用コードを実行しない ............................................ if(debug==1) { ... デバッグ用コード ... } .............................................. } * try() 関数。エラーが起きても立往生しないで次の作業ができるようにする [#i131435d] シミュレーション等でエラーが起こる可能性のある場合、続けて次のシミュレーションを続行するのに便利。 #書式 try(expr, silent = FALSE, first = TRUE) # expr は実行可能な表現 # silent =TRUE ならエラーコードを表示しない ## シミュレーションを100回実行、エラーにならなかった場合の結果を記録 set.seed(123) # 結果を再現できるように乱数種を指定 x <- rnorm(50) # 正規乱数50個 doit <- function(x) # シミュレーションコード { x <- sample(x, replace=TRUE) # x から重複を許して(x と同じ数だけ)無作為抽出 if(length(unique(x)) > 30) mean(x) # 重複を除いて 31 個以上なら平均計算 else stop("too few unique points") # さもなければ中断 } ## doit を100回実行し、結果をリスト res に連続記録 res <- lapply(1:100, function(i) try(doit(x), TRUE)) res [[1]] # 1回目は失敗 [1] "Error in doit(x) : too few unique points\n" attr(,"class") [1] "try-error" [[2]] ＃2回目は成功 [1] 0.1529146 [[3]] # 3回目は失敗 [1] "Error in doit(x) : too few unique points\n" attr(,"class") [1] "try-error" [[4]] # 4回目は成功 [1] 0.08090636 ------- 途中省略 ------------------------ [[100]] [1] 0.04618465 > unlist(res) # 非リスト化(文字列ベクトルになる) [1] "Error in doit(x) : too few unique points\n" [2] "0.152914625907379" [3] "Error in doit(x) : too few unique points\n" [4] "0.0809063649112537" ------- 途中省略 ------ [99] "-0.0810531720129378" [100] "0.0461846522112245" * その他 [#a51a1cc7] ** 丸括弧の不思議 (付値と表示を同時に行なう) [#ude97240] > x=1:4 # 付値のみで結果の表示はない > (x=1:4) # 丸括弧で囲むと、付値した結果を表示 [1] 1 2 3 4 > test1(x) <- function() y <- 2*x > test1(1:4) # 何も表示されない > test2 <- function(x) (y <- 2*x) > test2(1:4) # 丸括弧で囲むと付値した上で表示される [1] 2 4 6 8 ** メニューによる選択 menu関数 [#uab57e41] -書式 COLOR(red){menu(choices, graphics = FALSE, title = "")}~ -引数：COLOR(red){choices}: 選択肢を表す文字列ベクトル, COLOR(red){graphics}: 未使用, COLOR(red){title}: メニューのタイトルに使われる文字列~ -返り値：選択された項目の番号 1,2,...。もし未選択(0 を入力)なら 0 を返す。 ## 選択肢番号により switch 文で実行内容を変える > switch(menu(c("A","B"), title="My Menu (0 for no selection)") + 1, cat("Nothing is chosen\n"), cat("A is chosen\n"), cat("B is chosen\n")) My Menu (0 for no selection) 1:A 2:B Selection: 1 A is chosen > switch(menu(c("A","B"), title="My Menu") + 1, cat("Nothing is chosen\n"), cat("A is chosen\n"), cat("B is chosen\n")) My Menu (0 for no selection) 1:A 2:B Selection: A A is chosen > switch(menu(c("A","B"), title="My Menu") + 1, cat("Nothing is chosen\n"), cat("A is chosen\n"), cat("B is chosen\n")) My Menu (0 for no selection) 1:A 2:B Selection: 0 Nothing is chosen ** 論理値の数値化 (r-help 記事より) 2003.11.16 [#qb7ee0b5] R の論理値 TRUE, FALSE は数値が必要とされる文脈では整数値 1, 0 に強制変換されるが、明示的に数値に変換するには次のようにする。 > x <- matrix(1:4, 2, 2) > x [,1] [,2] [1,] 1 3 [2,] 2 4 > xx <- x>2.5 > xx # 論理値行列 [,1] [,2] [1,] FALSE TRUE [2,] FALSE TRUE > xx+0 # 数値に強制変換 [,1] [,2] [1,] 0 1 [2,] 0 1 > xx + 0:0 # 保管モードを整数にする [,1] [,2] # 大規模行列ではメモリ使用量が減るメリット [1,] 0 1 [2,] 0 1 ## storage.mode(xx) <- "integer" でも良い ## as.integer 関数では行列属性を除いてしまう > as.integer(xx) [1] 0 0 1 1 R では論理値 TRUE, FALSE は文脈に応じて整数 1,0 に強制変換されるので、表示の際気持ちが悪いということがなければ、整数として使いたい場合でも特に整数に明示的に変換する必要はない。それでも変換したければ次のようなトリックが使える。 > is.integer(TRUE+FALSE) # 足し算をすれば強制変換される [1] TRUE > is.integer(TRUE*TRUE) # かけ算をすれば強制変換される [1] TRUE > is.integer(TRUE+0:0) # 0:0 は整数 0 になるからそれを足せば整数に強制変換される [1] TRUE > is.integer(TRUE*1:1) # 1:1 は整数 1 になるからそれをかけると整数に強制変換される [1] TRUE > is.real(x+0) # 0 は実数とみなされるので結果も実数に強制変換される [1] TRUE > as.integer(c(TRUE, FALSE)) # 明示的な整数への変換例 [1] 1 0 > x = matrix(c(TRUE, FALSE, FALSE, TRUE), c(2,2)) # 論理値行列 > x [,1] [,2] [1,] TRUE FALSE [2,] FALSE TRUE > x+0:0 # 整数値行列に変換 [,1] [,2] [1,] 1 0 [2,] 0 1 **ダブルコロン演算子 (2003.11.25) [#af05fa5f] -書式 COLOR(red){pkg::name} -引数 --COLOR(magenta){pkg} パッケージ名 --COLOR(magenta){name} 変数・関数名(必要に応じ二重引用府で囲む) -説明 パッケージ pkg 中の変数 name の値を返す。もしパッケージを未ロードならばロードされる。パッケージの名前空間への代入はできない。 -例 > base::log # base パッケージ中の関数 log の情報を得る function (x, base = exp(1)) if (missing(base)) .Internal(log(x)) else .Internal(log(x, base)) <environment: namespace:base> > base::"+" # base パッケージ中の関数 "+" の情報を得る .Primitive("+") もしパッケージもしくは自己定義関数中で既存の関数と同じ名前のオブジェクトを重複定義しても、ダブルコロン演算子を使えば問題ない。インストール済みのパッケージをオンデマンドでロードするにも便利。 > base::log(10) [1] 2.302585 > log <- function(x) exp(x) # log 関数の再定義 > log function(x) exp(x) > log(1) # 再定義による値 [1] 2.718282 > base::log(1) # 本来の log 関数の値が得られる [1] 0 > x = rnorm(100) > medpolish(x) # eda パッケージ中の関数 medpolish (未ロードなのでエラー) Error: couldn't find function "medpolish" > eda::medpolish(x) # eda パッケージを自動ロードするのでエラーにならない Final: 0 Median Polish Results (Dataset: "x") Overall: -0.08048712 Row Effects: [1] 0.859951758 0.055019070 1.770125505 -0.700803046 0.792711958 # 以下省略 > eda::medpolish <- function(x) {median(x)} # eda パッケージ中への代入は不可 Error: Object "eda" not found ** cut 関数 ： 数値ベクトルを指定した区間に分類する (2003.12.05) [#b330a797] 値の分類区間を breaks ベクトルで与え、分類ラベルを labels で与える。ラベルは因子になるので、数値としてのあつかいはそのままでは不可能 > x <- 2000*runif(10) > y <- cut(x, breaks=c(-Inf, 500, 1000, 2000, Inf), labels=c(0,1,2,3)) > y [1] 2 2 2 2 2 2 2 2 1 1 Levels: 0 1 2 3 > as.integer(y) # 数値に変更 [1] 3 3 3 3 3 3 3 3 2 2 ## ラベルを FALSE にすると区間のかずに応じたラベルが数値として与えられる > y <- cut(x, breaks=c(-Inf, 500, 1000, 2000, Inf), labels=FALSE) > y # 結果は数値 0.1.2.3 [1] 3 3 3 3 3 3 3 3 2 2 ## labels 引数を省略すると該当区間を表す文字列がラベル（因子）として与えられる > y <- cut(x, breaks=c(-Inf, 500, 1000, 2000, Inf)) > y [1] (1e+03,2e+03] (1e+03,2e+03] (1e+03,2e+03] (1e+03,2e+03] (1e+03,2e+03] [6] (1e+03,2e+03] (1e+03,2e+03] (1e+03,2e+03] (500,1e+03] (500,1e+03] Levels: (-Inf,500] (500,1e+03] (1e+03,2e+03] (2e+03,Inf] ** 指定秒数プログラムを休止する Sys.sleep [#vca026cf] > p1 <- proc.time() # 現在の内部タイマーを得る > Sys.sleep(3.7) # 3.7 秒実行休止 > proc.time() - p1 # 経過時間を得る [1] 0.00 0.01 3.71 0.00 0.00 # 3.71 秒経過 ** あるディレクトリ中にあるファイルリストを得る (2004.2.10) [#o4f5f6ac] 書式： list.files(path = ".", pattern = NULL, all.files = FALSE, full.names = FALSE, recursive = FALSE) dir(path = ".", pattern = NULL, all.files = FALSE, # list.files の別名 full.names = FALSE, recursive = FALSE) 引数： -COLOR(red){path} 完全なパス名を与える文字列 -COLOR(red){pattern}: オプション。ファイル名を指定する正規表現 -COLOR(red){all.files}: 論理値。もし 'FALSE' なら隠しファイルは除く。もし 'TRUE' ならすべてのファイルが対象。 -COLOR(red){full.names}: 論理値。もし 'TRUE' ならディレクトリパスがファイル名の前に付け加えられる。もし 'FALSE' ならファイル名だけ。 -COLOR(red){recursive}: 論理値。サブディレクトリを再帰的にリストするか？ 返り値: ファイル名を表す文字列のベクトル~ ~ 注意：OSに依存する部分がある。 例： list.files(R.home()) ## 最初が文字 a-l と r (大文字を含む)だけをリストアップ(正規表現) dir("../..", pattern = "^[a-lr]", full.names=TRUE) ## 現ディレクトリ中のすべてのファイルについてある処理を行なう例 lapply(list.files(), function(x) paste('filename is ',x)) ## 同じことを for ループで実現 files <- list.files() for (each.file in files){ paste('filename is ',x) } ## 正規表現によるファイル名(a3b.txt といったファイル名)のリストアップと、 ## その内容を読み込んだデータファイルをリストへ一括格納。 tab <- lapply(list.files(pattern="^?[[:digit:]]?.txt"), read.table) # 各データフレームは tab[[1]],...,tab[[4]] で得られる ** オブジェクトに属性を指定 [#h92917a0] COLOR(red){structure(x, ...)} はオブジェクト COLOR(red){x} に COLOR(red){tag = value} の形で任意個数の属性を付与する。 > structure(1:6, dim = 2:3) # 次元属性 dim = 2:3 を与える [,1] [,2] [,3] [1,] 1 3 5 [2,] 2 4 6 ** コンソールへの出力の際、一行の幅を変える (r-help 記事より、2004.10.09) [#hed088e1] 例えばベクトルの内容を出力する際、一行あたりの要素の数を少なくしたいことがある。 width オプションの既定値は80文字である。 > x <- runif(10) > x [1] 0.3670568 0.7209854 0.9568679 0.4202735 0.1225961 0.2212387 0.3502483 [8] 0.1798455 0.8755473 0.5800494 > options(width=30) > x [1] 0.3670568 0.7209854 # 一行の文字数が30を越えないように折り返される [3] 0.9568679 0.4202735 [5] 0.1225961 0.2212387 [7] 0.3502483 0.1798455 [9] 0.8755473 0.5800494 > options(width=50) > x [1] 0.3670568 0.7209854 0.9568679 0.4202735 [5] 0.1225961 0.2212387 0.3502483 0.1798455 [9] 0.8755473 0.5800494 ** ファイル・ディレクトリを R から操作する低水準インタフェイス関数群 (2004.11.18) [#b9f6de1f] Unix-like OS では system() 関数を使い、Unix 命令を実行可能なので、これらの関数がなくてもとりあえず困らない。またこれらの関数がすべての OS でうまく動くかどうか不明。 OS 依存のディレクトリ、ファイル名の困難を避けるために、list.files() 関数(少し上を参照)と併用すると安心。 file.create(...) # ファイルを新規作成 file.exists(...) # ファイルが存在するかどうかチェック file.remove(...) # ファイルを削除 file.rename(from, to) # ファイル名を変更 file.append(file1, file2) # ファイルを合併 file.copy(from, to, overwrite = FALSE) # ファイルをコピー file.symlink(from, to) # シンボリックリンク dir.create(path, showWarnings = TRUE) # 新規ディレクトリを作成 **再帰呼び出し（2005.1.22） [#tf064f6b] ***再帰呼び出し [#z88a51c4] 以下のような入れ子になったリストを考えてみる． list(10, list(list(30, 20), 1:10), list(2:5)) ここで，このようなリストの末端にある数値を全て足す関数を定義してみる．このリストの個々の要素は同じ構造の（ただし入れ子が一段少ない）リストなので，個々の要素を自分自身を呼び出すことによって（再帰呼び出しによって）処理することが出来る． # リスト x の末端にある数値の総和を求める関数 myfunc <- function(x) { if (is.atomic(x)) return(sum(x)) # リストでなければ単に和を返す s <- 0 # 変数を初期化 for (i in x) { # リストの要素それぞれについて s <- s + myfunc(i) # 末端の数値の和を計算してそれを加算していく } return(s) # 総和を返す } 次は再帰プログラムの例ではありきたりの階乗を求めるプログラムである． myfunc <- function(n) { if (n <= 1) return(1) else n * myfunc(n-1) } myfunc(5) [1] 120 # 改良版：ベクトルの各要素ごとに階乗を求められる myfunc2 <- function(x) { recurse <- (x>1) x[!recurse] <- 1 if (any(recurse)) { y <- x[recurse] x[recurse] <- y*myfunc2(y-1) } return(x) } myfunc2(5) [1] 120 myfunc2(1:10) [1] 1 2 6 24 120 720 5040 40320 [9] 362880 3628800 Recall によって関数の名前に依存せずに再帰呼び出しを行なうことも出来る．先の例では関数名を変えると本体で自分自身を呼び出しているところも変更しなければらなかったが，Recall を使って再帰呼び出しを行なっていればその必要はなくなる． ***数値積分 [#y688af94] 再帰関数を用いることで，数値積分を行うことが出来る．以下は数値積分を行う関数 area(関数, 左端の点, 右端の点) である． area <- function(f, a, b, eps = 1.0e-06, lim = 10) { fun1 <- function(f, a, b, fa, fb, a0, eps, lim, fun) { ## 関数‘fun1’ は‘area’ の中からだけ見える d <- (a + b)/2 h <- (b - a)/4 fd <- f(d) a1 <- h * (fa + fd) a2 <- h * (fd + fb) if(abs(a0 - a1 - a2) < eps || lim == 0) return(a1 + a2) else { return(fun(f, a, d, fa, fd, a1, eps, lim - 1, fun) + fun(f, d, b, fd, fb, a2, eps, lim - 1, fun)) } } fa <- f(a) fb <- f(b) a0 <- ((fa + fb) * (b - a))/2 fun1(f, a, b, fa, fb, a0, eps, lim, fun1) } 上の関数 area() を定義すると，以下のように積分計算が出来る． area(sin,0,pi/2) [1] 0.999987
