知っているといつか役にたつ(?)関数達 (2)

とりあえず使いそうもないが、知っているといつか役にたつかも知れない関数達(およびその使いかた)のメモです。r-help を見ていると「へぇー、こんなのあったの」という関数が毎週のように見つかりますね。メモしておかないとすぐ忘れそうなので、という理由で設けた極めて個人的なページですが、皆さんの「へぇー」も勝手に付け加えて下さい。

アーカイブ (1) 知っているといつか役に立つ(?)関数達


(155) LaTex? ソース中にRコマンド列をメモしておく方法(2008.06.04)

Rの出力を含む LaTeX 原稿中に対応するRコマンド(特に画像作成出力)をメモとして残しておきたいことがよくある。もちろん LaTeX の行コメントマーク % を各コマンド行の先頭に置いておけばよい。しかし、これでは再実行の際、マウスでC&Pするのにとても不便。「Sweave 使ったら」という声が聞こえそうだが、個人的には正直あんなもの使う気にはならぬ。永らく妙案も思いつかなかったが、最近やっと気づいた方法を紹介する(あまりに単純な方法なので FAQ だったらご免)。要するにRコマンド列を LaTeX の(無意味な)マクロ定義の中に置いておくだけ。

\def\temp{ % 次の行から R コマンドを並べる。\temp はマクロ名で適当 
  x  <- as.matrix(mtcars)
  rc <- rainbow(nrow(x), start=0, end=.3)
  cc <- rainbow(ncol(x), start=0, end=.3)
  hv <- heatmap(x, col = cm.colors(256), scale="column",
                RowSideColors = rc, ColSideColors = cc, margins=c(5,10),
                xlab = "specification variables", ylab= "Car Models",
                main = "heatmap(<Mtcars data>, ..., scale = \"column\")")
  dev.ocpy2eps(file="heatmap.eps")
  dev.off()
} % 忘れず括弧を閉じておく 
  • 私は以下のようなマクロを定義しています(常用する自分用のスタイルファイルに書いておけばいちいち定義しなくともよいので)。R 等で言えば if(0) { ... } のようなものですね。これは結構便利です。せっかく書いた文章の一部分を取りあえず除外したり,エラーが出るけどそれがどこで出ているか範囲を絞るためとか。お試し下さいませ。
    \newcommand{\SKIP}[1]{}
    使用法は
    \SKIP{
      スキップすべき範囲
    }
  • やはり同じようなことを考える人がすでにいましたか。ところで LaTeX はマクロ定義の構文的正しさを、実際に使われるまではチェックしないので(マクロたる所以)こうした方法が有効なわけです。R の if(0){あるメモ} とか無味な関数定義化 temp <- function(){あるメモ} では構文的な正当性を詮索されてしまいますから、LaTeX ほど融通が効かないわけです。

(154) ベクトルからNA値を除いたものの個数を知る方法(初級Q&Aより、2008.02.27)

色々考えられるが、とりあえず4種類。じつは complete.cases(), na.omit() 関数をベクトルに適用することなど、今まで考えもしなかった。

> ( x <- sample(c(1:3,NA), 10, rep=TRUE) )
 [1]  1  2  1  1 NA NA  2  3  2  2
> !is.na(x)
 [1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE
> sum(!is.na(x))
[1] 8
 
> x > -Inf
 [1] TRUE TRUE TRUE TRUE   NA   NA TRUE TRUE TRUE TRUE
> sum(x > -Inf, na.rm=TRUE)
[1] 8

complete.cases() 関数は行列、データフレームに適用すると、各行(ケース)がNA値を含まないか否かを示す論理値ベクトルを返すが、ベクトルにも使える。

> complete.cases(x)
 [1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE
> sum(complete.cases(x))
[1] 8

おそらく、以下が正統派(?)。

> na.omit(x)
[1] 1 2 1 1 2 3 2 2
attr(,"na.action")
[1] 5 6
attr(,"class")
[1] "omit"
> length(na.omit(x))
[1] 8

NA 以外の例外値がもしあったらどうなるか(心配性)。

> is.na(c(Inf,-Inf,NaN, NA))          # NaN (Not a Number)はNA扱い
[1] FALSE FALSE  TRUE  TRUE
> complete.cases(c(Inf,-Inf,NaN,NA)) # Inf, -Inf は正当な数扱い
[1]  TRUE  TRUE FALSE FALSE
> c(Inf,-Inf,NaN,NA) > -Inf      # -Inf > -Inf は FALSE
[1]  TRUE FALSE    NA    NA
> na.omit(c(Inf,-Inf,NaN,NA)) 
[1]  Inf -Inf
attr(,"na.action")
[1] 3 4
attr(,"class")
[1] "omit"

おまけ。

> x
     [,1] [,2] [,3] [,4]
[1,]    1   NA    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
> na.omit(x)              # NA値を含まない行だけからなる副行列が得られる
     [,1] [,2] [,3] [,4]
[1,]    2    6   10   14
[2,]    3    7   11   15
[3,]    4    8   12   16
attr(,"na.action")
[1] 1
attr(,"class")
[1] "omit"

(153) ブロック対角行列を作る(初級Q&Aより、2008.02.22)

(小)正方行列からなるリストを元にブロック対角行列を作るには Matrix パッケージの関数 bdiag() が使える。但しこれは独自の疎(スパース、ゼロ成分を縮約)行列形式だから、通常のマトリックス形式にするには as.matrix() 関数で変換する(気持ち悪くなければ必ずしも変換する必要は無いが)。

> L <- list(matrix(seq(4),2,2),matrix(seq(9),3,3),matrix(seq(25),5,5))
> library(Matrix)
> bdiag(L)
10 x 10 sparse Matrix of class "dgCMatrix"
                             
 [1,] 1 3 . . . .  .  .  .  .
 [2,] 2 4 . . . .  .  .  .  .
 [3,] . . 1 4 7 .  .  .  .  .
 [4,] . . 2 5 8 .  .  .  .  .
 [5,] . . 3 6 9 .  .  .  .  .
 [6,] . . . . . 1  6 11 16 21
 [7,] . . . . . 2  7 12 17 22
 [8,] . . . . . 3  8 13 18 23
 [9,] . . . . . 4  9 14 19 24
[10,] . . . . . 5 10 15 20 25
> as.matrix(bdiag(L))
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    3    0    0    0    0    0    0    0     0
 [2,]    2    4    0    0    0    0    0    0    0     0
 [3,]    0    0    1    4    7    0    0    0    0     0
 [4,]    0    0    2    5    8    0    0    0    0     0
 [5,]    0    0    3    6    9    0    0    0    0     0
 [6,]    0    0    0    0    0    1    6   11   16    21
 [7,]    0    0    0    0    0    2    7   12   17    22
 [8,]    0    0    0    0    0    3    8   13   18    23
 [9,]    0    0    0    0    0    4    9   14   19    24
[10,]    0    0    0    0    0    5   10   15   20    25

小行列は必ずしも正方行列である必要はない。

> X <- list(matrix(seq(4),2,2), matrix(seq(6),2,3), matrix(seq(24),3,8))
> (Y <- bdiag(X))
7 x 13 sparse Matrix of class "dgCMatrix"
                                   
[1,] 1 3 . . . . . .  .  .  .  .  .
[2,] 2 4 . . . . . .  .  .  .  .  .
[3,] . . 1 3 5 . . .  .  .  .  .  .
[4,] . . 2 4 6 . . .  .  .  .  .  .
[5,] . . . . . 1 4 7 10 13 16 19 22
[6,] . . . . . 2 5 8 11 14 17 20 23
[7,] . . . . . 3 6 9 12 15 18 21 24

疎行列は普通の行列と同じと思って操作してよい。結果のクラスはケースバイケースで様々。

> Y %*% t(Y)
7 x 7 sparse Matrix of class "dgCMatrix"
                               
[1,] 10 14  .  .    .    .    .
[2,] 14 20  .  .    .    .    .
[3,]  .  . 35 44    .    .    .
[4,]  .  . 44 56    .    .    .
[5,]  .  .  .  . 1436 1528 1620
[6,]  .  .  .  . 1528 1628 1728
[7,]  .  .  .  . 1620 1728 1836
> Y %*% t(as.matrix(Y))
7 x 7 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]   10   14    0    0    0    0    0
[2,]   14   20    0    0    0    0    0
[3,]    0    0   35   44    0    0    0
[4,]    0    0   44   56    0    0    0
[5,]    0    0    0    0 1436 1528 1620
[6,]    0    0    0    0 1528 1628 1728
[7,]    0    0    0    0 1620 1728 1836

疎行列形式だからといって保管メモリサイズが小さくなるかというと、必ずしもそうではないので得したつもりになってはいけない。

> object.size(X)
[1] 552
> object.size(Y)
[1] 1244
> object.size(as.matrix(Y))
[1] 928

(152)列名ラベル付きの行列の添字操作で列範囲を指定するのに「:」記法を使う(2008.02.22)

データフレームの加工に便利な subset() 関数には変数名ラベルに対する「:」記法による範囲指定ができることは結構知られているが、同じことは列名ラベルが付いた行列でもできる(?subset とすればしっかり書いてあることだが)。

> X <- data.frame(abc=runif(3), xyz=runif(3), foo=1:3, bar=rnorm(3), ouch=10:12)
> XX <- as.matrix(X)
> XX
           abc        xyz foo        bar ouch
[1,] 0.1786102 0.01356814   1  0.9822413   10
[2,] 0.2946404 0.40075416   2 -1.5008734   11
[3,] 0.6441781 0.44519585   3 -0.6949562   12
> subset(XX, select=xyz:bar)   # XX[,2:4] とすれば良いわけではあるが
            xyz foo        bar
[1,] 0.01356814   1  0.9822413
[2,] 0.40075416   2 -1.5008734
[3,] 0.44519585   3 -0.6949562

(151) 関数のベクトル化(2008.02.15)

Vectorize関数はベクトル化されていない関数をベクトル化する。下の(126)を参照。

> foo <- function(x) 1   # 定数値関数
> foo(1:5)               # foo はベクトル化されていないから返り値は一つ  
[1] 1
> Foo <- Vectorize(foo)  # foo をベクトル化した関数 Foo
> Foo(1:5)               # ベクトル化定数値関数
[1] 1 1 1 1 1
> mapply(foo, 1:5)       # Vectorize は実は mapply 関数のラッパー関数
[1] 1 1 1 1 1

変数の一部だけベクトル化することもできる。

> bar <- function(x,y) if (y > 0) x else y*x  # if (y > 0) 構文はベクトル化されていない
> bar(1,(-1):5)                               # つまり if((-1) > 0) 1 else ((-1):5)*1        
[1] -1  0  1  2  3  4  5
Warning message:
In if (y > 0) x else y * x :
   条件が長さが2以上なので,最初の一つだけが使われます 
> Bar <- Vectorize(bar, "y")                  # y変数をベクトル化する
> Bar(1,(-1):5)
[1] -1  0  1  1  1  1  1

(150) try() 関数の silent 引数(初級Q&A記事より、2008.02.11)

エラーが起きても実行を中断しない try() 関数はシミュレーション等で便利であるが、エラーメッセージそのものはコンソール表示される。これを抑制するオプション引数が silent (既定は FALSE でエラー出力)。但し抑制することは一般には問題を引き起こす可能性がある。以下はある操作がエラーを起こせば FALSE さもなければ TRUE を返すコード。同様の関数に tryCatch() がある。

> x1 <- matrix(1:6,ncol=3) 
> x2 <- matrix(1:4,ncol=2)
> colnames(x1) <- c("A","B","C")
> colnames(x2) <- c("A","C") 
> ifelse(class(try(y <- x2[,"B"])) == "try-error", FALSE, TRUE)
Error in try(y <- x2[, "B"]) :  添え字が許される範囲外です 
[1] FALSE
> y
 エラー:  オブジェクト "y" は存在しません 
> ifelse(class(try(y <- x1[,"B"])) == "try-error", FALSE, TRUE)
[1] TRUE
> y
[1] 3 4
## silent オプション引数でエラーメッセージを抑制する
> ifelse(class(try(y <- x2[,"B"], silent=TRUE)) == "try-error", FALSE, TRUE)
[1] FALSE

(149) リスト x の長さを length(x) <- n で変更すると成分 NULL で拡大される。

学生に教えてもらった例(とうに既知?)。

> ( x <- as.list(c(1,2)) )
[[1]]
[1] 1

[[2]]
[1] 2
> length(x) <- 4
> x
[[1]]
[1] 1

[[2]]
[1] 2

[[3]]
NULL

[[4]]
NULL

> x <- data.frame(runif(3), 1:3) # データフレームの場合
> length(x) <- 4 
> x
$runif.3.
[1] 0.3874197 0.8127907 0.8422990

$X1.3
[1] 1 2 3

[[3]]
NULL

[[4]]
NULL

help(length) [#z8986722]

Description: Get or set the length of vectors (including lists) and factors, and
of any other R object for which a method has been defined.~
Details: The replacement form can be used to reset the length of a vector. If 
a vector is shortened, extra values are discarded and when a vector is lengthened, 
it is padded out to its new length with NAs (nul for raw vectors).
> x <- 1:10
> length(x) <- 5 # 切り詰めることもできる
> x
[1] 1 2 3 4 5

(148) 大量で時間のかかるオブジェクトを作る際の工夫(2007.10.27)

各々がそれなりに時間のかかる大きなオブジェクトを大量に作りたいとする.作業途中にパソコンを他の用に使いたいこともあるし、ブレーカーが飛ぶこともあるし(我が家では珍しくない)、複数のパソコンで分散処理したいこともあるし、何よりも全部のオブジェクトを少ないメモリに同時に抱えておくこともできない,とする。困った挙げ句に思い付いた一つの工夫.要するに一つのオブジェクトが計算できたら、連番ファイルに順にセーブしておき、本番では逆にそれらを順に読み込んで使う.

## 途中で中断しても、前回計算が終わった次から続行
## for(i in 3000:4000) 等として別々のパソコンで分散計算も可能
for(i in 1:1e4) {
  filename <- paste("result.", i, ".data", sep="") # 連番ファイル名文字列作成
  if (file.exists(filename)) next               # 該当ファイルが既にあれば、パス
  ......                                           # ある作業、結果はオブジェクト x
  save(x, file=filename)                           # i番目のオブジェクトをセーブ
}
## 全てのファイル"result.1.data",...,"result.10000.data"が準備できたとする
## 今度は、それらを順に読み込んで次の作業を行う  
for (i in 1:1e4) {
  filename <- paste("result.", i, ".data", sep="") # 連番ファイル名文字列作成
  load(file=filename)                     # i番目のxを読み込む 
  ......                                  # i番目のxを使った作業  
} 

(147) データフレームの列の指定 Part2(2007.09.17)

(146)と関連するが,以下の違いを覚えておくと吉。

> x <- data.frame(a=1:3, b=rnorm(3), c=runif(3), d=3:1)
> x
  a           b         c d
1 1 -0.74118966 0.6880276 3
2 2  0.01489056 0.3304210 2
3 3  0.80607979 0.8301477 1
> x[,2]
[1] -0.74118966  0.01489056  0.80607979
> class(x[,2])
[1] "numeric"
> x[2]
            b
1 -0.74118966
2  0.01489056
3  0.80607979
> class(x[2])
[1] "data.frame"
> x[,c(2,4)]
            b d
1 -0.74118966 3
2  0.01489056 2
3  0.80607979 1
> class(x[,c(2,4)])
[1] "data.frame"
> x[c(2,4)]
            b d
1 -0.74118966 3
2  0.01489056 2
3  0.80607979 1
> class(x[c(2,4)])
[1] "data.frame"

(146) データフレームの列の指定(2007.09.14)

久保さんのブログ9月12日で,

> # たとえばこういう data.frame があったとする
> d
     x    y
1  1.6  1.1
2  1.1  0.6
3  0.4 -0.7
4 -0.3  0.8
5  0.2 -0.1
> # apply() と sprintf() を組みあわせて新しい label 列が作れる
> apply(d, 1, function(r) sprintf("xy=%g:%g", r["x"], r["y"]))
[1] "xy=1.6:1.1"  "xy=1.1:0.6"  "xy=0.4:-0.7" "xy=-0.3:0.8" "xy=0.2:-0.1"
> # あるいは paste() を使ってもよい
> d$label <- apply(d, 1, function(r) paste("xy", r["x"], r["y"]))

というのがあり,applyを使わない方法でうっかり次のようにやってしまった。

> paste("xy=", d["x"], ":", d["y"], sep="")
[1] "xy=c(1.6, 1.1, 0.4, -0.3, 0.2):c(1.1, 0.6, -0.7, 0.8, -0.1)"

正しくは,こっちの方だ。カンマがあるとないじゃ大違い。(それぞれの指定法にはちゃんと意味があり,場合によってどちらを使うかが決まるわけだが)

> paste("xy=", d[,"x"], ":", d[,"y"], sep="")
[1] "xy=1.6:1.1"  "xy=1.1:0.6"  "xy=0.4:-0.7" "xy=-0.3:0.8" "xy=0.2:-0.1"

割り込みコメント:次の方法がより直感的でわかりやすい? (それにしても、なぜこんな欄を作る必要があるのか、そこがわからぬ.)

> ( transform(d, label=paste("xy=",x,":",y,sep="")) )
     x    y          label
1  1.6  1.1  xy=1.6:1.1
2  1.1  0.6  xy=1.1:0.6
3  0.4 -0.7 xy=0.4:-0.7
4 -0.3  0.8 xy=-0.3:0.8
5  0.2 -0.1 xy=0.2:-0.1

書式からわかるように transform 関数を使った例は d を展開した環境中で評価(?)し、変数を一々取り出さない.他の2例は隠されたループ毎にデータフレーム(リスト)から変数を一々取り出している(のでどうしても時間を食いがちになる)。もちろんこの効果が目にみえるのは相当大きなデータフレームの場合.

> invisible({rm(d);gc();gc()}); set.seed(1); d <- data.frame(x=1:1e5, y=rnorm(1e5))
> system.time(d <- transform(d, lable=paste("xy=",x,":",y,sep="")))
   ユーザ  システム      経過 
    2.325     0.028     2.431 
> invisible({rm(d);gc();gc()}); set.seed(1); d <- data.frame(x=1:1e5, y=rnorm(1e5))
> system.time(d$label <- apply(d, 1, function(r) paste("xy", r["x"], ":", r["y"],sep="")))
   ユーザ  システム      経過 
   10.020     0.244    10.515 
> invisible({rm(d);gc();gc()}); set.seed(1); d <- data.frame(x=1:1e5, y=rnorm(1e5))
> system.time(d$label <- apply(d, 1, function(r) sprintf("xy=%g:%g", r["x"], r["y"])))
   ユーザ  システム      経過 
    4.669     0.032     4.814 

以下でも良いが何故かやたらに時間がかかるのが奇妙。sprintf は短距離型? 名前からして sprinter に似ているし(笑)。

> invisible({rm(d);gc();gc()})
> set.seed(1); d <- data.frame(x=1:1e5, y=rnorm(1e5))
> system.time(d <- transform(d, label=sprintf("xy=%g:%g",x,y)))
   ユーザ  システム      経過 
  389.320     1.744   409.348 

割り込みコメント2:次の方法がより直感的でわかりやすいく早いと思われます

d$label <- paste("xy=", d$x, ":", d$y, sep="")
d$label <- sprintf("xy=%g:%g", d$x, d$y)

applyを使うのはベクトルで操作出来ない場合で, ベクトル操作が可能なものは そのまま使えば良いです. 習得につれてapply中毒症状を起こす方が多いと思いますが, 原点にたち帰り,

d$x * d$y

が何を返すか, 関数型?Rにおけるベクトルのなんたるかを思い返せば良いでしょう applyと書いたら, applyを外して出来ないか考え直すのが良いでしょう.

(145) オブジェクトのコンソールへの整形出力の微妙な違い (2007.09.10)

以下の例の微妙な差に注意。
知っている必要はほとんど無いが、もし将来R本を書こうとする人は知っておく必要がある(出版社の校正係を悩ませるから :-)。

> (x <- 1)
[1] 1
> (x <- 1:2)
[1] 1 2
> (x <- 1:3)
[1] 1 2 3
## 途中省略
> (x <- 1:9)
[1] 1 2 3 4 5 6 7 8 9
> (x <- 1:10) # ここから [1] の先頭に空白が入るようになる
 [1]  1  2  3  4  5  6  7  8  9 10
> (x <- 1:11)
 [1]  1  2  3  4  5  6  7  8  9 10 11
> (x <- 1:20)
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
> (x <- 1:30)  # これは当然
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30
  • そして,当然こうなる

    (x <- 1:1000)

       [1]    1    2    3    4    5    6    7    8    9   10   11   12
      [13]   13   14   15   16   17   18   19   20   21   22   23   24
       :
      [97]   97   98   99  100  101  102  103  104  105  106  107  108
     [109]  109  110  111  112  113  114  115  116  117  118  119  120
       :

(144) quote() 関数と eval() 関数 (2007.09.04)

中級Q&A記事に関連して power.prop.test のコードを眺めていて、感心した例。

> body <- quote(f(x,y,z)) # f(x,y,z) はある3変数関数
> foo <- function(x=NULL, y=NULL, z=NULL) {
           ## 引数 x,y,z のうちどれかふたつは実引数を与える
           if (sum(c(is.null(x),is.null(y),is.null(z))) != 1) 
             stop("just one of x,y,z should be NULL")
           ## 関数 x → foo(x,y,z) を用いた作業(y,z は given)
           if(is.null(x)) {bar <- function(x) eval(body);
                    .............................
                          } 
           ## 関数 y → foo(x,y,z) を用いた作業(x,z は given)
           else if(is.null(y)) {bar <- function(y) eval(body);
                                ............................ 
                               }
           ## 関数 z → foo(x,y,z) を用いた作業(x,y は given)
           else {bar <- function(z) eval(body);
                 .................................
                }
         }
> foo(,2,3) # 使用例
> foo(1,,3)
> foo(1,2,)

(143) ? 関数の隠し機能(?) (2007.09.02)

知らなくても困らない機能頁に help() 関数の省略形 ? の奇妙な動作例が紹介されていた。これは意外に知っていると役に立つかもしれないので転載。

> matrix(1:16, 4,3,2)    # うっかり変なコードを入力し注意をくらう
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12
Warning message:
 行列のデータ長 [16] が列数  [3] を整数で割った、もしくは掛けた値ではありません   in: matrix(1:16, 4, 3, 2) 
> ?matrix(1:16, 4,3,2) # ヒストリ機能を使いすかさず matrix 関数のヘルプを参照

しかし、いつでもうまくいくとは限らないのが不思議.以下の例でわかるように、? 関数はこの場合 .helpForCall? という呼出し式に関連するヘルプドキュメントを表示する関数を起動しているらしい。

> ?log(1)        
> ?gamma(1)
> ?sin(1)               # 三角関数族は何故か不調 
 以下にエラー.helpForCall(e1Expr, parent.frame()) : 
        no documentation for function 'sin' and signature 'x = “numeric”'
 追加情報: Warning message:
no method defined for function 'sin' and signature 'x = “numeric”' in:  .helpForCall(e1Expr, parent.frame()) 

もう一つ奇妙なこと。この場合 ? は help() の簡略形では無い.

> help(log(1))
No documentation for 'log(1)' in specified packages and libraries:
you could try 'help.search("log(1)")'

注:実は隠し機能ではなく,help(help) に説明があった。これによれば method?lm といった形式のもうひとつのヘルプ参照法もある.

(142) 引数として渡されたオブジェクト名を得る(2007.08.19)

deparse(substitute(foo)) とする。これを使うとデバグ時の変数の表示がわかりやすくなるかも。

> dmp <- function(x) cat(deparse(substitute(x)),"=",x,"\n")
> x <- 3
> dmp(x)
x = 3 
> y <- "abc"
> dmp(y)
y = abc 

(141) 数値 1,0 は論理値が要求される局面では TRUE,FALSE とされる(2007.08.18)

何を今さらと思われるかもしれないが、案外以下のような場合律義に replace=TRUE としていませんか(私のことです)。大文字の入力は面倒だし,FALSEはよく指がもつれるし。

> sample(1:10, 10, replace=1) # sample(1:10,10,r=1), sample(1:10,10,1) 等とすればもっと簡単
 [1] 1 8 6 1 9 8 5 1 6 7
> sample(1:10, 10, replace=0)
 [1]  5  4  7  6  3  8  1  2 10  9
  • 数値は全て論理値として解釈される。0 は FALSE,「0 以外」はTRUE。以下のような場合はトリビアルであるが,プログラムを書くときに便利に使える局面がある。ただ,プログラムの可読性は低くなる。(「関数達」から外れてきた)
    > sample(1:10, 10, replace=3.14)
     [1] 1 9 6 5 3 5 6 2 1 7
  • sample(1:10, ...) は sample(10, ...) と同じである。
  • 逆に,数値演算をするときなどは TRUE と FALSE はそれぞれ,1 と 0 とされる。
    > sum(c(TRUE, FALSE, TRUE))
    [1] 2
    > TRUE+1
    [1] 2
    > TRUE*10
    [1] 10

(140) 行列の一重括弧演算子には添字行列が使える (2007.08.04)

但し結果はベクトルになる。

> x
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1   11   21   31   41   51   61   71   81    91
 [2,]    2   12   22   32   42   52   62   72   82    92
 [3,]    3   13   23   33   43   53   63   73   83    93
 [4,]    4   14   24   34   44   54   64   74   84    94
 [5,]    5   15   25   35   45   55   65   75   85    95
 [6,]    6   16   26   36   46   56   66   76   86    96
 [7,]    7   17   27   37   47   57   67   77   87    97
 [8,]    8   18   28   38   48   58   68   78   88    98
 [9,]    9   19   29   39   49   59   69   79   89    99
[10,]   10   20   30   40   50   60   70   80   90   100
> i
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
[3,]    3    5    7
> x[i]          # つまり as.vector(x)[as.vector(i)] と同じ
[1] 1 2 3 3 4 5 5 6 7
  • 添え字として使われる行列が2列だと,もっと役に立つよ。
    三次元配列に対しては3列を持つ行列...,以下同様
    > x <- matrix(1:100, 10)
    > i <- matrix(c(1,1, 2,2, 3,3, 9,2), byrow=TRUE, ncol=2)
    > i
         [,1] [,2]
    [1,]    1    1
    [2,]    2    2
    [3,]    3    3
    [4,]    9    2
    > x[i] # x[1,1], x[2,2], x[3,3], x[9,2] を意味することになる
    [1]  1 12 23 19

(139) 無駄でない知識、論理和と論理積の優先度 (2007.08.02)

多くの計算機言語(本来のS言語を含む)と異なり、Rでは論理積は論理和よりも演算子の優先度が高い。

> c( T|F&F, (T|F)&F, T|(F&F) ) 
[1]  TRUE FALSE  TRUE          # T|F&F は(T|F)&F でなく T|(F&F)  
  • 殆どの言語はRと同じく論理積の方が高いと思いますが・・・
  • 紛らわしい(自信がない)場合は,積極的に括弧を使いましょう・・・(プログラムが読みやすくなりますし)
  • help("|") より抜粋。
    See Syntax for the precedence of these operators: unlike many
        other languages (including S) the AND and OR operators do not have
        the same precedence (the AND operators are higher than the OR
        operators).
  • S では確かに論理演算子(| & || &&)は同順位になっている。しかし,ほかの多くのプログラミング言語では &,&& は |,|| より優先されるので,help の記述が不適切ということでは?

(138) 無駄な知識、関数と引数リストの間には任意個数の半角空白をおける (2007,07.31)

if 文の書式は if(cond) でも if (cond) でも良いことは良く知られているが、もしかしてと思い確かめてみたら、これは一般関数でも同じだった。

> log (1)
[1] 0
> log  (1)
[1] 0
> order (3:1)
[1] 3 2 1
> order                  (3:1)
[1] 3 2 1
  • 配列の [ ] も同じですよ
    > x <- sqrt(1:10)
    > x    [6]
    [1] 2.449490
    > "["    (x, 6) # ついでに
    [1] 2.449490
    > "order"      (3:1) # もひとつ,ついでに
    [1] 3 2 1
    > ":"         (1,3) # 最後に,もひとつ
    [1] 1 2 3
    > ifelse  (1:3%%2, letters[1:3], LETTERS[1:3]) # ifelse は関数だから
    [1] "a" "B" "c"
    > "ifelse"  (1:3%%2, letters[1:3], LETTERS[1:3]) # これはエラーではないが
    [1] "a" "B" "c"
    > if    (3 > 2) print(TRUE) # if は「関数ではない」から
    [1] TRUE
    > "if"    (3 > 2) print(TRUE) # こんなことするとエラーになる
     エラー:syntax error
  • if 文は関数です。
    > "if"(3 > 2, print(TRUE))
    [1] TRUE
    > "if" (3 > 2, print(TRUE))
    [1] TRUE
  • そうでした(_ _)

(137) Rはコードの意図を見抜く?(2007.07.24)

次の行列 x を行列 x0 で初期化する二つのコードの計算時間比はどう解釈したらよいのでしょう。まるで R は二つ目のコードは無駄な繰り返しと見抜いているような(!?)。

> x <- x0 <- matrix(0, 1e3, 1e3)
> system.time(for (n in 1:1e4) {x[1,1] <- 1; x <- x0})
   ユーザ  システム      経過 
   49.787   165.662   217.995 
> system.time(for (n in 1:1e4) x <- x0)
   ユーザ  システム      経過 
    0.004     0.000     0.003 
  • プロファイル見るとgc(ガベージコレクション)が走ってます.(ループ回数) この条件を内部で最適化出来ると, 嬉しいように思います...(けど,難しい...)
    # gcの対象を小さくするのにlocalを使う実験
    > x <- x0 <- matrix(0, 1e3, 1e3)
    > system.time(for (n in 1:1e3) {x[1,1] <- 1; x <- x0})
       ユーザ  システム      経過
       13.213     6.705    19.926
    > system.time(for (n in 1:1e3) local({x[1,1] <- 1; x <- x0}))
       ユーザ  システム      経過
        7.861     6.592    14.453
    # "[<-"(subassign)がgcを誘発している
  • ということは最初の二番目のコードの実行時間は文字通りにうけとってよいということでしょうか。いくらなんでも早すぎるような気がしたので不審におもったのですが。
    > system.time(for (n in 1:1e4) x <- 0)
       ユーザ  システム      経過 
        0.004     0.000     0.040 
  • 今どきのコンピュータならそんなものです. 演算自体は水面の白鳥でgcはさしずめ水掻きのようなものです. 通常はgcが余計な物を消して性能に貢献してるんですが, 反復が多すぎると欠点として表れる事もあります.
    # localを使ってgcの影響を排除すると, localの影響になる
    > system.time(for (n in 1:2e2) local({x[1,1] <- 1; x <- x0}))
       ユーザ  システム      経過
        2.676     1.389     4.114
    > system.time(for (n in 1:2e2) local({x[1,1] <- 1; y <- x0}))
       ユーザ  システム      経過
        2.816     1.256     4.090
    # ,subAssignはgcを誘発する. Rはxの値を捨てていいのか判らないが, リプレスされれば気づく.
    > system.time(for (n in 1:2e2) {x[1,1] <- 1; x <- x0})
       ユーザ  システム      経過
        5.668     1.376     7.086
    > system.time(for (n in 1:2e2) {x[1,1] <- 1; y <- x0})
       ユーザ  システム      経過
        0.004     0.012     0.017
    通常はもっとコストの高い演算を行うケースの方が多いので, gcによる負荷は考えなくても良いはずです.

(136) apply の第二引数はベクトルも取れる(2007.07.12)

> x <- 1:24
> dim(x) <- c(3,4,2)
> x
, , 1

     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

, , 2

     [,1] [,2] [,3] [,4]
[1,]   13   16   19   22
[2,]   14   17   20   23
[3,]   15   18   21   24

> apply(x, 1, sum)
[1]  92 100 108
> apply(x, 2, sum)
[1]  48  66  84 102
> apply(x, 3, sum)
[1]  78 222
> apply(x, 1:2, sum)
     [,1] [,2] [,3] [,4]
[1,]   14   20   26   32
[2,]   16   22   28   34
[3,]   18   24   30   36
> apply(x, c(1,3), sum)
     [,1] [,2]
[1,]   22   70
[2,]   26   74
[3,]   30   78
> apply(x, 2:3, sum)
     [,1] [,2]
[1,]    6   42
[2,]   15   51
[3,]   24   60
[4,]   33   69
> apply(x, 1:3, sum)
, , 1

     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

, , 2

     [,1] [,2] [,3] [,4]
[1,]   13   16   19   22
[2,]   14   17   20   23
[3,]   15   18   21   24

(135) データフレームの変数の成分を取り出す処理速度の比較(2007.07.01)

同値な4種類の操作の速度比較。結構違う!変数名でアクセスするのが一番早いらしい。

> z <- data.frame(x=runif(1e3), y=sample(letters[1:10], 1e3, replace=TRUE))
> system.time(for (n in 1:1e4) z$x[1:100])
   ユーザ  システム      経過 
    0.088     0.004     0.091 
> system.time(for (n in 1:1e4) z[[1]][1:100])
   ユーザ  システム      経過 
    0.164     0.004     0.171 
> system.time(for (n in 1:1e4) z[[c(1,1)]][1:100])
   ユーザ  システム      経過 
    0.196     0.004     0.197 
> system.time(for (n in 1:1e4) z[1][[1]][1:100])
   ユーザ  システム      経過 
    1.664     0.008     1.705 

因子変数は少し様子が違う?

> system.time(for (n in 1:1e4) z[[c(2,1)]][1:100])
   ユーザ  システム      経過 
    0.188     0.000     0.190 
> system.time(for (n in 1:1e4) z$y[1:100])
   ユーザ  システム      経過 
    0.464     0.000     0.469 
> system.time(for (n in 1:1e4) z[[2]][1:100])
   ユーザ  システム      経過 
    0.592     0.000     0.590 
> system.time(for (n in 1:1e4) z[2][[1]][1:100])
   ユーザ  システム      経過 
    2.024     0.000     2.053 

(134) 因子の操作は時間がかかる (r-help 記事より。2007.06.24)

リスト(データフレームを含む)の操作はベクトル(行列、配列を含む)より格段に時間がかかることは R FAQ であるが、もう一つ要注意なのは、因子の操作は元のベクトルの操作に比べ時間がかかること。これは、因子 y は、水準ベクトル levels(y) で間接的に表現されていることを理解すれば当然。データフレーム中の変数はしばしば(知らぬ間に)因子化されてしまう。さらにデータフレームはリストだから、余計に時間を食う。数値データフレームの処理に時間がかかりすぎる際は、数値行列に変換(t(t(DF)) で変換できる)して操作すると相当早くなることがある。

> x <- sample(0:100, 1e3, replace=TRUE)
> y <- as.factor(x)
> str(y)
 Factor w/ 101 levels "0","1","2","3",..: 100 58 60 14 29 27 74 40 88 34 ...
> levels(y)  # 因子の水準ベクトル(文字列)
  [1] "0"   "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "9"   "10"  "11" 
 [13] "12"  "13"  "14"  "15"  "16"  "17"  "18"  "19"  "20"  "21"  "22"  "23" 
 (途中省略)
 [85] "84"  "85"  "86"  "87"  "88"  "89"  "90"  "91"  "92"  "93"  "94"  "95" 
 [97] "96"  "97"  "98"  "99"  "100"
> y[1]  # 表示の際は整形されるので、水準ベクトルがくっつく他は x[1] (=2) と同じに見えるが
 [1] 2
 101 Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ... 100
> str(y[1]) # 実際は 3 (水準ベクトルの3番目の意味)と記録されている
 Factor w/ 101 levels "0","1","2","3",..: 3
> levels(y)[y[1]]  # つまり x[1] (の文字版)
[1] "2"
> as.numeric(levels(y)[y[1]]) # これでやっと x[1] が復元できる
[1] 2
> all.equal(x, as.numeric(levels(y)[y]))
 [1] TRUE
## x と y の操作時間の比較(約70倍)
> system.time(for (n in 1e4) for(i in 1:1e3) x[i])
   ユーザ  システム      経過
    0.000     0.000     0.001
> system.time(for (n in 1e4) for(i in 1:1e3) y[i])
   ユーザ  システム      経過
    0.064     0.004     0.067 

===== 先日初級Q&に投稿された件について

> df <- data.frame(table(sample(1:10*5, 2000, replace=TRUE)))
> system.time(for (i in 1:1e4) as.numeric(levels(df[,1]))[df[,1]])
   ユーザ  システム      経過 
    1.922     0.037     2.311 
> system.time(for (i in 1:1e4) as.numeric(as.character(df[,1])))
   ユーザ  システム      経過 
    1.466     0.027     1.761 
推奨されるやり方は時間が掛かる
> f <- df[,1]
> system.time(for (i in 1:1e4) as.numeric(levels(f))[f])
   ユーザ  システム      経過 
    0.569     0.022     0.705 
> system.time(for (i in 1:1e4) as.numeric(as.character(f)))
   ユーザ  システム      経過 
    0.730     0.013     0.863 
いったんベクトルに代入してから変換すると,推奨されるやり方が早い
しかしいまあ,それは相対的な話に過ぎなくて,絶対的な速度差は無視できる

=== 上のコメント例はむしろデータフレーム操作は時間を食う、という通則に該当しそうです。データフレームを二度参照するのと一度参照するの差では。
=== 注意:因子は値でグループ化されている。したがって sort(), order(), unique() 等グループとして操作可能な処理はかえって早くなる可能性がある。

(133) is.integer 関数(2007.05.31)

整数を表すときには,L を添えよう。

> is.integer(9.8)
[1] FALSE
> is.integer(9)
[1] FALSE
> is.double(9)
[1] TRUE
> is.integer(9L)
[1] TRUE

(132) テキストコネクション (2007.05.30)

テキストコネクションは文字列(ベクトル)をあたかもテキストファイルであるかのように扱うことを可能にする。これと、ファイル読み込み関数を組み合わせると結構便利に使える。

(1) 数値列を表す文字列を、数値ベクトルに変換(専用関数 type.convert() もあるが)。

> x <- "123, 12.3, 1.23, 1.23e-1, 1.23e1"    # 数値列を表す文字列
> scan(textConnection(x), sep=",")   
[1] 123.000  12.300   1.230   0.123  12.300

(2) スラッシュ記号で区切られた文字列ベクトルをデータフレームにする。

> t <- c("a/a/a","b/bb/bbb","ccc/cc/c")
> read.table(textConnection(t), sep = "/")
   V1 V2  V3
1   a  a   a
2   b bb bbb
3 ccc cc   c
> t <- c("a/a/a","b/bb/bbb","ccc/cc") # 不揃いな場合
> read.table(textConnection(t), sep = "/", fill=TRUE)
   V1 V2  V3
1   a  a   a
2   b bb bbb
3 ccc cc    

(3) 文字列を単語に分解(空白で分離)する。

> x <- "To R or not to R, that is a question."
> scan(textConnection(x), what="")
Read 10 items
 [1] "To"        "R"        "or"        "not"       "to"        "R,"      
 [7] "that"      "is"        "a"         "question."

(4) 書き込みの例(?textConnection より)

> zz <- textConnection("foo", "w") # テキストコネクションを書き込みモードで開く
> ctl <- c(4.17, 5.58, 5.18, 6.11, 4.5, 4.61, 5.17, 4.53, 5.33, 5.14)
> trt <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69)
> group <- gl(2, 10, 20, labels = c("Ctl", "Trt"))
> weight <- c(ctl, trt)
> sink(zz)  # テキストコネクションへの sink 
> anova(lm.D9 <- lm(weight ~ group))
> cat("\nSummary of Residuals:\n\n")
> summary(resid(lm.D9))
> sink()   # sink 終了
> close(zz)  # テキストコネクションを閉じる
> foo    # 文字列ベクトル foo が作られている  
 [1] "Analysis of Variance Table"                                        
 [2] ""                                                                  
 [3] "Response: weight"                                                  
 [4] "          Df Sum Sq Mean Sq F value Pr(>F)"                        
 [5] "group      1 0.6882  0.6882  1.4191  0.249"                        
 [6] "Residuals 18 8.7293  0.4850               "                        
 [7] ""                                                                  
 [8] "Summary of Residuals:"                                             
 [9] ""                                                                  
[10] "      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. "
[11] "-1.071e+00 -4.938e-01  6.850e-02  7.109e-18  2.462e-01  1.369e+00 "

(131) 時々勉強し直さないと (2007.05.17)

単に気づかなかっただけか、それとも最近になって付け加えられたのか、見知らぬ機能を幾つも発見(オールドユーザー)。
関数 seq() の仲間:(2.4.0からだそうです)

> seq_len(5)               # seq(length.out=) の整数高速版
[1] 1 2 3 4 5
> seq_along(letters[1:5])  # seq(along.with=) の整数高速版
[1] 1 2 3 4 5
> seq.int(3, 10.5, by=2)   # seq() の高速版(int は内部関数 internal の意味)
[1] 3 5 7 9

関数 system.time() は既定で時間計測前にガベージコレクションを行う。 (TRUEになったのは2.2.0から。2.3.0からは測定精度がmsになった。)

system.time(expr, gcFirst = TRUE)

(130) ifelse 文も実際は関数だから (from r-help, 2007.04.18)

ifelse 文は実際は関数 "ifelse"(cond, expr1, expr2) で、選択された式を返り値にもつ。そんなわけで次のような一見奇妙な現象が起きる。

> ifelse(TRUE, print("yes"), print("no"))
[1] "yes"  選択された式 print("yes") の実行結果
[1] "yes"  返り値 print("yes") が実行された結果
> ifelse(TRUE, "yes", "no")
[1] "yes"  この際は返り値は文字列 "yes" 
> x <- ifelse(TRUE, print("yes"), print("no"))
[1] "yes"
> x
[1] "yes"

(129) ComputerModern?フォントでEPSファイルを作る (2007.04.17)

EPS形式の出力にもあるように、EPSファイルを作る際に"family=ComputerModern?"を指定するとCMフォントで代替してくれるので、TeXに取り込む際に少し嬉しい。ただし以下の点には注意

  • 余程凝った数式などを入れたい場合には、別途PSfragを使うなりする方が良い
  • 学会によってはAdobe純正のPSファイルじゃないと受け付けてくれない場合がある(電子情報通信学会とか)
  • 作られたEPSファイルには、フォントは埋め込まれておらず(指定されるだけ)、dvips/dvipdfmxを経て初めて実フォントが埋め込まれる
  • 埋め込み前の段階でGhostscript/GhostView?で見たいのなら、GS側にもCMフォントを登録する必要あり

(128) クイズ (2007.03.29)

公式マニュアル「The R Language Definition」中の例を確認中に出会ったRの一見不可解な挙動(マニュアルには foo(x <- y) という構文は使うなと注意があります)。

> foo <- function(x) x
> foo1 <- function(x) x^2
> foo2 <- function(x) (x)
> foo3 <- function(x) 4*x
> rm(x); y <- 3
> foo(x <- y)  # 何故か反応無し
> foo1(x <- y)
[1] 9
> foo2(x <- y)
[1] 3
> foo3(x <- y)
[1] 12
> foo(3)
[1] 3
> foo(y)
[1] 3
> foo((x <- y)) # こうすると答えが得られる
[1] 3
> x; y
[1] 3
[1] 3

(127) 関数のコメント (2007.03.23, from r-help)

Rにおけるコメントについては、一行ならば #、複数行ならば if(0){...} とすることは良く知られているが、いまひとつの簡単な方法がある。これは関数中におかれた文字列は実質上無意味であることを利用している.

> foo <- function(x) {
+ "これは二行にわたる、長い
+ コメントです。"
+ x^2}
> foo(3)
[1] 9
> foo
function(x) {
"これは二行にわたる、長い
コメントです。"
x^2}

注:if (0) { } でコメントアウトする場合は, { } の部分も文法に従っていないとシンタックスエラーになることに留意のこと(つまり,なんでも書ける訳ではない。解釈されるけど実行されないだけである。当然ながら," " のわざも,途中に " が出てくるような注釈を書くときには,気をつけないといけない。
注:直前の注にある問題は、一重・二重引用符を併用すると防げる。

> foo <- function(x) {
 'cat("This is the function foo\n") # 二行コメント化
  cat("argument x is",x,"\n")'
 x^2}
> foo(3)
[1] 9

(126) ベクトル化定数値関数、初級Q&A記事より (2007.03.01)

ベクトル化された関数を引数に要求する関数(例えば一次元数値積分関数 integrate())に定数値関数をあたえるには?もちろん function(x) 1 等では駄目.

> X <- runif(10)
> f0 <- function(x) 1
> f1 <- function(x) 1 + 0*x # ベクトル化定数値関数
> f0(X)
[1] 1
> f1(X)
 [1] 1 1 1 1 1 1 1 1 1 1

max() 関数等ベクトル引数をとるが、値がスカラーになる関数には、親切にもベクトル化版が用意されている.

> g0 <- function(x) max(x,1)
> g1 <- function(x) pmax(x,1)
> g0(X)
[1] 1
> g1(X)
 [1] 1 1 1 1 1 1 1 1 1 1

(125) browser() 関数は実行を一時中断し、それを呼び出した環境をチェックすることを可能にする (2007.02.26)

関数呼び出しが複雑な入れ子になっている場合のデバッグに使う?

> f0 <- function() {
          f1 <- function() browser()
          f2 <- function() f1()
          f2()
        }
> f0()
Called from: f1()  # 環境ブラウザーがf1()から呼び出された
Browse[1]> where   # その時の環境スタックを表示
where 1: f1()
where 2: f2()      # f1()はf2()から呼び出された
where 3: f0()      # f2()はf0()から呼び出された
Browse[1]> Q       # ブラウズ中止

(124) 環境オブジェクトは参照渡しが原則 (2007.02.25)

環境(変数名とその値の対の集合であるフレームと、その親(囲み)環境へのポインタであるエンクロージャからなる)はRのオブジェクトであり、関数引数に使ったり、変数に付値できる。しかし、他のRオブジェクトと異なり、コピーされない(つまり参照渡しが原則)。

> e <- new.env()      # 新しい環境 e
> assign("x", 2, env=e)   # e のフレームに変数 x を値 2 で登録
> e1 <- e         # e を e1 に付値  
> get("x", env=e1)
[1] 2   
> assign("x", 3, env=e1) # e1 のフレーム中の変数 x の値を変える
> get("x", env=e)     # e 中の変数 x の値は?
[1] 3            # 変化している、つまり e と e1 は同一! 
> e2 <- .GlobalEnv        # 大局的環境を e2 に付値
> assign("z", letters[1:5], env=e)
> zエラー:オブジェクト "z" は存在しません
> assign("z", letters[1:5], env=e2)
> z
[1] "a" "b" "c" "d" "e"

(123) 新しい総称的関数 foo() とそのメソッド関数 foo.default(), foo.mean(), foo.median() を定義する例 (2007.02.25)

> foo <- function(x, ...) UseMethod("foo", x) # 総称的関数fooの定義 
> foo.default <- function(x) print(x)         # 既定メソッド関数
> foo.mean <- function(x) print(mean(x))      # クラス"mean"に対するメソッド関数
> foo.median <- function(x) print(median(x))  # クラス"median"に対するメソッド関数
> x <- runif(10); foo(x)                      # 明示的クラス属性がなければfoo.default()が使われる
[1] 0.6013047 0.1863832 0.3672611 0.2846233 0.6553732 0.4497602 0.7911005
[8] 0.8818249 0.8193772 0.2082053
> class(x) <- "mean"; foo(x)                  # foo.mean()が使われる
[1] 0.6644322
> class(x) <- "median"; foo(x)                # foo.median()が使われる
[1] 0.6781848

(122) Rの関数定義とスコープ規則に関する一見奇妙な例 (2007.02.22)

Rの公式マニュアル中の例(これも頭を悩ますバグの種になりそうな)

> foo <- function() {   # 関数を返す関数
+ y <- 10
+ bar <- function(x) x+y
+ return(bar)}
> h <- foo()            
> h(1:10)
 [1] 11 12 13 14 15 16 17 18 19 20
> y <- 3                # y の値を変えてみる        
> h(1:10)        # しかし結果は同じ
 [1] 11 12 13 14 15 16 17 18 19 20
> y <- 10
> hh <- function(x) x+y # 一見 h と同じ関数を定義
> hh(1:10)        # y は 10 だから h(1:10) と同じ   
 [1] 11 12 13 14 15 16 17 18 19 20
> y <- 3         # y の値を変えてみると?
> hh(1:10)        # 結果が変わる! 
 [1]  4  5  6  7  8  9 10 11 12 13
> h           # h と hh の定義の微妙な違いに注意  
function(x) x+y     # h の y は定義時に 10 に(関数定義環境中で)固定されている
<environment: 0x84b1be0>
> hh           # hh の y は特定の値に拘束されておらず、その値はスコープ規則で決められる
function(x) x+y

(121) 永続付値 <<- は必ずしも「永続」でない例 (2007.2.18)

永続付値演算子 x <<- value, value ->> x は,変数 x を検索パスに沿って探索し,それに付値する.もしそれが無ければ,大局的環境 .GlobalEnv? 中の変数 x (存在しなければ作成したうえで)に値を付値する。もし、大局的環境にいたる前に x という名前の変数が見つかればどうなる?

foo <- function(x) { x <- 0; bar <- function(y) x <<- y; bar(x); cat(x, fill=TRUE) }
> foo(3)
3                                        # これは foo() の実行環境中の x の値
> x 
エラー:オブジェクト "x" は存在しません  # foo() 実行後は x は存在しない!
  • 「永続」という呼び方に問題があるのではないでしょうか。大域変数への代入ですから,順次上を見ていって,存在したらそれに作用する。無ければトップレベルで作成の上,付値される。bar の上が foo で,そこにある x に付値する。そして,その x は foo の中の局所変数なので,関数から出たら無くなる。bar も同じく,なくなる。
    なお,foo(3) の結果は 0 になりました。
    > foo <- function(x) { x <- 0; bar <- function(y) x <<- y; bar(x); cat(x, fill=TRUE) }
    > foo(3)
    0
  • 失礼。以下のように訂正します。
    > foo <- function(x) { x <- 0; bar <- function(y) x <<- y; bar(3); cat(x, fill=TRUE) }
    > foo(1)
    3
    > x 
    エラー:オブジェクト "x" は存在しません
    もともとが permanent assignment ですから訳「永続」は仕方がないかと。大域変数でも誤解の程度は同じ位では(?)。恐らくこれは、本来のS言語では環境(S言語ではフレームというらしい)には関数の実行環境と、大局的環境の二種類しかない(らしい)ことから来た(その限りでは問題の無い)命名だったんだろうと考えます。

(120) リスト・データフレームに同一の変数名ラベルを与えると? (2007.2.11)

(119) の記事に関連して偶然(うっかりして)気付いた例。

> ( x <- list(a=1:2, b=3:4, c= runif(2), c=c("abc",'def')) )
$a
[1] 1 2
$b
[1] 3 4
$c
[1] 0.4171226 0.6537818
$c                       # 同一の変数名が共存
[1] "abc" "def"
> x[["c"]]               # 二重鈎括弧演算子では最初のものが取り出される 
[1] 0.4171226 0.6537818
> x[[4]]
[1] "abc" "def"

データフレームでは同じ変数名を与えると自動的に調整される

> x <- data.frame(a=1:2, b=3:4, c= runif(2), c=c("abc",'def'))
> x
  a b         c c.1
1 1 3 0.9148885 abc
2 2 4 0.1790133 def
> x[["c"]]
[1] 0.9148885 0.1790133
> x[["c.1"]]
[1] abc def
Levels: abc def

(119) 引用符あれこれ (2007.2.11)

Rではデータ、変数ラベル、関数名、関数仮引数名等として文字列(オブジェクト・ラベル名)が頻繁に登場する。文字列は二重引用符 " で囲むのが基本であるが、一重引用符 ' で囲んでも良い。オブジェクト表示では二重引用符で表現される。ベクトル、リスト、データフレームの変数ラベルは引用符で囲まなくても良いが、成分操作では引用符で囲む必要がある。特殊な引用符としてバックチック(backtick)記号 ` がある。これはオブジェクト名、関数仮引数ラベル、データフレームの変数名ラベルに使える(ほとんど冗長であるが)。バックチック記号が本質的に有用な例があったら教えてください。

>  c("a", 'b')     # 一重・二重引用符混用
[1] "a" "b"       # 表示では二重引用符になる
> list("a"=1:2, 'b'=c("abc",'def'))
$a
[1] 1 2
$b
[1] "abc" "def"
> "foo" <- function(x) x
> foo
function(x) x
> 'foo' <- function(x) x
> 'foo'(3)
[1] 3
> foo <- function(x=1, `y`=4) x*y
> foo()
[1] 4
> list(a=1:2, "b"=3:4, `c`= runif(2), 'd'=c("abc",'def'))
$a
[1] 1 2
$b
[1] 3 4
$c
[1] 0.8238691 0.5051710
$d
[1] "abc" "def"
> x[['a']]       # x[[`a`]] は不可
[1] 1 2
> x[["a"]]
[1] 1 2

バックチック記号の使用例。

> `x`                               # オブジェクト名
$a
[1] 1 2
$b
[1] 3 4
$c
[1] 0.8238691 0.5051710
$c
[1] "abc" "def"
> `foo` <- function(x=1, `y`=4) x*y # 関数名、仮引数ラベル
> foo()                             # "foo", `foo` でも可
[1] 4
> `foo`(`x`=2, y=3)
[1] 6
> foo <- function(x=1, "y"=4) x*y    # 仮引数名には一重・二重引用符は使えない 
エラー:"foo <- function(x=1, "y"" に構文エラーがありました

役に立つ(?)例が有りました。次の関数は value of y is という名前の仮引数をもちます(笑)。

> foo <- function(x=1, value of y is=3) x*y
エラー:"foo <- function(x=1, value of" に構文エラーがありました
> foo <- function(x = 1, `value of y is` = 3) x*{`value of y is`}
> foo()
[1] 3
> foo(3)
[1] 9
> foo(2, `value of y is` = 5)
[1] 10

要するにオブジェクト・ラベル名に句(複合文字列)を使うことを許す。実はRの言語そのものを扱う関数の内部で結構使われているようです。

> `My favourite thing` <- "R"
> `My favourite thing`
[1] "R"
> "A B C" <- 3  # エラーにはならない!
> "A B C"    # これでは文字通りの文字列     
[1] "A B C"
> `A B C`    # こうすれば大丈夫
[1] 3

(118) 関数の引数既定値の印象的な例 (2007.01.30)

Rの関数の引数既定値の与え方は融通無碍でしばしば驚かされるが、次も印象的な例。これは既定値というよりは既定動作。

> foo <- function(x, y=stop("you should provide y!"))  x*y
> foo(1,2)
[1] 2
> foo(1) 
以下にエラーfoo(1) : you should provide y!

こんなのもあり。

> foo <- function(x, y={bar <- function(z) log(z); bar(x)})  x*y
> foo(2,1)
[1] 2
> foo(2)
[1] 1.386294

(117) データセット用万能関数、その名も data() (2007.01.30)

データファイルの書式は色々あり、それぞれ読み込む関数も違う、ああ面倒、と思った人がいたのか、万能型データセット関数 data() がいつの間に(しばらく前にはなかったはず)か登場していた。

拡張子 .R か .r のファイル。source() 関数で読み込まれる
拡張子 .RData か .rda のファイル。load() 関数で読み込まれる
拡張子 .tab, .txt, .TXT のファイル。read.table(..., header = TRUE) を使い読み込まれる
拡張子 .csv か .CSV のファイル。read.table(..., header = TRUE, sep = ";") を使いデータフレームとして読み込まれる
data()                        # 利用可能なすべてのデータセットを表示
try(data(package = "rpart") ) # パッケージ rpart 中のデータセットを表示
data(USArrests, "VADeaths")   # データセット USArrests, VADeaths を読み込む

(116) 関数引数の遅延評価の教訓的例 (2007.01.28)

R関数の実引数の本当の値は「それが実際に必要になった時点で決まる」。だから、以下の例はまことにごもっともだが。

> foo <- function(x, y=x) {x <- x+1; y}
> foo(10, 1)
[1] 1
> foo(10) # y の値が必要になった時点では x の値は 11 だから y=11 とされる
[1] 11

(115) R 2.4.0 で登場の実験的機能 readNEWS (2007.01.28)

論ずるより打ち込むが易し、readNEWS() と してみてください。

(114) 行列に対する ifelse 関数 (2007.01.28)

r-help 記事を見ていて驚いた ifelse 関数の使い方。ifelse 関数はベクトル化されているの知っていたが、さらに行列(配列)化されている! つまり、次元属性を保存する。

> X
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> Y
     [,1] [,2] [,3]
[1,]   10   13   16
[2,]   11   14   17
[3,]   12   15   18
> Z
          [,1]      [,2]      [,3]
[1,] 0.4087897 0.4984087 0.1219828
[2,] 0.8372826 0.4903184 0.2826118
[3,] 0.4453304 0.8820485 0.3584058
> ifelse(Z <= 0.5, X, Y)
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]   11    5    8
[3,]    3   15    9

(113) テーブルのテーブル (2007.01.27)

データ中の項目の登場回数を集計するのが table(x) 関数だが、table(table(x)) とすると面白いことがわかる(学生に教えられ感心した例)。

> x <- sample(letters[1:10], 100, rep=TRUE)
> table(x)
x                 
 a  b  c  d  e  f  g  h  i  j  # "a","b","c",... がそれぞれ 14,11,8,... 回登場
14 11  8 13  9  7 10  9 10  9
> table(table(x))

 7  8  9 10 11 13 14           # 9回登場した項目が3個、等 
 1  1  3  2  1  1  1

(112) ターミナル(*term)のタイトルバーにカレントディレクトリを表示しておく

.Rprofileに以下を追加

if(interactive()){
  utils::assignInNamespace("setwd",function(x){
      .Internal(setwd(x))
      cat(paste("\x1b]0;!",
                R.Version()$version.string,
                "!",
                Sys.info()["login"],
                "@",
                Sys.info()["nodename"],
                ":",
                getwd(),
                "\x07",
                sep=""
                ))
  },"base")
  setwd(getwd())
}

(111) RGuiのタイトルバーにカレントディレクトリを表示しておく

.Rprofile.siteに以下を追加

utils::setWindowTitle(base::getwd())
utils::assignInNamespace("setwd",function(dir){
  .Internal(setwd(dir));setWindowTitle(base::getwd())},"base")

このwikiで昔似たようなものを見かけたような気もするんだけど・・・

Rguiを違うディレクトリでいくつも起動して作業する人には、便利だと思う。

(110) 関数を使いこなす ifelse

vector の中である条件をみたす値だけを別の vector の値で置き換えたい

> new <- c(10, NA, 21, 16, NA, 13)
> old <- c(9,   8, 20, 15, 12, 11)
> replacement <- is.na(new) # どこが欠側値か
> new[replacement] <- old[replacement] # 欠側値のみ置き換え (条件つきコピー!)
>  new # 欠側値を置き換えられた vector
[1] 10  8 21 16 12 13

以下のようにすると,わかりやすくもある

> new <- c(10, NA, 21, 16, NA, 13)
> old <- c(9,   8, 20, 15, 12, 11)
> new <- ifelse(is.na(new), old, new)
> new
[1] 10  8 21 16 12 13

(109) sapply の多変数バージョン: mapply

do.call と良く似たものかな?

> junk <- mapply(cat, list(1:3,4:6), fill=TRUE)
1 2 3
4 5 6

何かに付値しないと(例では junkに付値),ゴミ(不要なもの)が出力される。適用する関数によっては,do.call と異なった出力になることもある。 (89) も参照。

> a <- as.list(NULL)
> a <- c(list(1:3), list(4:9))
> a
[[1]]
[1] 1 2 3

[[2]]
[1] 4 5 6 7 8 9

> do.call(function(x) sum(1/x), a)
以下にエラーfunction (x)  : 使われていない引数 (4:9)
> mapply(function(x) sum(1/x), a)
[1] 1.833333 0.995635

(108) 関数を引数に適用する関数:do.call

(Schemer/Lisperのために)
Rのapply系関数はScheme/Lispのmap*関数です。
Scheme/LispのapplyはRではdo.callです。

一例

> cat(list(1:3,4:6))
以下にエラーcat(list(...), file, sep, fill, labels, append) : 
        引数 1 (タイプ 'list') は 'cat' で取り扱えません
>; do.call(cat,list(1:3,4:6))
1 2 3 4 5 6

最後の呼び出しは

> cat(1:3,4:6)

と等価です。

結果がlistで返される関数の戻り値をいろいろやりたいとき、lapply/sapplyを使う前に、do.callが使えないかを考えてみるといいかも。

===== 以下と,どう違うか?

> cat(unlist(list(1:3,4:6)))
1 2 3 4 5 6

===== 全然違いますよ(catの出力は一緒ですが)。

> f<-function(x1,x2=1)x1+x2
> f(unlist(list(1:3,4:6)))
[1] 2 3 4 5 6 7
> do.call("f",list(1:3,4:6))
[1] 5 7 9

ということです。

===== いや,同じだといっているのではなく,また,処理・関数によっては do.call と そのほかの候補はいつも同じではないですねといっているだけで,「do.call も検討したら?」というのに逆らっているわけではありません。

(107) 英単語の先頭だけ大文字にする

> ToUpper <- function(s) paste(toupper(unlist(strsplit(s,NULL))[1]), 
                                 paste(unlist(strsplit(s,NULL))[-1], collapse=""), 
                               sep="")
> ToUpper("abcdef")
[1] "Abcdef"
> ToUpper2 <- function(s) {substr(s, 1, 1) <- toupper(unlist(strsplit(s,NULL))[1]); return(s)}
> ToUpper2("abcdef")
[1] "Abcdef"
> ToUpper3 <- function(s) {substr(s, 1, 1) <- toupper(substr(s, 1, 1)); return(s)}
> ToUpper3("abcdef")
[1] "Abcdef"

(106) データフレーム操作に関するある問題(北大久保さんのブログより, 2007.01.22)

以下のようなデータフレーム data の group 変数の種類ごとに、y 変数が最大である行を取りだし、結果を data 類似のデータフレーム data2 にする。久保さんの解はブログ参照。他のエレガントな解を求む。なお、unique(data$group) を sort(unique(data$group)) に変えれば行はアルファベット順に並ぶ。

> data <- data.frame(group=sample(letters[1:3], 10, rep=TRUE), 
                     x=sample(1:10, 10, rep=TRUE), 
                     y=runif(10))
> data
   group x          y
1      a 5 0.11266816
2      c 9 0.36755007
3      b 4 0.01489094
4      b 8 0.84579738
5      c 1 0.43642891
6      a 6 0.50199929
7      a 8 0.12206555
8      c 5 0.49085408
9      a 9 0.29270500
10     b 5 0.37572725
> data2 <- NULL
# lapply 関数内部から data2 に永続付値し、lapply 関数自体の返り値は無視
> lapply(unique(data$group), 
         function(g) {
           temp <- data[data$group == g,]
           data2 <<- rbind(data2, temp[which.max(temp$y),]) } 
        )
> data2
  group x         y
6     a 6 0.5019993
8     c 5 0.4908541
4     b 8 0.8457974

より巨大なデータフレームの場合。

> data <- data.frame(group=sample(letters[1:10], 1e4, rep=TRUE), 
                     x=sample(1:10, 1e4, rep=TRUE), 
                     y=runif(1e4))
> data2 <- NULL
> lapply(unique(data$group), 
         function(g) {
           temp <- data[data$group == g,]
           data2 <<- rbind(data2, temp[which.max(temp$y),]) } 
        )
> data2
     group  x         y
2955     i  8 0.9994580
6527     f  1 0.9988888
6102     b 10 0.9987106
159      e  3 0.9999912
4747     c  7 0.9992973
5063     j  3 0.9996161
4429     g  1 0.9983561
9113     h  8 0.9997426
9473     d  4 0.9992060
6208     a  4 0.9970925

===== どういうのをエレガントと言えばよいのかわからないが,トリック的なものを

> set.seed(12345) # 追試のために
> data <- data.frame(group=sample(letters[1:10], 1e4, rep=TRUE), 
+                       x=sample(1:10, 1e4, rep=TRUE), 
+                       y=runif(1e4))
> # (106) の例解
> data2 <- NULL
> lapply(unique(data$group), 
+           function(g) {
+             temp <- data[data$group == g,]
+             data2 <<- rbind(data2, temp[which.max(temp$y),]) } 
+          )
途中省略
> data2
     group  x         y
793      h  4 0.9999448
4294     i  5 0.9996329
6025     e  7 0.9997141
9697     b 10 0.9990397
6514     d  4 0.9995094
4792     f  6 0.9983010
4348     j  3 0.9998260
4153     a  7 0.9974178
4219     g  3 0.9988353
8096     c  4 0.9997417
> # 別解
> g <- unique(data$group) # 要素の種類のすべてについて
> o <- sapply(g, function(i) which.max(data$y[data$group == i])) # 要素の種類ごとに y の最大値の位置
> y2 <- diag(sapply(g, function(i) data$y[data$group == i][o])) # その最大値
> x2 <- diag(sapply(g, function(i) data$x[data$group == i][o])) # そのときの X の値
> data.frame(group=g, x=x2, y=y2) # g,x2,y2 をデータフレームにする
   group  x         y
1      h  4 0.9999448   # 上の解法による結果と同じなのでホッとした
2      i  5 0.9996329
3      e  7 0.9997141
4      b 10 0.9990397
5      d  4 0.9995094
6      f  6 0.9983010
7      j  3 0.9998260
8      a  7 0.9974178
9      g  3 0.9988353
10     c  4 0.9997417

別解

by(data,data$group,function(i)data2<<-rbind(data2,i[which.max(i$y),]))

でも永続付置は邪道

別解お見事。永続付値が気持ち悪ければ tempenv <- new.env() で自前の専用環境を作り、そこへ, assign("data2", NULL, env=tempenv), get("temp2", env=tempenv) 等とすれば、気持ち悪くない(けど面倒くさい)?

別解の修正(副作用がない関数プログラミング版)。 list->data.frameの変換にはdo.callを使えばいいということに気づいたので。

do.call("rbind",by(data,data$group,function(i)i[which.max(i$y),]))

ただ、rownamesが変わっちゃうのが嫌な感じ。

(105) R の文字列ベクトルは実は整数値ベクトル!? (2007.1.21)

次のながぁーい文字列三種類を10万個並べたベクトルの保管メモリ量を計ってみると、意外なことに整数10万個を並べたベクトルのそれと実質同じ(成分あたり約4バイト)。ということは、R の文字列ベクトルは内部的には密かに因子として表現されている!? つまり三種類の文字列の何番目かという整数値(as.integer(f) で得られる)で表現されているらしい。

> x <- sample(c("abcdefghijkl", "ijklmnopqrst", "qrstuvwxyz1"), 1e5, rep=TRUE)
> f <- factor(x)
> object.size(x)
[1] 400144
> object.size(f)
[1] 400360
> object.size(1:1e5)
[1] 400024
> levels(f)
[1] "abcdefghijkl" "ijklmnopqrst" "qrstuvwxyz1"
> object.size(x)/1e5
[1] 4.00144
> object.size(f)/1e5
[1] 4.0036
> object.size(as.integer(f))/1e5
[1] 4.00024

===== っていうか,logical, character などは,実際の値が保管されているアドレス(ポインター)が保管されている?

> object.size(rep(TRUE,1e5))
[1] 400024
> object.size(rep(pi,1e5))
[1] 800024
> object.size(rep(3.14,1e5))
[1] 800024
> object.size(rep("a",1e5))
[1] 400056
> object.size(rep("abcdef",1e5))
[1] 400056

==== object.sizeは概算見積.

> A=sample(c("ABCDEFGH","12345678","abcdefgh"),1e5,replace=T)
> length(serialize(A,con=NULL))
[1] 1600022
> object.size(A)
[1] 400144
> B<-data.frame(A=sample(c("ABCDEFGH","12345678","abcdefgh"),1e5,replace=T))
> length(serialize(B,con=NULL))
[1] 400272
> object.size(B)
[1] 400728

書いては無いけど,実際のメモリー使用量ならserializeするのが一番近いと思われます.
==== (最初のコメントへのコメント)なるほどその可能性がありますね。整数の4バイトでなく、ポインタの4バイトですか。help(factor) を読んでいたら、文字列ベクトルは factor(の水準)で表現するとメモり使用量が減る、と書いてあったので、上のようなことかなと思いました。すると、因子表現も実態は水準文字列へのポインタ表現、ということなのでしょうか。
==== 追加コメント。help(scan) より。やはり、R の文字列ベクトルには文字列をそのまま記録するのではなく、共通部分集合を参照することにより、メモりを節約する工夫がされているようです。

    'scan' attempts to share storage with character strings that have
    already been read in the call. If an upper bound on the number of
    character strings cannot be deduced from 'nmax' or 'n', sharing is
    used for the first 10000 unique strings which are read in.

(104) カテゴリの re-order(中澤さんHPより)

「ロジスティック回帰で独立変数に3水準以上のカテゴリ変数があるとき,そのまま投入すると1番目のカテゴリがreferenceになるのだが,2番目以降のカテゴリをreferenceにしたい場合がある」場合の話。例えば

> warpbreaks$tension
 [1] L L L L L L L L L M M M M M M M M M H H H H H H H H H L L L L L L L L L M M
[39] M M M M M M M H H H H H H H H H
Levels: L M H
> summary(lm(breaks ~ wool + tension, data=warpbreaks))

Call:
lm(formula = breaks ~ wool + tension, data = warpbreaks)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.500  -8.083  -2.139   6.472  30.722 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   39.278      3.162  12.423  < 2e-16 ***
woolB         -5.778      3.162  -1.827 0.073614 .  
tensionM     -10.000      3.872  -2.582 0.012787 *  
tensionH     -14.722      3.872  -3.802 0.000391 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 11.62 on 50 degrees of freedom
Multiple R-Squared: 0.2691,     Adjusted R-squared: 0.2253 
F-statistic: 6.138 on 3 and 50 DF,  p-value: 0.001230 

となっているが、tensionの "M" を reference にしたい場合は

> warpbreaks$tension2 <- relevel(warpbreaks$tension, ref="M")
> warpbreaks$tension2
 [1] L L L L L L L L L M M M M M M M M M H H H H H H H H H L L L L L L L L L M M
[39] M M M M M M M H H H H H H H H H
Levels: M L H
> summary(lm(breaks ~ wool + tension2, data=warpbreaks))

Call:
lm(formula = breaks ~ wool + tension2, data = warpbreaks)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.500  -8.083  -2.139   6.472  30.722 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   29.278      3.162   9.260    2e-12 ***
woolB         -5.778      3.162  -1.827   0.0736 .  
tension2L     10.000      3.872   2.582   0.0128 *  
tension2H     -4.722      3.872  -1.219   0.2284    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 11.62 on 50 degrees of freedom
Multiple R-Squared: 0.2691,     Adjusted R-squared: 0.2253 
F-statistic: 6.138 on 3 and 50 DF,  p-value: 0.001230 

と、関数 relevel() を使用すればよい。ちなみに、関数 factor() でも可。

> warpbreaks$tension3 <- factor(warpbreaks$tension, levels = c("M", "L", "H"))
> warpbreaks$tension3
 [1] L L L L L L L L L M M M M M M M M M H H H H H H H H H L L L L L L L L L M M
[39] M M M M M M M H H H H H H H H H
Levels: M L H
> summary(lm(breaks ~ wool + tension3, data=warpbreaks))

Call:
lm(formula = breaks ~ wool + tension3, data = warpbreaks)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.500  -8.083  -2.139   6.472  30.722 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   29.278      3.162   9.260    2e-12 ***
woolB         -5.778      3.162  -1.827   0.0736 .  
tension3L     10.000      3.872   2.582   0.0128 *  
tension3H     -4.722      3.872  -1.219   0.2284    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 11.62 on 50 degrees of freedom
Multiple R-Squared: 0.2691,     Adjusted R-squared: 0.2253 
F-statistic: 6.138 on 3 and 50 DF,  p-value: 0.001230

(103) Rのgmpパッケージによる任意精度数演算 (2007-01-18)

中澤先生の所から, そう言うことも出来ると言うことで双子素数の計算

> library(gmp)
> "+"("*"(2003663613, "^"(as.bigz(2),195000 )),1)
> "-"("*"(2003663613, "^"(as.bigz(2),195000 )),1)

これではわからん! と言う人は

> 2003663613 * as.bigz(2)^195000 + 1
> 2003663613 * as.bigz(2)^195000 - 1

(102) 何度も attach する時はご用心 (2007.1.13)

リスト・データフレームを attach すると、それ(のコピー)が環境として検索パスに登録され、スコープ規則を使い、その(名前つき)変数にアクセスできる。しかし、同じ変数名成分をもつ複数のリストを attach すると?

> x <- list(a=1:3, b=3:5)
> y <- list(a=runif(4), c=letters[3:5])
> attach(x)
> a
[1] 1 2 3
> attach(y) # 既に同じ名前の変数があるよと警告

        The following object(s) are masked from x :

         a

> a         # y の成分にアクセス
[1] 0.0188327 0.3558205 0.1908796 0.9801353
> search()  # y は検索パス上 x より大局的環境に近い位置だから 
 [1] ".GlobalEnv"        "y"                 "x"
 [4] "package:stats"     "package:graphics"  "package:grDevices"
 [7] "package:utils"     "package:datasets"  "package:methods"
[10] "Autoloads"         "package:base"
> detach(y)
> a         # y を detach したので、x の成分が見えるようになった 
[1] 1 2 3
> detach(x)
> a         # もう a という名前の変数はない
エラー:オブジェクト "a" は存在しません
> a <- NA   # もし大局的環境に既に変数 a があったら?
> attach(x) # やはり警告、しかしメッセージがさっきと微妙に違う
 
       The following object(s) are masked _by_ .GlobalEnv :
 
        a
# 大局的環境中の変数は隠蔽されない! (環境 x は常に検索パスで大局的環境より遠い位置におかれるから)
> a         
[1] NA
> search()
 [1] ".GlobalEnv"        "x"                 "package:stats"
 [4] "package:graphics"  "package:grDevices" "package:utils"
 [7] "package:datasets"  "package:methods"   "Autoloads"
[10] "package:base"

(101) リストに対する二重鈎括弧演算子にベクトル添字を与える (2007.1.13)

リストの成分を取り出す二重鈎括弧演算子にベクトル添字を与えることなど考えたこともなかったが、じつはとても便利な機能がある。(ヘルプ文章を良く読むと確かにそう書いてある。)

## x[[c(i,j,k...)]] は x[[1]][[j]][[k]]... という意味
> x <- list(a=1:5, b=list(c=5:6, d=letters[1:5]), e= runif(5))
> x[c(1,2)]  # 一重鈎括弧演算子に対するベクトル添字は範囲指定の意味
$a
[1] 1 2 3 4 5
$b
$b$c
[1] 5 6
$b$d
[1] "a" "b" "c" "d" "e"
> x[[c(1,2)]]     # x[[1]][[2]] つまり x[[1]][2] と同じ
[1] 2
> x[[c(2,2,3)]]   # x[[2]][[2]][[3]] つまり x[[2]][[2]][3] と同じ
[1] "c"

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1719d)