R tips (最終改訂:2003/03/23)
(間瀬ウェッブページに戻る)


以下は R のメイリングリスト r-help や R-jp に投稿された記事や、個人的につまずいた経験から まとめた R 使用上のヒントを脈絡無くまとめたものです(内容の正確さは保証できません)。出典は 一々述べませんが、まとめて感謝の意を表します。

なお充実したPaul Johnson による R tips 集も参考になります。
R を含む統計一般に関する情報が青木繁伸さんと 中澤港さんのウェッブページから得られます。


目次グループ

[ベクトル、行列、配列、リスト], [グラフィックス],[プログラミング],[有用な関数]
[ヘルプ],[数値解析],[インタフェイス],[その他]

目次


    === ベクトル、行列、配列、リスト ===  
  1. (実際は整数だけからなる)実数型の配列・行列の保管モードを整数型にする
  2. 行列からその一部分を(ベクトルとして)取り出す
  3. 行(列)ベクトルをあたえて行列を作る
  4. 行列のベクトル化
  5. ベクトルの大小順並べ変え order(),rev()
  6. 複数回の同一試行の結果をリストに一括記録 lapply(), try()
  7. 行列の列(行)和からなるベクトルを計算
  8. ベクトルの成分がある性質をもつような添字を取り出す
  9. ベクトルの外積 outer()
  10. 行列の添字操作における結果のベクトル化を防ぐオプション drop=FALSE
  11. 配列のあるマージンに関数を適用 apply()
  12. 行列をある行・列に関してソート sort()
  13. 行列の特定要素の添字の取り出し which(., arr.ind=TRUE)
  14. 配列の添字順序の変更 aperm()
  15. (二重)リストの初期化 as.list()
  16. 同じ範囲に対する多重ループをリストを範囲に取ることにより簡略化
  17. 成分を与えて行列を作る matrix()
  18. パターンを持つベクトルの生成 sequence()
  19. ベクトルからある値に最も近い要素の添字を求める

  20. === グラフィックス ===  (目次先頭に戻る)
  21. 幹葉表示 stem 関数の scale パラメータ
  22. グラフの一部に影を付ける
  23. プロットの x 軸ラベルを45度回転
  24. 行列のグラフィカル表示 image()

  25. === プログラミング ===  (目次先頭に戻る)
  26. 自前の関数を保存し、次回に使えるようにする
  27. 数値ベクトルの対話的入力
  28. 空のオブジェクトを作る <- NULL, <- numeric(0)
  29. キー入力があるまでプログラム、プロットの実行を停止 readline(),par(ask=TRUE)
  30. 通し番号付きのオブジェクト名をプログラム中で生成 assign()
  31. 1<=i<j<=n に対する二重のループ
  32. 日付・時刻 date(), Sys.time()
  33. 条件分岐 ifelse()
  34. R の命令を表す文字列を評価実行 eval(paste(text=...))
  35. コード中の長い行の一括注釈化 if(0){...}
  36. 関数から複数の値を返す (リスト返り値)
  37. 関数内部での関数定義
  38. 再帰的関数定義
  39. 永続付値 <<- の秘密
  40. 一次的作業環境の使用 new.env(TRUE,NULL)

  41. === 有用な関数 ===  (目次先頭に戻る)
  42. ランダムサンプル、置換、ブートストラップ sample()

  43. === ヘルプ ===  (目次先頭に戻る)
  44. R のドキュメントやメイリングリストからキーワード検索
  45. R ヘルプ文章のキーワード検索関数 help.search()
  46. ヘルプ機能 help(),?,help.search(),help.start(),apropos()

  47. === 数値解析 ===  (目次先頭に戻る)
  48. 浮動小数演算の怪 machine epsilon
  49. 固有ベクトルの計算とチェック eigen()
  50. 正方行列の逆行列 solve()
  51. 行列の成分の自乗の総和
  52. 数値の表示桁数 options(digits=...)

  53. === インタフェイス ===  (目次先頭に戻る)
  54. R のバッチモード実行(Unix の場合)
  55. gzip を使って圧縮されたファイルの読み込み
  56. R の全コンソール入力・出力のログファイル化(Unix の tee 命令)
  57. 現在のディレクトリを R 中から変更 getwd(),setwd()

  58. === その他 ===  (目次先頭に戻る)
  59. 平面の Delaunay 三角化と Voronoi 分割用パッケージ
  60. 時間・時刻・日付関連 Sys.time(),proc.time(),system.time()
  61. パターンによるオブジェクト抹消 rm(list=pat("...")

  62. (目次先頭に戻る)

■ (実際は整数だけからなる)実数型の配列・行列の保管モードを整数型にする (目次に戻る)
8バイト実数から4バイト整数に変換(必要メモリが半分になる)
storage.mode(Array) <- "integer" 
または
Array[] <- as.integer(Array)
両者とも dim, dimnames 等の属性を保つ。 (注:Array 要素は実際は整数型に強制変換されている)
■ 行列からその一部分を(ベクトルとして)取り出す (目次に戻る)
> x <- matrix(runif(5*5),ncol=5,nrow=5)
> y <- x[cbind(c(1,2,3),c(5,4,3))]
y はベクトル c(x[1,5],x[2,4],x[3,3]) になる。

注:cbind(1:3,5:3) は実際は次の行列をあらわす

     [,1] [,2]
[1,]    1    5
[2,]    2    4
[3,]    3    3
特に x[cbind(1:5),cbind(1:5)] は x の対角成分からなるベクトル diag(x) と一致する。
■ R のドキュメントやメイリングリストからキーワード検索する方法 (目次に戻る)

(1) CRAN には R 関連のドキュメントや、メイリングリストの 過去記事中からキーワード検索するエンジンへのリンクがある。

(2) Google のキーワード検索を用い次の書式で検索(これはキーワード rpart と NA を含む該当文書を全てリストアップする)

site:r-project.org rpart NA

(3) インターネットダウンロードソフト wget を使いメイリングリスト r-help の過去のメイルを取得

書式: wget --mirror http://cran.r-project.org/doc/mail-archives/r-help

(注意:文章の総量は莫大でありサーバーに多大の負担をかけるので要注意)この中からキーワード検索するには grep を使うか 全文検索ソフト(例えば Namazu)をさらに使う。

(4) ペンシルバニア大学の心理学科の運営する R およびその他の関連文書の検索ツール htdigを使用


■ 行(列)ベクトルをあたえて行列を作る (目次に戻る)
>mat <- rbind(c(1,2,3), c(4,5,6)) # 行ベクトルをあたえて行列を作る
>mat
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
>cbind(c(1,2,3),c(4,5,6)) # 列ベクトルをあたえて行列を作る
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6

■ 行列のベクトル化 (目次に戻る)
>mat <- rbind(c(1,2,3), c(4,5,6))  #2x3 行列 
>mat
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
縦ベクトル t(c(1,2,3,4,5,6)) にするには
> vec <- as.vector(t(mat))
> vec
[1] 1 2 3 4 5 6
とする。行列は列主導で配置されるから転置演算 t(mat) が必要。転置しないと
> vec <- as.vector(mat)
> vec
[1] 1 4 2 5 3 6
もしベクトルでなく 1x6 行列が必要なら次のようにする: (行列は dim 属性 c(nrow,ncol) を持つベクトルに他ならないことを注意):
> matrix(t(mat), ncol=1)
     [,1]
[1,] 1
[2,] 2
[3,] 3
[4,] 4
[5,] 5
[6,] 6

■ ベクトルの大小順並べ変え order(),rev() (目次に戻る)
> x <- c(2,4,3,1)
> x
[1] 2 4 3 1
> ox <- order(x)
> ox
[1] 4 1 3 2
> x[ox]
[1] 1 2 3 4
逆順に並べ変えるには
> rox <- rev(order(x))
> rox
[1] 2 3 1 4
> x[rox]
[1] 4 3 2 1
長さの等しい二つのベクトル x,y を x の大小順にプロット:
> oo<-order(x)
> plot(x[oo],y[oo],type="l")

■ 浮動小数演算の怪 machine epsilon (目次に戻る)

(1) 計算機における浮動小数演算は有限桁の2進小数で行われるため数学的な計算とは異なる微妙な差異 (2進小数から10進法小数に最終的に変換する誤差を含め)が現われる。(1+x)-1 = x となる正の最小実数 x ((1-y)-1 = -y となる正の最小実数 y) が存在し machine epsilon (machine negative epsilon) と呼ぶ。IEEE 754 という浮動小数点演算の規約では x=2^(-52) (y=2^(-53)) となる。R の浮動小数点演算に関する基本変数 リスト .Machine にはこの値が収められている。

> .Machine$double.neg.eps
[1] 2.220446049250313080847e-16
> .Machine$double.neg.eps
[1] 1.110223024625156540424e-16
簡単に言えばこれらの数の半分以下の数を 1 より大きな数に加える・引く操作を行っても桁落ちで無効になる。
> options(digits=22)
> 1+(1:9)/10*2^{-52}
[4] 1.000000000000000000000 1.000000000000000000000 1.000000000000000222045
[7] 1.000000000000000222045 1.000000000000000222045 1.000000000000000222045
> 1-(1:9)/10*2^{-53}
[1] 1.0000000000000000000000 1.0000000000000000000000 1.0000000000000000000000
[4] 1.0000000000000000000000 1.0000000000000000000000 0.9999999999999998889777
[7] 0.9999999999999998889777 0.9999999999999998889777 0.9999999999999998889777
一連の演算結果が数学的には整数値(例えば 0)になるはずでも、詳しく見るとそうなっていないことがしばしば 起きる理由の一つが演算途中で machine epsilon 未満のオーダーの誤差が紛れ込むことにある。したがって、 たとえば実数 x に対し if (x == 0) という判定条件は危険であり、代わりに例えば if (abs(x) <= 5*.Machine$double.epsilon) 程度の判定条件にとどめた方が良い。これは収束の判定条件で特に問題となる。

(2) 実数の計算機における浮動小数点表現の問題から、大きさが極端に異なる実数を多数回加算すると時として 思いもかけぬ結果を得ることがある。 例:数列 1:10 をスチューデント化したものを3乗して加える。答えは正しく零になるはずであるが

[1] -4.440892e-16 
丁度 4*.Machine$double.neg.eps に等しいことに注意(但しこの例はOSやプラットフォーム依存らしく、正しく零になる こともある!)。この問題に対する根本的な解決策は多倍長演算を用いるしかないが、上の例のように意外に無邪気 そうな例でも起きかねないことは常に心すべきである。比較的簡単に実行できる対処法は絶対値が小さい順に加える ことである。例えば次のような"用心深い"和関数
my.sum <- function(x) {sum(x[order(abs(x))])}
を用意して使う。但しこれでも万全ではなく、結果が数学的に丁度ある数になるはずという思い込みは避けるのが安全。
■ 複数回の同一試行の結果をリストに一括記録 lapply(), try() (目次に戻る)
注:ある試行でエラーが起きても残りを続行するために try 関数を使用(help(try)参照)
> doit <- function() {mean(rnorm(30))}      # 30個の正規乱数の平均値を求める関数
> res <- lapply(1:3, function(i) try(doit()))
> res
[[1]]
[1] -0.359225966054957
[[2]]
[1] 0.0249291122975083
[[3]]
[1] -0.0972657631255351

■ 行列の列(行)和からなるベクトルを計算 (目次に戻る)
> 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
> drop(x %*% rep(1,ncol(x)))    # 行和 (drop は 1x4 行列をベクトル化する)
[1] 15 18 21 24
> drop(t(x) %*% rep(1,nrow(x))) # 列和
[1] 10 26 42
もう一つのやり方
> apply(x,2,sum)         # 列和
[1] 10 26 42
> apply(x,1,sum)         # 行和
[1] 15 18 21 24          

■ 幹葉表示 stem 関数の scale パラメータ (目次に戻る)
stem(z, x) で幹の長さが既定値(=1)の x 倍になる
> z <- rnorm(100)
> stem(z)

   The decimal point is at the |

    -2 | 310
    -1 | 8666555
    -1 | 4333110
    -0 | 99988888777777655555
    -0 | 4444443322221100
     0 | 00011111223333444444
     0 | 5666777899
     1 | 000233344
     1 | 5788
     2 | 001
     2 | 6
  
  > stem(z,10)

    The decimal point is 1 digit(s) to the left of the |
  
    -22 | 8
    -20 | 6
    -18 | 6
    -16 | 841
    -14 | 9528
    -12 | 8321
    -10 | 420
     -8 | 5364
     -6 | 7665443187
     -4 | 7887751
     -2 | 877551931
     -0 | 983731
      0 | 235667146
      2 | 582238
      4 | 0114436
      6 | 12602
      8 | 4676
     10 | 236
     12 | 90489
     14 | 9
     16 | 27
     18 | 3
     20 | 048
     22 | 
     24 | 9

■ R ヘルプ文章のキーワード検索関数 help.search("...") (目次に戻る)
  > help.search("%*%")  # 文字列 "%*%" を含むヘルプ項目から検索(関連項目が以下のように表示される)
  
  Help files with alias or title matching `%*%',  
  type `help(FOO, package = PKG)' to inspect entry `FOO(PKG) TITLE':
  
  Arithmetic(base)        Arithmetic Operators
  kronecker(base)         Kronecker products on arrays
  match(base)             Value Matching
  matmult(base)           Matrix Multiplication
  outer(base)             Outer Product of Arrays

  > help(mutmult)       # 項目 mutmult のヘルプを表示

  matmult                 package:base                 R Documentation
  
  Matrix Multiplication
  
  Description:
  
       Multiplies two matrices, if they are conformable. If one argument
       is a vector, it will be coerced to a either a row or column matrix
       to make the two arguments conformable. If both are vectors it will
       return the inner product.
  
  Usage:
  
       a %*% b
  
  ....... (以下省略) .......

■ ベクトルの成分がある性質をもつような添字を取り出す (目次に戻る)
  > x <- 12:1
  > x
   [1] 12 11 10  9  8  7  6  5  4  3  2  1
  > x[x>5]
  [1] 12 11 10  9  8  7  6
  > y <- (1:length(x))[x>5]        # y は x[i]>5 であるような添字のベクトル
  > y                
  [1] 1 2 3 4 5 6 7
  > x[y]
  [1] 12 11 10  9  8  7  6

■ 空のオブジェクトを作る <- NULL, <- numeric(0) (目次に戻る)
  > x <- NULL         # または x <- numeric(0)
  > x
  NULL                # x は空のオブジェクト
  > x[3] <- 3         # x[1], x[2] 抜きで x[3] に値を代入
  > x
  [1] NA NA 3         # x[1], x[2] はまだ不定
  > x[1:2] <- c(2,3)
  > x
  [1] 1 2 3

■ 自前の関数を保存し、次回に使えるようにする save(), source() (目次に戻る)
(1) 関数の定義をテキストファイル(例えば "myfun.R")に書き source("myfun.R") 関数で読み込む。 R の初期設定ファイル .Rprofile に行 source("myfun.R") を入れておけばいつも自動的に 読み込まれる。

(2) R セッション中に関数定義を save() 関数で Rdata ファイルに保存しておく。関数 attach() で 定義ファイルを読み込むことができる。

    
> myfun <- function(x) cat(x,"\n")
> save(myfun, file="myfun.Rdata")
   .......... 次の R セッション............
> attach("myfun.Rdata")

■ 数値ベクトルの対話的入力 readline() (目次に戻る)
readline 関数で指定プロンプトを提示し、対話的に文字列を入力(任意の分離記号を指定できる) 結果 Str は文字列 "10,20,30"。strsplit 関数で文字列リストに変換し、unlist 関数で文字列ベクトル化、 as.numeric 関数で数値化する
> Str <- readline("Enter x: ")               
Enter x: 10,20,30                            
> x <- as.numeric(unlist(strsplit(Str, ",")))  
> x
[1] 10 20 30                        # x は数値ベクトル c(10,20,30)

■ ベクトルの外積 outer() (目次に戻る)
z <- outer(x,y,FUN="afun") は次元 c(dim(x),dim(y)) は 要素 z[i,j] が afun(x[i],y[j]) である行列になる。 関数 afun は任意、ベクトルにタグ名があれば保存される。x%o%y は outer(x,y,FUN="*") の省略形
> x <- runif(3)
> y <- 1:3
> outer(x,y,FUN="^")
            [,1]        [,2]         [,3]
 [1,] 0.56709640 0.321598331 0.1823772568  # [1,2] 成分は 0.56709640^2
 [2,] 0.06480775 0.004200045 0.0002721955
 [3,] 0.06183690 0.003823802 0.0002364520

■ 行列の添字操作における結果のベクトル化を防ぐオプション drop=FALSE (目次に戻る)
> x_1:12
> dim(x) <- c(4,3)
> y <- x[c(2,4),]
> dim(x[c(2,4),])
[1] 2 3
> dim(x) <- c(12,1)           # x は 12x1 行列
> dim(x[c(2,4),])
NULL                          # 結果はもはや行列でない!
> dim(x[c(2,4),,drop=FALSE])  # オプション drop=FALSE
[1] 2 1                       # 結果は 2x1 の行列!

■ グラフの一部に影を付ける (目次に戻る)
例:標準正規分布の密度関数のグラフを -4<=x<=4 の範囲でプロット。そのうち -2<=x<=2 の範囲に灰色の影を 付ける(実際は多数の多角形に分割し polygon 関数で塗りつぶし)
> plot(dnorm,-4,4)
> xvals <- seq(-2,2,length=50)  # 領域をx軸方向に50個の多角形(台形)に等分割
> dvals <- dnorm(xvals)         # 対応するグラフの高さ   
> polygon(c(xvals,rev(xvals)),c(rep(0,50),rev(dvals)),col="gray")
   # xvals,rep(0,50) の組合せでx軸上の辺を結ぶ
   # rev(xvals),rev(dvals) の組合せでグラフに沿った辺を結ぶ

■ gzip を使って圧縮されたファイルの読み込み (目次に戻る)
(1) 関数 zip.file.extract()
(2) 関数 zip.unpack() (Windows版)
(3) 関数 shell() または system() を使い予め解凍した上で通常ファイルとして読み込み
(4) 関数 gzfile() を用いる。これは gzip 化された出力を行うことも出来る

■ プロットの x 軸ラベルを45度回転 (目次に戻る)
> plot(x, y, xaxt="n")
> axis(1, labels=FALSE)
> pr <- pretty(x)
> par(xpd=TRUE)
> text(pr, par("usr")[3], pr, adj=0, srt=-45)  # srt=-45 で時計回りに45度回転

■ 配列のあるマージンに関数を適用 apply() (目次に戻る)
> 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

■ キー入力があるまでプログラム、プロットの実行を停止 readline(), par(ask=TRUE) (目次に戻る)
キー入力があるまでプログラム、プロットの実行を停止。指定メッセージを出力し、何かをキー入力したら以下を実行
> 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) を提示する。 これは example 関数によるデモを実行する際、ひき続くグラフィックス表示を 一つずつ眺めるためにも有効である。par(ask=FALSE) を実行すれば元に戻る。
■ 通し番号付きのオブジェクト名をプログラム中で生成 assign() (目次に戻る)
文字ベクトルや数値を用いてあるオブジェクト名をプログラム中で生成し、値を付値するには
> assign(paste("file", 1, "max", sep=""), 1)
> ls()
[1] "file1max"
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)) 
  }
}

■ 1<=i<j<<=n に対する二重のループ (目次に戻る)
(1) for(i in 1:(n-1))
      for(j in (i+1):n))           # 但し n=1 で失敗
         ...
(2) for( i in seq(n) )
       for( j in seq(n)[-seq(i)] ) # n=1 でもOK
         ...

■ 行列のグラフィカル表示 image() (目次に戻る)
行列を画像と考え image 関数で表示する。option の col で色調を変えることが出来る (col の例: rainbow,heat.colors,topo.colors,terrain.colors, gray)
> x <- y <- 1:10
> z <- matrix(rnorm(100),ncol=10)
> image(x,y,z, col=topo(12))      # 地形図色調を12段階で使用

■ 行列をある行・列に関してソート order() (目次に戻る)
> x_matrix(runif(9),ncol=3)
> x
          [,1]      [,2]      [,3]
[1,] 0.6999342 0.4610615 0.3838415
[2,] 0.5060148 0.9964767 0.5274574
[3,] 0.1510918 0.9524839 0.5363432
> order(x[,2])
[1] 1 3 2
> x[,order(x[,2])]                 # 第2列でソート
          [,1]      [,2]      [,3]
[1,] 0.6999342 0.3838415 0.4610615
[2,] 0.5060148 0.5274574 0.9964767
[3,] 0.1510918 0.5363432 0.9524839
> order(x[1,])
[1] 3 2 1
          [,1]      [,2]      [,3]
[1,] 0.1510918 0.9524839 0.5363432
[2,] 0.5060148 0.9964767 0.5274574
[3,] 0.6999342 0.4610615 0.3838415

■ 平面の Delaunay 三角化と Voronoi 分割用パッケージ (目次に戻る)
CRAN にある R パッケージ ``deldir''。statlib にある Splus 用のパッケージ ``delaunay'' を R に移植したもの。
■ R のバッチモード実行(Unix の場合) (目次に戻る)
(1) R プログラムファイル myfile.R のバッチ実行(結果はファイル myfile.Rout に記録される)
R BATCH myfile.R & 
(2) 次のような Unix のシェルスクリプトを使う(csh の場合)
#!/csh -f
#
# Something for preparation
#
R --vanilla << EOF
   # Put R Code here, it's convinient that you can use macro 
   # substitutions here.
EOF
#
# Something for finalize
#

■ 行列の特定要素の添字の取り出し which(., arr.ind=TRUE) (目次に戻る)
ある条件を満たす行列の特定要素の行列添字の全部を取り出すには which 関数を オプション arr.ind=TRUE で使用する。既定の arr.ind=FALSE では行列をベクトル化 した上で該当要素の添字ベクトルを返す。
> x <- matrix(1:12,3,4)
> x
   [,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> y_which(x%%3==0, arr.ind=TRUE)  ## 3 の倍数である要素の行列添字を取り出す 
> y
     row col
[1,]   3   1
[2,]   3   2
[3,]   3   3
[4,]   3   4
> x[y[1,1],y[1,2]]
[1] 3
> x[y[4,1],y[4,2]]
[1] 12
単に該当要素を全部取り出すならベクトル化した方がわかりやすい。R の行列・ 配列は次元属性を持つベクトルであり、ベクトルとしてアクセスできる。
> x[x%%3==0]          # 3の倍数である x の全要素 
[1]  3  6  9 12 
> x[x%%2==0] <- 0     # 行列 x の偶数要素を全て 0 に置き換え(行列をベクトルとして添字操作)
> x
   [,1] [,2] [,3] [,4]
[1,]  1  0  7   0
[2,]  0  5  0  11
[3,]  3  0  9   0

■ コード中の長い行の一括注釈化 if(0){...} (目次に戻る)
R コード中の # はそれ以降行末までをコメントとみなし無視する(日本語もOK)。 コード中の各行は先頭に # をつければ注釈できるが、長い行をひとかたまりに注釈化するには次の構文で囲む。 つまり if 文の条件が常に FALSE だから括弧内はすべて無視される
if(0) {(長い行........)}

■ 配列の添字順序の変更 aperm() (目次に戻る)
x <- array(runif(5^4),c(5,5,5,5))	    # 三次元配列
y <- aperm(x,perm=c(2,1,3,4))     	    # 添字順序の変更
z <- aperm(x,perm=match(1:4,c(2,1,3,4)))   # 添字順序の変更
zz <- aperm(y,perm=c(2,1,3,4))             # 添字順序の変更
とすると x[i1,i2,i3,i4] は y[i2,i1,i3,i4] となる。z[i1,i2,i3,i4] は x[i2,i1,i3,i4] となる。 zz と z は同じ。

perm 引数を省略すると添字順序が逆順になる。つまり perm=rev(seq(length(dim(x)))) としたのと同じ。

x <- matrix(runif(5^2), ncol=5)
y <- aperm(x,c(2,1))
と置くと y は t(x) と一致。
■ (二重)リストの初期化 <- as.list() (目次に戻る)
リストは属性の異なる任意のオブジェクトを一括保持でき、添字で操作できるため便利なデータタイプである。 リスト成分は二重かぎ括弧演算子 [[,]] でアクセスする。名前付き成分の場合には [["name"]] もしくは $ 演算子でアクセスする。Lst[[1]] がベクトル・配列の時は Lst[[1]][1,2] 等でその要素にアクセスできる。

一重リストの初期化は

> Lst <- as.list(NA)  # リストオブジェクト Lst の初期化
> Lst
[[1]]
[1] NA
> Lst[[1]] <- c(1,2,3) # リストの第一成分(ベクトル)
> Lst
[[1]]
[1] 1 2 3
> Lst[[3]] <- matrix(1:4,ncol=2)  # リストの第3成分(行列)
> Lst                            
[[1]]
[1] 1 2 3
[[2]]                             # 自動的に途中の成分が生成される
NULL
[[3]]
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> Lst[[2]] <- "second character element" # 空の第二成分に文字列を代入
> Lst
[[1]]
[1] 1 2 3
[[2]]
[1] "second character element"
     [,1] [,2]
[1,]    1    3
[2,]    2    4
最初から要素数を指定して空のリストを作るには
> Lst <- as.list(rep(NA,2))  # 要素数2
> Lst
[[1]]
[1] NA
[[2]]
[1] NA
リストを具体的要素を与えて構成するには list 関数を使う。各要素には名前タグを付けることが出来、 要素番号以外に名前タグで操作(二重かぎ括弧以外に $ 演算子を使える)できる。
> Lst <- list(sex=c(1,2),name=c("Tarou","Hanako"))
> Lst
$sex
[1] 1 2
$name
[1] "Tarou"  "Hanako"
> Lst[[1]]
[1] 1 2
> Lst[["sex"]]
[1] 1 2
> Lst$sex
[1] 1 2
二重リスト(リストの各成分がまたリスト)を初期化するには
> Lst <- as.list(c(NA,NA))
> Lst[[1]] <- as.list(c(NA,NA))
> Lst[[2]] <- as.list(c(NA,NA))
> Lst
[[1]]
[[1]][[1]]
[1] NA
[[1]][[2]]
[1] NA
[[2]]
[[2]][[1]]
[1] NA
[[2]][[2]]
[1] NA
> Lst[[1]][[2]] <- c(1,3)
> Lst[[1]]
[[1]]
[1] NA
[[2]]
[1] 1 3
> Lst[[1]][[2]]  # Lst[[1,2]] は使えない
[1] 1 3

■ 同じ範囲に対する多重ループをリストを範囲に取ることにより簡略化 (目次に戻る)
for ループのループ範囲にはベクトル以外にリストを与えることが出来る。 例えば2重ループ
> for (i in 1:n) { for (j in 1:n) f(i,j)}
> rangelist <- multiplerange(1:n)
> for (i in rangelist) f(i[1],i[2])
と出来る(同じ多重範囲のループが瀕出する際は、コードが簡潔になり、実行時間も 多少早くなる)。ここで multiplerange(1:n) は全ての組合せ c(i,j), 1<=i,j<=n, からなるリストを作る関数:
multiplerange <- function(range) {
  mplerange <- as.list(rep(NA,length(range)^2))
  ct <- 0 
  for (i1 in range) { for (i2 in range) { 
    mplerange[[ct<-ct+1]] <- c(i1,i2)
  }} 
  return(mplerange)
}

■ 日付・時刻 date(), Sys.time() (目次に戻る)
 
> 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年元旦からの経過秒数

■ R の全コンソール入力・出力のログファイル化(Unix の tee 命令) (目次に戻る)
R 起動時にコンソールへのキー入力、命令等の出力を全て一つのファイルに落しておき後で点検・編集することは 有用である。Unix ベースのOSでは R 起動時に
R | tee ログファイル名  (例 R | tee R.20030321.log)
とすればログファイル(名前はファイル名として可能なものであれば任意)にコンソールに表示されたものが全て 保存される。Unix の tee 命令はコンソールへの出力を分岐してファイルへも出力するためのツールである。

既存のファイルへの追加出力数には tee 命令にオプション -a を加える

R | tee -a ログファイル名  (例 R | tee -a R.20030321.log)

■ ランダムサンプル、置換、ブートストラップ sample() (目次に戻る)
長さ n の与えられたベクトルの要素から長さ m の部分ベクトルをランダムに取り出す
> x <- 1:n                  # 1:n の代わりに単に n としても良い
> sample(x,m,replace=TRUE)  # 同じ要素が選ばれても良い
> sample(x,m,replace=FALSE) # 同じ要素は選ばれ無い(既定動作)

> sample(x)                 # 特に m=n とするとランダムな置換になる

> sample(c(0,1), 100, replace = TRUE) # 100個のベルヌイ試行
オプションの確率ベクトルを与えると各要素が選ばれる確率を指定できる(既定値は等確率 1/n)
> p <- runif(n); p <- p/sum(p)
> sample(x,m,replace=TRUE,prob=p)

■ 成分を与えて行列を作る matix() (目次に戻る)
成分ベクトル x を行列に変換するには
> m <- matrix(x,nrow=i,ncol=j,byrow=TRUE, dirnames=z) # 一般形 
> M <- matrix(1:6,,3,byrow=TRUE) # ncol=3 (自動的に nrow=2 になる)
> M
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
> M <- matrix(1:6,,3)        # matrix(1:6,,3,byrow=FALSE) と同じ
> M
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> M <- matrix(1:6,3,)        # nrow=3 (自動的に ncol=2 になる)
> M
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6
     [,1] [,2]
[1,]    1    1               # ベクトル 1:3 がリサイクル使用されている
[2,]    2    2
[3,]    3    3


■ パターンを持つベクトルの生成 sequence() (目次に戻る)
ベクトル x に対し sequence(x) は seq(x[1]),seq(x[2]),... をこの順に並べたベクトルを生成する。
> sequence(3:1)
[1] 1 2 3 1 2 1
> sequence(1:3)
[1] 1 1 2 1 2 3
> sequence((-1):1)
> sequence(c(3,3,2))
[1] 1 2 3 1 2 3 1 2

■ ベクトルからある値に最も近い要素の添字を求める (目次に戻る)
> x=rnorm(1000)                           # 1000個の正規乱数
> which(abs(x-1.1) == min(abs(x-1.1)))    # 1.1 に最も近い値を持つ要素の添字
[1] 657
> x[which(abs(x-1.1) == min(abs(x-1.1)))] # 該当要素値
[1] 1.112788

■ 条件分岐 ifelse() (目次に戻る)
論理値ベクトル 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

■ ヘルプ機能 help(),?,help.search(),help.start(),apropos() (目次に戻る)
R の機能のヘルプ文章は help("...")、または ?... でオンライン表示される。
例: help("plot")、または ?plot

help.start() は Netscape を起動(もし使用可能なら)し、以後のヘルプ文章はブラウザに表示される。

help.search("...") はロードされている全てのパッケージから該当項目を検索する。
例: help.search("plot")

help.search は正規表現を受け付ける。
例: help.search("*trans*") (途中に trans を含む全ての項目)

検索すべきキーワードが曖昧な場合は apropos 関数にキーワードの一部を与えて検索する。
例: apropos("files") (タイトル・名前に文字列 files を含む項目を表示)


■ 時間・時刻・日付関連 Sys.time(),proc.time(),system.time() (目次に戻る)
> Sys.time()                           # OS の時間情報
[1] "2003-03-22 13:49:24 JST"
> ? Sys.time
> format(Sys.time(), "%a %b %d %X %Y") # 整形済みの時間情報
[1] "土  3月 22 13時49分42秒 2003"
> date()
[1] "Sat Mar 22 13:53:57 2003"

> ptm <- proc.time()                   # R 起動後の経過時間
> for (i in 1:50) mad(runif(500))
> proc.time() - ptm                    # 上記コードの実行時間
[1] 0.05 0.00 0.05 0.00 0.00

> system.time(for (i in 1:50) mad(runif(500))) # コードの実行時間
[1] 0.04 0.01 0.05 0.00 0.00
時間関数の返すベクトルは ユーザーCPUタイム、システムCPUタイム、経過時間、 副プロセスの使用ユーザーCPUタイム、副プロセスの使用システムCPUタイム (最後の二つは普通ゼロ)

時間の解像度はシステム固有であり、普通 1/100 秒。経過時間は最も近い時間に 1/100 秒単位で丸められる。


■ 固有ベクトルの計算とチェック eigen() (目次に戻る)
m は次のような対称行列:
 0.4015427  0.08903581 -0.2304132
 0.08903581 1.60844812  0.9061157
-0.23041322 0.9061157   2.9692562

> sm <- eigen(m, sym=TRUE)          # 関数 eigen を適用        
> sm$values                 
[1] 3.4311626 1.1970027 0.3510817   # 三つの固有値

>sm$vectors                         # それぞれに対応する固有ベクトル
     [,1]        [,2]       [,3] 
[1,] -0.05508142 -0.2204659  0.9738382 
[2,]  0.44231784 -0.8797867 -0.1741557 
[3,]  0.89516533  0.4211533  0.1459759

> V <- sm$vectors 
> t(V) %*% V                          # (計算誤差内で) V は直交行列
     [,1]          [,2]          [,3] 
[1,]  1.000000e+00 -1.665335e-16 -5.551115e-17 
[2,] -1.665335e-16  1.000000e+00  2.428613e-16 
[3,] -5.551115e-17  2.428613e-16  1.000000e+00

> V %*% diag(sm$values) %*% t(V)      # (計算誤差内で) m が復元されている
      [,1]       [,2]       [,3] 
[1,]  0.40154270 0.08903581 -0.2304132 
[2,]  0.08903581 1.60844812  0.9061157 
[3,] -0.23041320 0.90611570  2.9692562

■ R の命令を表す文字列を評価実行 eval(paste(text=...)) (目次に戻る)
String2Eval を R の命令を表す文字列とすると eval(parse(text = String2Eval)) はそれを解釈実行する。
> eval(parse(text="ls()"))
> eval(parse(text = "x[3] <- 5"))

■ 関数から複数の値を返す (リスト返り値) (目次に戻る)
一つの関数から複数の値を返すにはベクトル、配列等のオブジェクトを返り値とすれば良い。return 関数に 複数の引数を与えると、それらは自動的にリストとして返される。リストの各成分には元の変数名が名前 タグとして自動的に付加される。
> 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
[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

■ 関数内部での関数定義 (目次に戻る)
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)
  }

■ 正方行列の逆行列 solve() (目次に戻る)
正方行列の逆行列を solve() 関数で求めることができる。しかしながら逆行列を求めることは数値精度の点で 問題を生じることが多く、また時間を食うという意味で"本当に必要なときだけ"使うことが数値計算の常識 (実際逆行列そのものが必要になることは殆んど無いはず)。
> A=array(runif(9),c(3,3))
> B=solve(A)
> A%*%B                                         # 検査(数値精度の範囲内で単位行列
             [,1]          [,2]          [,3]
[1,] 1.000000e+00 -1.103176e-17 -2.745742e-17 
[2,] 1.083863e-16  1.000000e+00 -1.039106e-16
[3,] 3.718813e-17  1.237075e-16  1.000000e-00
> B%*%A                                         # 検査(数値精度の範囲内で単位行列
              [,1]          [,2]          [,3]
[1,]  1.000000e+00  5.573382e-17 -2.075354e-16
[2,] -2.778268e-17  1.000000e+00  2.326698e-16
[3,] -7.285839e-17 -5.995638e-17  1.000000e+00

■ 再帰的関数定義 (目次に戻る)
例: 階乗関数の定義
> test <- function(x){ifelse(x==0,1,x*test(x-1))}  # 関数 test の定義中で test を呼び出している。
> test(0)
[1] 1
> test(1)
[1] 1
> test(10)
[1] 3628800
> test(26)
[1] 4.032915e+26
> test(27)
Error in rep(yes, length = length(ans)) : evaluation is nested too deeply: infinite recursion?

■ 行列の成分の自乗の総和 (目次に戻る)
> x <- array(runif(8),c(2,4))
> sum(diag(t(x)%*%x))
[1] 3.185785
> sum(diag(x%*%t(x)))
[1] 3.185785

■ パターンによるオブジェクト抹消 rm(list=pat("...")) (目次に戻る)
不要な変数・関数・オブジェクトを作業スペースから抹消するにはパターン照合が使える。パターン照合には 正規表現や部分照合が使える。
> rm(list=objects(pattern="^lm.*"))   # lm. で始まる全てのオブジェクトを抹消
> rm(list=ls(pat="g"))                # オブジェクト名に g を含む全てのオブジェクトを抹消
ls(pat="...") 中には次のような正規表現が使える:
ls(pat="^foo")  名前が foo で始まるオブジェクト
ls(pat="foo$")  名前が foo で終るオブジェクト
ls(pat="foo")   オブジェクト名の一部に文字 foo を含む
注意: ワイルドカード * は(何故か)使えない。ls(pat="*"), ls(pat="g*") 等は全てのオブジェクトに該当する ので rm(list=ls(pat="g*") は全てのオブジェクトを抹消するので危険。
> ls()
[1] "f" "x" "y" "z"
>> ls(pat="*")        # 隠しファイル以外全てを表示
[1] "f" "x" "y" "z"
> ls(pat="*x")
[1] "x"
> ls(pat="x*")
[1] "f" "x" "y" "z"   # 隠しファイル以外全てを表示
ドット . で始まるファイルは隠しファイルとみなされ ls 命令では表示されない。特に表示するためには
> ls(all=TRUE)
[1] ".Op" ".Random.seed" ".Traceback" "f" "x" "y" "z" "x.temp"
とする。この種のオブジェクトは R 自身が自動的に生成するものがあり、無闇に抹消するのは危険。

関数中での使用には注意がいる。関数中で使用された ls 関数はその関数の実行環境中のオブジェクトのみを 探すので、もし関数を実行した環境中で探すにはオプション envir= で該当環境を指定する必要。

rm(list = ls(pat = "^g", envir = .GlobalEnv), envir = .GlobalEnv)  # 大域的環境中で探し、抹消する

■ 現在のディレクトリを R 中から変更 getwd(),setwd() (目次に戻る)
R からファイルを操作するには完全なディレクトリを込めて指定するのが一番安全だが、現在のディレクトリを R から変更することもできる。但し R から(Unix の)OS命令を実行する基本である system("命令名")は必ずしも希望 の結果を得られないこともあるらしい。
> getwd()             # 現在のディレクトリを表示
[1] "/home/mase"
> system("cd ..")     # Unix 命令 "cd .." を実行
> getwd()
[1] "/home/mase"      # しかし実際は変更はされていない
> getwd()
[1] "/home/mase/"
> setwd("..")
> getwd()
[1] "/home"           # 変更されている
> setwd("~")
> getwd()
[1] "/home/mase"

■ 数値の表示桁数 options(digits=...) (目次に戻る)
R の浮動小数の表示桁数 num を変えるには options(digits=num) を使う。それ以降の全ての表示に影響。 但し num は 22 を越えることはできない。しかし 16 ビットの浮動小数演算で実際に保証される桁数は おおよそ15桁であり、num=22 にしてもそれだけ正確になるわけではない。print() 関数も桁数指定オプ ションを持つ。
> options(digits=10); log(3.1)
[1] 1.131402111
> options(digits=15); log(3.1)
[1] 1.1314021114911
> options(digits=22); log(3.1)
[1] 1.131402111491100592744
> options(digits=23)
参考:
> .Machine$integer.max           # 最大可能正浮動小数 (=2^(31))
[1] 2147483647
> .Machine$double.xmax
[1] 1.797693134862315708145e+308 # 最大可能正整数

■ 永続付値 <<- の秘密 (目次に戻る)
永続付値演算子はRの大局的環境(.GlobalEnv という名前の隠しオブジェクト)中のオブジェクトに対する操作であり、 特に関数中で用いるときは注意がいる。一方で大局的環境中のオブジェクトは任意の関数から直ちに参照でき、特に 副関数に複雑な引数を引き渡すとき便利である。永続付値演算子を関数中で用いるとき陥りやす間違いは次の例でみられる:
> v <- 1
> test <- function () {v <- 2; v <<- 3; cat("v=",v,"\n")}
> test()
v= 2
> v
[1] 3
この例では関数内部での変数 v と永続付値で生成された大局的環境(つまり関数 test を実行した環境 .GlobalEnv 中の変数 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 の値(いつの間にかベクトルになっていることに注意)
教訓: 大局的変数は便利だが使用には注意がいる。使うときはできるだけ個性的な名前を用いる。決して大局的環境、 もしくは関数中の局所的変数と紛らわしい名前は使わない!
■ 一次的作業環境の使用 new.env(TRUE,NULL) (目次に戻る)
永続付値の項で述べた同じ名前でも環境(簡単にいえばメモリースペース)が異なれば違った実体を 持つという点は次の積極的に仮の環境を使う例で一層はっきりする。多少面倒だが多数の関数が 共有する変数はこうした一次的作業環境に置いておくと混乱が無いかも知れない。
> 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
(間瀬ウェッブページに戻る)