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

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

注意:(101)以降は 知っているといつか役にたつ(?)関数達(2) へ続く。


(100) (百回記念、というほどのものではないが) 北大久保さんのブログで見かけた問題の解の例 (2007.1.11)

rep 関数のベクトル化機能と、(あまり使い道がわからなかった) mapply 関数(apply 関数の多変量版)を使った例を紹介する(このページの (12)、数列の並行生成参照)。ある種のパッケージの関数は 0,1 型データを前提としているそうで、次のようなデータフレーム

  a b c
1 2 5 a    # 5回の内1が2回
2 0 5 b    # 5回の内1が0回 
3 5 6 c    # 6回の内1が5回

を、次のような 0,1 型データフレームにわざわざ直す必要があるらしい。

   a b c
1  1 5 a
2  1 5 a
3  0 5 a
4  0 5 a
5  0 5 a
6  0 5 b
7  0 5 b
8  0 5 b
9  0 5 b
10 0 5 b
11 1 6 c
12 1 6 c
13 1 6 c
14 1 6 c
15 1 6 c
16 0 6 c

以下は rep 関数と、mapply 関数を使った例。

temp <- function() {
    data <- data.frame(a=c(2,0,5), b=c(5,5,6), c=letters[1:3]) # テスト例
    Data <- data[rep(1:nrow(data), data$b),]
    Data$a <- ifelse(unlist(mapply(seq, data$a, 1 + data$a -data$b)) > 0, 1, 0)
    rownames(Data) <- 1:nrow(Data) # Data の列ラベルを通し番号に付け替え
  } 

何をしているか段階を追って解説:

> rep(1:3, c(5,5,6)) # ベクトル引数に注意、1を5回、2を5回、3を6回繰り返す
 [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 3
> data[rep(1:3, c(5,5,6)),]
    a b c
1   2 5 a
1.1 2 5 a
1.2 2 5 a
1.3 2 5 a
1.4 2 5 a
2   0 5 b
2.1 0 5 b
2.2 0 5 b
2.3 0 5 b
2.4 0 5 b
3   5 6 c
3.1 5 6 c
3.2 5 6 c
3.3 5 6 c
3.4 5 6 c
3.5 5 6 c
> mapply(seq, data$a, 1 + data$a -data$b)
[[1]]
[1]  2  1  0 -1 -2
[[2]]
[1]  0 -1 -2 -3 -4
[[3]]
[1] 5 4 3 2 1 0
> unlist(mapply(seq, data$a, 1 + data$a -data$b))
 [1]  2  1  0 -1 -2  0 -1 -2 -3 -4  5  4  3  2  1  0
> ifelse(unlist(mapply(seq, data$a, 1 + data$a -data$b)) > 0, 1, 0)
 [1] 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 0

===~ 久保さんご本人からのコメント。上の例の第二行は

 Data$a <- unlist(mapply(seq, data$a, 1 + data$a -data$b)) > 0

でOK. TRUE, FALSE じゃ見栄えが悪い(別に問題無いが)という人はつぎのようにする。TRUE,FALSE は0を加えるとそれぞれ1,0にかわることを注意。もちろん as.numeric() で数値に強制変換しても可。

 Data$a <- 0 + (unlist(mapply(seq, data$a, 1 + data$a -data$b)) > 0) 

== 別解 == (匿名子投稿)
rep を2重に使う(今回は当てはまらないが,each 引数が有効に働くこともある)

func1 <- function()
{
data <- data.frame(a=c(2,0,5), b=c(5,5,6), c=letters[1:3]) # テスト例
Data <- data[rep(1:nrow(data), data$b),]
Data$a <- ifelse(unlist(mapply(seq, data$a, 1 + data$a -data$b)) > 0, 1, 0)
dimnames(Data)[[1]] <- 1:nrow(Data)
}
func2 <- function()
{
data <- data.frame(a=c(2,0,5), b=c(5,5,6), c=letters[1:3]) # テスト例
Data <- data[rep(1:nrow(data), data$b),]	# 百回記念を参考に
Data$a <- rep(rep(1:0, nrow(data)), rbind(data$a,  data$b-data$a)) # 1列目だけを置き換える(rbind は中澤さんの案を拝借)
rownames(Data) <- 1:nrow(Data)
}
===== 実行結果 =====
> system.time(for (i in 1:10^4) func1())
   user  system elapsed 
 29.567   0.245  29.593 
> system.time(for (i in 1:10^4) func2())
   user  system elapsed 
 24.954   0.230  25.047 
data の行数や data$b が大きくなると,両者の時間差はだんだん大きくなる

お節介にも、別解を解説(なるほど):

> rbind(data$a,  data$b-data$a)
      [,1] [,2] [,3]
 [1,]    2    0    5
 [2,]    3    5    1
> rep(1:0, nrow(data))
 [1] 1 0 1 0 1 0
> rep(c(1,0,1,0,1,0), c(2,3,0,5,5,1))
[1] 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 0
# func2 の第3行は結局こういう意味(行列はベクトル化されて使われる)
> Data$a <- rep(c(1,0,1,0,1,0), c(2,3,0,5,5,1))

(99) リストとしての環境 (2007.1.6)

環境はリストとしてアクセス出来る! おそらく内部的にはリストとして構成されている?

> e1 <- new.env() # 新しい環境
> e1$a <- 1       # リスト風操作
> e1$b <- 2
> e1              # これでは中身は見えない
<environment: 0x848ca40>
> ls(env=e1)      # 環境中の変数を見る、変数 a,b が作られている
[1] "a" "b" 
> as.list(e1)     # リストに強制変換
$b                # 変数は成分とされる
[1] 2
$a
[1] 1
> assign("c", 3, env=e1) # 環境への変数の登録の正式な仕方
> as.list(e1)
$c
[1] 3
$b
[1] 2
$a
[1] 1
> as.list(e1)[[2]] # 正式の仕方は get("b", env=e1)
[1] 2
> as.list(e1)[2]   # 成分をリストとして得る
$b
[1] 2
 
> as.list(.GlobalEnv)  # 大局的環境をリスト化!
$e1                       # e1 は大局的環境中の環境オブジェクト
<environment: 0x848ca40>
$x                         # たまたま存在したベクトル
[1]  2 -3 -2 -3  3  0  0 -2  2 -2
> .GlobalEnv$y <- 1:4  # y <- 1:4  のまわりくどい仕方 
> y
[1] 1 2 3 4

(98) sink 関数のオプション split=TRUE (2007.1.4 2010.10.02追記)

Unix では R | tee R.log とすると、R セッションのコンソール出力を全部ファイル R.log に記録できます。これが使えない(未確認) Windows ユーザーは可哀想だなと思い続けてきました。ところがしかし、sink("R.log", split=TRUE) とすると、コンソールへの出力はそのままに、ファイル R.log へも逐一記録されることを今日知りました。ただしエラー・警告メッセージ出力は駄目のようです。お試しあれ。一度 sink() で中断しても sink("R.log", split=TRUE, append=TRUE) とすれば再開追加記録できますから、tee 命令よりも柔軟です。
失礼、キーボード入力はファイルに記録されませんから、やはり tee 命令の方が便利です。何とかならないものだろうか?

Mac OS X だけかも知れませんが,Mac では,コンソールウインドウのキャパシティーは相当大きいらしく,かなりの仕事をした後でも,先頭から残っています。
「ファイル」-->「保存」または「別名で保存」することにより,キー入力と結果出力の完全な記録が残ります。

いやいや,Windows だって,同じくです。「ファイル」-->「ファイルへ保存」でコンソールに表示されたものは全部ファイルへ保存できます。 追記 実際に試しましたがWindowsのGUI版RのキャパシティーはMac版と比べるとかなり小さいです。キャパシティーを変更するコマンドがあるのかは知りません-のの

いやいや, Unix(以外でもいいけど) emacs上でessを使えば, C-x C-w でバッファ上のデータは全部ファイルへ保存出来ます. (新春大喜利?)

理想的なのは、入出力・メッセージすべてを、同じファイルに、必要なところだけを見繕って保存出来るような機能(出来ればR内部から操作できる)なんですが、心あたりありませんか(要望掲示版向きになってしまいそうな)。
追記 一番可能性があるのはTeachingDemosパッケージのtxtStart()です。こちらはWindowsでは上手く動きますが、MacのGUI版Rに問題があるらしく、Macでは上手く動きません。ただし、ターミナル上のRからは上手く動きます。 -のの

(97) cat 関数の便利なオプション (要望掲示版記事より, 2007.1.3)

コンソールへの出力に良く用いられる cat 関数にはいくつかの便利なオプションがある。 sep = "," とすれば、出力項目間に分離記号 , をはさむことができる。既定の fill=FALSE では改行文字 "\n" を置かないと改行され無いが、改行文字 "\n" を置くと暗黙のうちに最後に空白が一つ置かれる。fill = 30 等とすれば、出力は幅30文字毎に改行され行末に不要な空白は置かれない。fill = TRUE ならば、改行幅はシステムオプションである options()$width (普通80)文字毎に改行される。labels オプションで出力各行の先頭に順に置かれるラベル文字文字列ベクトルを指定できる。

> cat("iteration = ", 9, "\n")
iteration =  9
> cat(letters, fill=30)  # 出力幅30文字で改行
a b c d e f g h i j k l m n o
p q r s t u v w x y z
> cat(letters, fill=40)  # 出力幅40文字で改行
a b c d e f g h i j k l m n o p q r s t
u v w x y z
# 既定の出力幅80文字(最後に一つの空白が追加される)
> cat(letters, "\n")
a b c d e f g h i j k l m n o p q r s t u v w x y z
> cat(letters[1:5], "\n", sep=",")      # 分離記号 "," をはさむ
a,b,c,d,e,  # 既定(fill=FALSE)では行末に空白が加わるため最後にもカンマ
> cat(letters[1:5], fill=TRUE, sep=",") # fill=TRUE とすれば良い
a,b,c,d,e
# 行頭ラベル、実際の行数しか使われないことを注意
> cat(paste(letters, 100 * 1:26), fill = TRUE,
      labels = paste("{", 1:10, "}:", sep = ""))
{1}: a 100 b 200 c 300 d 400 e 500 f 600 g 700 h 800 i 900 j 1000 k 1100 l 1200
{2}: m 1300 n 1400 o 1500 p 1600 q 1700 r 1800 s 1900 t 2000 u 2100 v 2200
{3}: w 2300 x 2400 y 2500 z 2600

(96) リストも attach, detach できる (2007.1.2)

データフレームを attach すると、その成分変数を変数名だけでアクセスできることは周知の事実。しかし、実は(すべての成分が名前を持つ)リストも attach できる。データフレームは行列風の外観を持つ(そしてすべての成分が名前を持つ)リストの一種だから当然といえば当然かも知れないが。

> foo <- list(A=1:4, B=letters[1:4])
> foo$A      # 通常のアクセス法
[1] 1 2 3 4
> A 
エラー:オブジェクト "A" は存在しません
> attach(foo) # foo を atatch する
> A            # 変数名 A だけでアクセスできる   
[1] 1 2 3 4
> B
[1] "a" "b" "c" "d"
> detach(foo)  # detach する
> A
エラー:オブジェクト "A" は存在しません
> B
エラー:オブジェクト "B" は存在しません

? attach によれば,第一引数は,“database”. This may currently be a data.frame or a list or a R data file created with save or NULL or an environment.
save で保存したデータファイルの使用例は以下のようなものか。使い道は…。巨大なデータベースなどは,メモリーに読み込まなくて使えるとか??

> data <- data.frame(x=1:5, y=rnorm(5))
> save(data, file="saved.data")
> rm(data)
> data$x
NULL
> attach("saved.data")
> data$x                # data は明示的には読み込んでいない
[1] 1 2 3 4 5
> data
  x           y
1 1 -1.23017225
2 2 -1.32405869
3 3  1.26124227
4 4  1.31923172
5 5 -0.08075376
> data$y[3]             # 要素も引き出せる
[1] 1.261242

(95) コンソール出力の抑制 (2006.12.28)

Q&A(初級)の中間さんのコメントで知った関数 suppressWarnings() は警告のコンソールへの出力を抑制する。確かに警告ばかり出されては不快になるであろうからこうした関数を用意したくなる気持ちはわかる(嘘)。同様に関数の最後に return 文の代わりに使われる「コンソール出力無しの返り値関数」invisible() もコンソールへの自動出力を行なう命令と一緒に使えば気分が良い、例 invisible(女房の小言())。また try 関数もオプション付きで try(EXPR, TRUE) とすると途中の警告メッセージの出力を抑制してくれる。

# 警告メッセージの出力を消す
> suppressWarnings(log((-1):2))
[1]       NaN      -Inf 0.0000000 0.6931472
> warnings() # 警告は後から見ることができる 
警告メッセージ: 
計算結果が NaN になりました in: log(x)
# 不要な出力を消す(t <- gc() とダミー変数に付値しても良い)
> gc()
         used (Mb) gc trigger (Mb) max used  (Mb)
Ncells 312298  8.4     597831 16.0   467875  12.5
Vcells 361111  2.8     857161  6.6 74307219 567.0
> invisible({gc(); gc()}) # 結果出力が消える

(94) gc() は二度繰り返すべし。(2006.12.27)

「Rの基礎とプログラミング技法」、U. Ligges 著、石田基広訳、によるとガベージコレクション関数 gc() は二度続けて行なう方が良いそうである。理由はもっともで、直近の計算結果を保存する .Last.value というオブジェクトが気づかぬままに巨大になっていることがあるが、一度目の gc() でこれが消え、二度目の gc() 実行で解放されたメモリーが整頓される、からだそうです。効果はしばしば劇的だそうですから、お試しあれ。

(93) 関数!、関数!、関数! (2006.12.27)

Rは関数型言語だというのは知ってはいたが、ここまでとは今の今まで気づかなかった。

> "+"(1,3)
[1] 4
> "*"(1,3)
[1] 3
> "^"(4,3)
[1] 64
> "/"(4,3)
[1] 1.333333
> "-"(4,3)
[1] 1

ここまでは知る人ぞ知る二項演算子型関数の別の表記。以下は私には驚き(先刻承知という人はえらい)
 
> "if"(TRUE,3)             # if 文は実は if 関数
[1] 3
> "if"(FALSE, 3,4)
[1] 4
> ":"(1,4)                 # 等差数列演算子
[1] 1 2 3 4
> "["(matrix(1:4,2,2),1,1) # 行列に対する括弧演算子
[1] 1
> "("(exp(1))              # これが例の計算結果をすぐに表示する () 記法の正体(?)
[1] 2.718282              
> "<-"(x,4)                # 付値演算子もこの通り
> x
[1] 4
> "="(x, NA)               # もう一つの付値演算子 x = NA 
> x
[1] NA
> (x <- for(i in 1:4) i)   # for 文だって実は値を返す立派な関数(?)
[1] 4

他にも「こんなのあるぞ」という方は教えて下さい(知ってもどういうことはないけれど)。


(92) 二つの多角形の交差部分の面積を求める、パッケージ gpclib (2006.11.18)

北大の久保さんのぎょーむ日誌より C&P。ちなみに、このブログは R の使い方の参考としても一級のソース。私見では、久保さんは日本一の R へビーユーザー。

library(gpclib)
x1 <- c(0, 2, 2, 0) # 第一の多角形(正方形)の四隅のx座標
y1 <- c(0, 0, 2, 2) # 同じくy座標
x2 <- c(1, 3, 3, 1) # 第二の多角形(正方形)の四隅の座標
y2 <- c(1, 1, 3, 3) # 同じくy座標
xy1 <- structure(c(x1, y1), .Dim = c(length(x1), 2))     # 多角形頂点座標からなる行列
xy2 <- structure(c(x2, y2), .Dim = c(length(x2), 2))
> xy1                                                    # こんな感じ
     [,1] [,2]
[1,]    0    0
[2,]    2    0
[3,]    2    2
[4,]    0    2
p1 <- as(xy1, "gpc.poly")                                # 多角形オブジェクトに変換
p2 <- as(xy2, "gpc.poly")
plot(append.poly(p1, p2))                                # 二つ併せてプロット
cat("# area of p1 =", area.poly(p1), "\n")
p12 <- intersect(p1, p2)                                 # 二つの多角形の交差部分からなる多角形
plot(p12, poly.args = list(col = "#ff8000"), add = TRUE) # 交差部分を色付で重ね描き
cat("# area of p12 =", area.poly(p12), "\n")             # area.poly() で交差部分の面積計算

(91) 環境および環境内の変数のロック (2006.11.18)

R掲示板で既に話題になった R.2.4 で登場の新機能。環境および環境内の変数(とその値を)ロックし、改変できないようにする。本家 S-plus に勝る R の長所が複数の環境(オブジェクトとその値を記録する領域)を持てること。ただし、逆にどの環境のオブジェクトを操作してるかを常に意識していないと、面倒なことにもなる。ロック機能は不注意によるオブジェクトの誤操作を防ぐ。以下はヘルプドキュメント:

bindenv                 package:base                 R Documentation

Binding and Environment Adjustments 

Description:

   これらの関数は環境と、環境の調整と、環境内のバインディング(拘束)に対する 
   実験的なインタフェイスである。環境や個々の拘束をロックしたり、変数をある関
   数にリンクすることを可能にする。  

Usage:

    lockEnvironment(env, bindings = FALSE)
    environmentIsLocked(env)
    lockBinding(sym, env)
    unlockBinding(sym, env)
    bindingIsLocked(sym, env)
    makeActiveBinding(sym, fun, env)
    bindingIsActive(sym, env)

Arguments: 

    拘束が既にロックされていない限り、'lockBinding' は特定の環境中の拘束を
    ロックする。ロックされた拘束の値は変更できない。しかし、環境がロックされ
    ていない限り、ロックされた拘束は環境から取り除くことができる。

    'makeActiveBinding' は、'sym' が引数無しの 'fun' の呼出しになり、'sym' へ
    の値の付値が、その値を単独の引数とする'fun' の呼出しになるように、関数 
    'fun' を定義する。これは、R 変数にリンクされた C 変数や、データベースに
    リンクされた変数、といった物の定義を可能にする。これはまた、システムの大域
    オブジェクトの thread-safe なバージョンを作るのにも有益である。 ~

以下は example(bindenv) の実行内容:

> e <- new.env()           # 新しい環境 e を作る
> assign("x", 1, env = e)  # 環境 e に変数 x を値 1 で拘束
> get("x", env = e)        #
[1] 1
> lockEnvironment(e)       # 環境 e をロック
NULL
> get("x", env = e)     
[1] 1
> assign("x", 2, env = e)  # すでに拘束されている変数の値は変えられる
> try(assign("y", 2, env = e)) # 新しい変数は拘束できない
以下にエラー assign("y", 2, env = e) : ロックされた環境にバインディングを
追加することはできません
> e <- new.env()           # 新しい環境 e
> assign("x", 1, env = e)  # 変数 x を環境 e に値 1 で拘束
> lockBinding("x", e)      # 変数 x をロック
NULL
> try(assign("x", 2, env = e)) # ロックされた変数の値は変更できない
以下にエラー assign("x", 2, env = e) : 
ロックされたバインディングの値は変更できません
> unlockBinding("x", e)     # 変数 x のロック解除
NULL
> assign("x", 2, env = e)   # 値の変更ができるようになる
> get("x", env = e)
[1] 2

次の機能は面白そうだが、本来の使い道は?

> f <- local({  # 関数環境に拘束された変数 x を持つ関数 f を定義
              x <- 1
              function(v) {
                if (missing(v))
                  cat("get\n")
                else {
                  cat("set\n")
                  x <<- v
                }
              x
             }
       })
> makeActiveBinding("fred", f, .GlobalEnv)  # "fred" を関数 f の同義語(synonym)とする
NULL
> bindingIsActive("fred", .GlobalEnv)
[1] TRUE
> fred  # 引数無しの呼出し f() と同効果  
get
[1] 1
> fred <- 2  # 呼出し f(2) と同効果
set
[1] 2
> fred  # f 環境に拘束された変数 x の値が 2 に変更されていることに注意
get
[1] 2

(90) R session を個別に保存・回復する、パッケージ session (2006.11.18)

パッケージ session は現在の R セッションを指定したファイルに保存し、回復する関数 save.session(), restore.session() を提供する。同様の機能(ファイル .Rdata に自動的に保存される)と似ているが、特定のセッション毎に別のファイルに保存でき、セッション中にロードしたライブラリ情報も記録されるので、便利かも。以下は別々のセッションを save.session("RSession.1"), save.session("RSession.2") で保存し、それを回復する例。画像情報等は記録できない。

> rm(list=ls()) # 掃除
> library(session) # パッケージ session をロード
> restore.session("RSession.1") # RSession.1 を回復
Loading all data...             # データを復帰
Loading packages...             # 使用されたライブラリも復帰
Restoring search path...        # 検索パスも復帰
Done.
> ls()
 [1] "birthwt.glm"   "birthwt.step"  "birthwt.step2" "bwt"
 [5] "cpus.lm"       "cpus.lm2"      "cpus.samp"     "cpus0"
 [9] "cpus1"         "ftv"           "ptd"           "quine.hi"
[13] "quine.nb"      "quine.nb2"     "quine.nxt"     "quine.stp"
[17] "race"          "v"             "x"
> rm(list=ls()) # 再び掃除
> restore.session("RSession.2") # 今度は別のセッション記録を復帰
Loading all data...
Loading packages...
Restoring search path...
Done.
> ls()
[1] "Iris"  "iris3" "train" "z"     "z1"

(89) リストの逆数の和

以下のようなリストを作成し、

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

a1 と a2 ごとに逆数の和を求めたい。sapplyを使えばいいとは思ったが、逆数和を求める関数がない。自分で関数を定義するしかないのか?
解1 個別に定義 …… 実は,リストの要素はたくさんあるので。。。

c(sum(1/a[[1]]),sum(1/a[[2]]))

解2

lapply(lapply(a, function(x) 1/x),sum)

解3 こっちの方がいいか

mapply(function(x) sum(1/x), a) 

実行例

>  a <- as.list(NULL)
>  a <- c(list(1:3), list(4:9))
>  c(sum(1/a[[1]]),sum(1/a[[2]]))
[1] 1.833333 0.995635
>  lapply(lapply(a, function(x) 1/x),sum)
[[1]]
[1] 1.833333

[[2]]
[1] 0.995635

>  mapply(function(x) sum(1/x), a) 
[1] 1.833333 0.995635

(88) Unix 風の(ファイル名の)ワイルドカード指定を正規表現に変換する関数 glob2rx (初級Q&A より、2006.11.07)

次の実行例は Rplot の後に数字がならび、最後が .jpg で終る文字列に対する R 風正規表現を得る。

> glob2rx("Rplot[0-9]*.jpg")
[1] "^Rplot[0-9].*\\.jpg$"  # 得られた正規表現。^ と $ はそれぞれ先頭・末尾を表す

応用例。作業ディレクトリにあるそうした名前のファイルの中から最新のタイムスタンプ (生成時刻)を持つものの名前を取り出す。解説すると、glob2rx で Unix 風のファイル名のワイルドカード指定を正規表現に変換し、list.files でそれに合致するファイル名の情報からなるデータフレームを得、その生成時刻列のみを ["mtime"] で取り出し a に付値、一番最新(大)の行番号を which.max([a[,1]) で求め、さらに対応する行ラベル名を row.names ベクトル(該当ファイル名文字列ベクトル)から得ている。同様のことを行なう別の方法は元記事参照。

> row.names(a <-  file.info(list.files(pattern=glob2rx("Rplot[0-9]*.jpg")))["mtime"])[which.max(a[,1])]
[1] "Rplot002.jpg"

(感想:R ユーザーの広がりを反映してか、最近 Q&A コーナにも答がいのある質問が増えてきましたね。)

(87) 既存関数の一部の仮引数の既定値を変更した関数を定義する (初級 Q&A より、2006.11.05)

pie 関数の引数 clockwise の既定値 FALSE (反時計回り)を TRUE (時計回り)に変更した関数を定義するには、という質問に対する回答。
面倒だが次のようにすれば、確実:

> pie2 <- function(x, labels = names(x), edges = 200, radius = 0.8,
                   clockwise = TRUE, init.angle = if(clockwise) 90 else 0,
                   density = NULL, angle = 45, col = NULL, border = NULL,
                   lty = NULL, main = NULL, ...)
          { pie(x, labels, edges, radius, clockwise, init.angle,
                density, angle, col, border, lty, main, ...)
          }

しかし、その他大勢引数 ... を使えば次のように簡単に定義できる:

pie2 <- function(...) pie(...,clockwise=TRUE) 

ただし、後の例では pie2(rep(1:10, 5), clockwise=FALSE) とするとエラーになる。最初の例は大丈夫。
更に追加コメントあり、おそらくこれが決定版。

pie2 <- function(x, clockwise=TRUE, ...) pie(x, clockwise=clockwise, ...)

(86) 関数 local (from r-help, 2006.10.22)

関数 local は(既定で)新しい環境の中で R expression を評価する。r-help に

> n <- 3
> f <- function(x) x^n
> f(2)
[1] 8
> n <-1
> f(2)
[1] 2

という例をあげ、ここで f 中の n を関数定義時の n に固定するには?、という質問があった。実際には、n は「ある複雑な計算により得られた特別な値」との解説。こうした時は local 関数を使い、関数定義時の n を関数固有の環境に固定すると良いらしい。

> n <- 2
> f <- local({nn <- n; function(x) x^nn}) # 単に x^n では駄目(親環境の n だから)
> f(2)
[1] 4
> n <- 3
> f(2)     # 依然 x^2 のまま
[1] 4
> f
function(x) x^nn
<environment: 0x8909d88>

ヘルプ文章からの抜粋:つまり変数 nn は関数 f の附属環境中に登録隠蔽されて、外からは見えなくなる(したがって変更されなくなる)

    'local' evaluates an expression in a local environment.  It is
    equivalent to 'evalq' except that its default argument creates a
    new, empty environment.  This is useful to create anonymous
    recursive functions and as a kind of limited namespace feature
    since variables defined in the environment are not visible from
    the outside.

(85) データフレームから分割表を作る (初級者Q&Aより, 2006.10.21)

以下のようなデータフレーム

year type value
2001 A    5
2002 A    10
2003 A    15
2001 B    1
2002 B    3
2003 B    5

を次のような二元分割表(まがい)に集計する方法は?

      A   B
2001  5   1
2002  10  3
2003  15  5

例えば year=2001 で type=A のものが二つ有れば例の様にはまとめられないことを注意。最初のデータフレームを x とすると y <- table(x)は year,type,value に対する標準的三元分割表を計算し、y[,,"15"] は value=15 の値に対する year,type の組合せごとの重複度を込めた二元分割度数分布テーブルを計算してくれるが、これは希望のものではないし...
この問題に対するエレガントで意表をついた回答(お見事)

> # 重複したデータがないことを、lengthで確かめて、、
> tapply(df$z,list(df$x,df$y),length)     
      A B 
 2001 1 1 
 2002 1 1 
 2003 1 1
> # sum (または mean)で集計。
> tapply(df$z,list(df$x,df$y),sum)     
       A B 
 2001  5 1 
 2002 10 3 
 2003 15 5 

第二引数がリストだと tapply 関数の結果がデータフレームに整形されて出力されるところが妙。
さらに xtabs 関数を使った正当派回答(xtabs 関数結構使えそう)

> df <- data.frame(year = c(2001, 2002, 2003, 2001, 2002, 2003),
+   type = c("A", "A", "A", "B", "B", "B"),
+   value = c(5, 10, 15, 1, 3, 5))
> xtabs(value~year+type,df)
      type
year    A  B
  2001  5  1
  2002 10  3
  2003 15  5

もし重複値があるとどうなるか。重複する値は加算されてしまう。

> df <- data.frame(year=c(2001,2002,2003,2001,2002,2003,2001),
+  type=c("A","A","A","B","B","B","A"),
+ value=c(5,10,15,1,3,5,100))
> xtabs(value~year+type,df)
      type
year     A   B
  2001 105   1
  2002  10   3
  2003  15   5

(84) Rprof (from r-help, 2006.8.22)

(81) の話題に関して次のような Rprof の使い方がでていたのでついでに紹介。Rprof は一定時間毎に使用関数をサンプリングし各関数がどれくらい時間を使用しているか報告する。以下の例では strsplit 関数が全実行時間の99%を占有していて、これがボトルネックであることがわかる(マッチングを行なっているから当然)。

> repeated.measures.columns <- paste(1:100000, collapse = ",") 
> library(utils)                # Rprof はライブラリ utils 中にある
> Rprof(tmp <- tempfile())      # Rprof 開始、結果をtmp に記録
> res1 <- as.numeric(unlist(strsplit(repeated.measures.columns, ",")))
> Rprof()                       # Rprof 終了
> summaryRprof(tmp)             # プロファイリング結果要約
$by.self                        # 関数 strsplit が実行時間の99.7%を占めている
                   self.time self.pct total.time total.pct
"strsplit"              23.68     99.7      23.68      99.7
"as.double.default"      0.06      0.3       0.06       0.3
"as.numeric"             0.00      0.0      23.74     100.0
"unlist"                 0.00      0.0      23.68      99.7 

$by.total
                   total.time total.pct self.time self.pct
"as.numeric"             23.74     100.0      0.00      0.0
"strsplit"               23.68      99.7     23.68     99.7
"unlist"                 23.68      99.7      0.00      0.0
"as.double.default"       0.06       0.3      0.06      0.3 

$sampling.time
[1] 23.74
> Rprof(tmp <- tempfile())
> res3 <- eval(parse(text=paste("c(", repeated.measures.columns, ")")))
> Rprof()

> summaryRprof(tmp)
$by.self
       self.time self.pct total.time total.pct
"parse"      0.42     87.5       0.42      87.5
"eval"       0.06     12.5       0.48     100.0 

$by.total
       total.time total.pct self.time self.pct
"eval"        0.48     100.0      0.06     12.5
"parse"       0.42      87.5      0.42     87.5

$sampling.time
[1] 0.48

(83) 最適化パッケージ trust (2006.8.21)

ついでに発見したパッケージ trust も多変数関数の極値を求める。導関数が必要。制約無し最適化だが、探索領域を指定できる(しかし境界上の最適値は見つけられるとは限らないとのこと)。導関数が利用できる場合は効率・性能が良いらしい。以下パッケージの解説引用:

Description:

     This function carries out a minimization or maximization of a
     function using a trust region algorithm.  See the references for
     details.

Usage:

     trust(objfun, parinit, rinit, rmax, parscale,
         iterlim = 100, fterm = sqrt(.Machine$double.eps),
         mterm = sqrt(.Machine$double.eps),
         minimize = TRUE, blather = FALSE, ...)

(82) 最適化パッケージ powell (2006.8.21)

新パッケージ powell は多変数関数の極値を求める。最適化アルゴリズムは色々ある、ということはどれも一長一短ということ。道具はたくさんあるほうが良い。このパッケージはPowell による UObyQA (Unconstrained Optimization by Quadratic Approximation)法を移植したもの。説明によると、目的関数を探索範囲で2次関数近似することにより最適値を求める方法らしい。以下パッケージの解説引用:

Description:

    Optimizes a function using Powell's UObyQA algorithm.

Usage:

     powell(par, fn, control = powell.control(), check.hessian = TRUE, ...)

Arguments:

     par: Starting values for objective function

      fn: A function to be optimized. The function takes the parameters
          ('par') as its first argument.

 control: A list of control parameters

check.hessian: logical; if 'TRUE' an eigenvalue decomposition is used
          to check the hessian for positive definiteness.

     ...: Additional arguments to be passed to 'fn'

Details:

     This function seeks the least value of a function of many
     variables, by a trust region method that forms quadratic models by
     interpolation. The algorithm is described in "UOBYQA:
     unconstrained optimization by quadratic approximation" by M.J.D.
     Powell, Report DAMTP 2000/NA14, University of Cambridge.

(81) 文字列を数値ベクトルに変換するいくつかの方法 (r-help, 2006.8.20)

文字列 "3,6,10" を数値ベクトル c(3,6,10) に変えるには、という質問に対する R guru 達の答えあれこれ.

repeated.measures.columns <- "3,6,10"

# Roger Bivand
as.numeric(strsplit(repeated.measures.columns, split = ",")[[1]])

# Peter Dalgaard & G. Grothendieck
scan(textConnection(x), sep=",")     # scan 関数の始めてみる使い方
closeAllConnections()                # textConnection を閉じる

scan(con <- textConnection(x), sep = ",") # もしくは
close(con)

# Marc Schwartz
as.numeric(unlist(strsplit(repeated.measures.columns, ",")))

# Brian D. Ripley
# a bit more general (e.g. allows spaces, works with complex numbers)

eval(parse(text=paste("c(", repeated.measures.columns, ")")))

# Marc Schwartz showed that Professor Ripley's suggestion is much faster
# than the competition with some system.time trials.

# 解説すると
# paste("c(", repeated.measures.columns, ")") で三つの文字列 "c(","3,6,10",")" をつないで文字列 "c(3,6,10)" を作り
# parse 関数でそれを R の構文として解釈し
# 最後に eval 関数でそれを評価する、という仕組み

既に,2006/08/03 の「中級Q&A」にもあったが,一番最初のと同じだけど unlist の方がスマート?

> s <- "3,6,10"
> (x <- as.numeric(unlist(strsplit(s, ","))))
[1]  3  6 10

(80) 日付を数値に変換する。(2006.08.15), 追加 (2006.11.13)

複数の機械で同じシミュレーションをする際、異なる乱数種を得るのに使用((78)参照)

> date2num <- function() {
   x <- date()
   .day    <- as.numeric(substr(x,9,10))
   .hour   <- as.numeric(substr(x,12,13))
   .minute <- as.numeric(substr(x,15,16))
   .second <- as.numeric(substr(x,18,19))
   .year   <- as.numeric(substr(x,21,25))
   return(.day*1e6+.hour*1e4+.minute*1e2+.second)
 }
> date();date2num()
[1] "Tue Aug 15 15:03:04 2006"
[1] 15150304

以下のような関数を作ってみた(2006.11.13)

date2num2 <- function() {
	as.numeric(gsub("[ :]", "", substr(date(),9,19)))
}

> date()
[1] "Mon Nov 13 18:54:06 2006"
> date2num2()
[1] 13185406
> date2num() # 最初に示された関数による結果
[1] 13185406

速度の測定

> system.time(for (i in 1:100000) date2num2())
           user          system           total   user.children system.children 
         10.276           0.095          10.295           0.000           0.000 
> system.time(for (i in 1:100000) date2num())
           user          system           total   user.children system.children 
         27.260           0.253          27.324           0.000           0.000 

(79) クラッシュに備え、一時間毎に必要な途中結果をセーブする (r-help 記事より, 2006.08.14)

John Morrow a 醇Pcrit :

>> Hello fellow R'ers, I have a simple calculation with a very large data set
>> being generated (34.9 million values) on a somewhat unreliable XP box that
>> will likely take ~ 74hrs.  I wanted to know if there is a way to have my
>> script automatically "save.image()" throughout the calculation in case of a
>> crash.  This could be on the basis of output generated or time elapsed.  I
>> checked the archive, and only got a hint of it from:
>> https://stat.ethz.ch/pipermail/r-help/1997-May/001611.html
>> Any quick suggestions would be greatly appreciated,
>> John Morrow

# make a save every hour
minutes = substr(date(),...);
if (minutes=='00') save_my_stuff;

(78) 計算が終わったら結果をメイルで送る(貧乏人の並列処理) (中級Q&A より、2006/07/21)

高価なグリッドシステムには手がでなくても、遊んでいるパソコンなら回りにいくつもあるご時世。同じことを繰り返すシミュレーションなら、複数のパソコンで分けて計算すれば真の分散処理。しかし、いちいち計算がおわったかどうか確認し、結果を手動で回収するのも、長時間かかる計算では億劫。次は、計算が終わると結果ファイルをメイルで送る工夫。但し、これは Unix でないと無理(?)。

x <- matrix(runif(100), 10,10)
dump("x",file="test.data")   # 結果をファイル test.data にかきこむ。
system("cat test.data | mail -s HereIsYourResult nanashinogonbei@foo.com")

膨大なファイルは圧縮して送りたいが、バイナリファイルはそのままではおくれないので uuencode 命令でテキスト化して送りたい(後で uudecode 命令で元に戻す)。ただし、これはもしかすると元のファイルよりも長くなるかもしれない。

system("cat test.data | gzip | uuencode test.data.gz | mail -s HereIsYourResult whoareyou@bar.com")

応用として、複数のファイルを tar 命令でまとめた上、uuencode してテキストファイルとして送ることも出来る。

system("tar czf - ./test* | uuencode comp.results.tgz | mail -s comp.results maseshigeru@gmail.com")

もっとありそうな使いかたは、計算が終わったことをメイルで知らせること。

system("echo \"Your job at Machine A is done.\" | mail -s Notice whoareyou@foo.com") 

(77) R 関数で実引数の参照渡しを実現するパッケージ proto (2006/06/23)

R 関数での引数は値渡し(関数内で操作されるのは実引数のコピー)が基本である。オブジェクトを関数で変更するには、返り値を付値するか、永続付値をつかう。ライブラリ proto は(その機能の一部として)オブジェクトの関数への参照渡し(実引数オブジェクトそのものを操作し、関数実行後はオブジェクトそのものが自動的に変更される)を実現する。巨大なオブジェクトを操作するときは有用かもしれない。R が複数の環境を持てることを利用してこうした機能を実現しているらしい(?)。

> library(proto)
> p <- proto(mat = matrix(1:25, 5))
> p # p の構造は何やら不可解
<environment: 0x85158ac>
attr(,"class")
[1] "proto"       "environment"
> str(p)
Classes 'proto', 'environment' length 3 <environment>
> p$mat  # p の mat 成分(?)にはリスト風にアクセスする
    [,1] [,2] [,3] [,4] [,5]
[1,]    1    6   11   16   21
[2,]    2    7   12   17   22
[3,]    3    8   13   18   23
[4,]    4    9   14   19   24
[5,]    5   10   15   20   25

> incr <- function(x) with(x, mat[,1] <- mat[,1] + 1) # テスト用関数、一見終了後は何も変化は無いように見える
> incr(p) # proto オブジェクト p に関数 incr を実行
> p$mat  # p の mat 成分が何時の間にか変更されている!
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    6   11   16   21
[2,]    3    7   12   17   22
[3,]    4    8   13   18   23
[4,]    5    9   14   19   24
[5,]    6   10   15   20   25
 
> x <- matrix(1:25,5) # 説明用のオブジェクト
> incrb <- function(y) y[,1] <- y[,1]+1
> incrb(x)
> x          # 関数内で操作されるのは x のコピーだから、当然 x 自体は変化しない
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    6   11   16   21
[2,]    2    7   12   17   22
[3,]    3    8   13   18   23
[4,]    4    9   14   19   24
[5,]    5   10   15   20   25
 
> incrc <- function(x) {x[,1] <- x[,1]+1; return(x)} # 変更した x のコピーを返す
> x <- incrc(x)   # x -> x のコピー -> x という順序で x が変更される
> x               # 変更後
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    6   11   16   21
[2,]    3    7   12   17   22
[3,]    4    8   13   18   23
[4,]    5    9   14   19   24
[5,]    6   10   15   20   25

(76) 可変長(その他大勢)引数の内容を取り出す(初級Q&A記事より, 2006/06/22)

その他大勢引数が実際は複数の引数になるとき、k 番目の項目を取り出すには list(...)とリスト化すれば良い。その他引数が全てスカラーなら c(...) としても良い。

> test <- function(i,...) print(list(...)[[i]])
> test(1,1,"abc",list(1:4),matrix(1:4,2,2))
[1] 1
> test(2,1,"abc",list(1:4),matrix(1:4,2,2))
[1] "abc"
> test(3,1,"abc",list(1:4),matrix(1:4,2,2))
[[1]]
[1] 1 2 3 4
> test(4,1,"abc",list(1:4),matrix(1:4,2,2))
     [,1] [,2]
[1,]    1    3
[2,]    2    4

特に ... 引数が実際は幾つの引数からなるかを知るには length(list(...)) とすれば良い。


(75) 日付関数(初級Q&A記事より, 2006/06/09)

R 及びその貢献パッケージには日付を扱う関数群が結構ある。プログラムではある起点年月日からの経過日数の差分で表す(Julian date と一般に呼ばれている)流儀が基本であるが、人間には "mm/dd/yyyy" という文字列で表すほうがわかりやすいので、両者を互いに変換する関数が必要になる。例えば R のパッケージの date には次の関数がある。

as.date                 Coerce Data to Dates
date.ddmmmyy            Format a Julian date
date.mdy                Convert from Julian Dates to Month, Day, and Year
date.mmddyy             Format a Julian date
date.mmddyyyy           Format a Julian date
is.date                 Date Objects
mdy.date                Convert to Julian Dates

一方まさに日付・時間専用のパッケージ chron には

chron                   Create a Chronological Object
cut.dates               Create a Factor from a Chron or Dates Object
dates                   Generate Dates and Times Components from Input
day.of.week             Convert between Julian and Calendar Dates
days                    Return Various Periods from a Chron or Dates Object
hours                   Return Hours, Minutes, or Seconds from a Times Object
is.holiday              Find Weekends and Holidays in a Chron or Dates Object
seq.dates               Generate Chron or Dates Sequences
trunc.times             Truncate times Objects

といった関数がある。chron の dates 関数は元祖 S の同名の関数に由来するらしい。ここで注意すべきは

> dates(0,out="yyyy/mm/dd") # chron パッケージの起点年月日
[1] 1970/Jan/01
> date.mmddyyyy(0) # date パッケージの起点年月日
[1] "1/1/1960"

と起点年月日が異なること。google で少し調べたら Julian date とは本来紀元前4713/1/1の正午を起点とする主に天文学で使われる日付単位とか(そして Julian とはユリウス暦の Julius、つまりシーザー、とは無関係)。Unix 界では1970/1/1/0:00:00 からの経過秒で表す慣習が。さらに SASでは 1960/1/1 からの経過日数を使うようで、date パッケージはこれを継承しているのか。ややこしい。Julian date を扱うときはまず起点年月日が何か確認する必要があるようだ。なお、年月日をどういう順序で並べるかも国により微妙にちがうらしい。 --


(74) 度数分布表を作る方法 (初級Q&A 記事より, 2006.04.11)

R には度数分布表を作る専用関数 table があるが、該当するデータが無い値は無視される。該当データが無いときは頻度ゼロとした度数分布ベクトルを作る専用機能はなさそうである。以下にそうした二つの比較的簡単な方法を紹介する。

> x <- sample(0:10, 20, replace=TRUE)
> x
 [1]  0  2  0 10  0  6  0  8 10  3 10  3  5  6  2  4  3  1  4  6
> table(x)  # 欠損した値 7,9 は無視される
x
 0  1  2  3  4  5  6  8 10
 4  1  2  3  2  1  3  1  3
> sapply(0:max(x), function(i) sum(x==i)) # 第一法
 [1] 4 1 2 3 2 1 3 0 1 0 3
> table(factor(x, levels=0:max(x)))       # 第二法
 0  1  2  3  4  5  6  7  8  9  10
 4  1  2  3  2  1  3  0  1  0   3
> apply(outer(x, 0:10, "=="), 2, sum)     # 第三法 いささかトリッキーだが,あえて追加
 [1] 4  1  2  3  2  1  3  0  1  0   3

これらは 0:max(x), levels=0:max(x) の部分をそれぞれ適宜変更すれば、負の値 も取る場合、非数値データの場合も使える。たとえば以下のような例。

> x <- c("a", "b", "d", "aa", "b")
> table(x)
x
 a aa  b  d 
 1  1  2  1 
> table(factor(x, levels=c("a", "aa", "b", "bcd", "c", "d")))

  a  aa   b bcd   c   d 
  1   1   2   0   0   1 

これを用いて例えば二項分布の対数尤度関数を定義すれば

N <- 20                   # データの可能な最大値 
X <- rbinom(100, N, 0.55) # 人工データ
loglikelihood <- function(p) {
                   log.prob <- dbinom(0:N, N, prob=p, log=TRUE) # 対数二項確率
                   freq <- table(factor(X, levels=0:N))
                   return(sum(freq*log.prob))
                 }

この関数の場合は,table 関数を使う必要はないので問題は簡単になる。

loglikelihood2 <- function(p) {
                    return(sum(dbinom(0:N, N, prob=p, log=TRUE)[X+1]))
                  } 

計算時間も短いようだ。

> system.time(for(i in 1:1000) loglikelihood(0.2))
[1] 4.39 0.13 5.07 0.00 0.00
> system.time(for(i in 1:1000) loglikelihood2(0.2))
[1] 0.20 0.02 0.27 0.00 0.00

二つのプログラムが同じ答えを出すことを確認

> loglikelihood(0.2)
[1] -884.4943
> loglikelihood2(0.2)
[1] -884.4943

お見事!(まさかFAQではないですね) 元記事投稿者

(73) R オブジェクト(参照情報)を BibTeX, LaTeX 形式に変換表示する toLatex() (2006.02.17)

推奨使用例は (72) を参照

toLatex                package:utils                R Documentation

Converting R Objects to BibTeX or LaTeX

Description:

     These methods convert R objects to character vectors with BibTeX
     or LaTeX markup.

Usage:

     toBibtex(object, ...)
     toLatex(object, ...)
     ## S3 method for class 'Bibtex':
     print(x, prefix="", ...)
     ## S3 method for class 'Latex':
     print(x, prefix="", ...)

(72) セッション情報 sessionInfo() (2006.02.16)

比較的最近登場した関数(R-2.0.0以降搭載) 恐らく使用環境(R バージョン、OS, 使用パッケージ)を明らかにせず質問をする初心者に業を煮やした常連回答者が付け加えた機能。R を使った研究論文の引用情報にも使うことを想定しているよう。例えば私の場合次のようになる:

> sessionInfo()
R version 2.2.0, 2005-10-06, i486-pc-linux-gnu

attached base packages:
[1] "methods"   "stats"     "graphics"  "grDevices" "utils"     "datasets"
[7] "base"

うぃんどうず版だとこんな風になる

 > sessionInfo()
R version 2.2.1, 2005-12-20, i386-pc-mingw32

attached base packages:
[1] "methods"   "stats"     "graphics"  "grDevices" "utils"     "datasets"
[7] "base"

help(sessionInfo) を引用すると

sessionInfo              package:utils              R Documentation

Collect Information About the Current R Session 

Description:

     Print version information about R and attached packages.

Usage:

     sessionInfo(package=NULL)
     ## S3 method for class 'sessionInfo':
     print(x, ...)
     ## S3 method for class 'sessionInfo':
     toLatex(object, ...)

Arguments:

 package: a character vector naming installed packages.  By default all
          attached packages are used. 

       x: an object of class '"sessionInfo"'.

  object: an object of class '"sessionInfo"'.

     ...: currently not used.

Examples:

     sessionInfo()
     toLatex(sessionInfo())

(71) 関数行列・ベクトルの定義と使用の例 (2006.2.3)

以下は学生のために考えた工夫の簡略版です(甘い教師)。すこし複雑ですが、参考のために紹介します。要点は8個の変数を持つ6個の方程式系の解を求めることですが、元の問題の性質上、変数は二つのベクトル th, r と一つの行列 Th という構造を持っています。また、関数自身も、ベクトルと行列の構造を持っています。この構造をできるだけ活かしてコーディングしないと、手に負えなくなります(本来の問題では変数の数は数十・数百!)。解を求めるには関数の自乗和を nlm 関数で最小化(ゼロで最小になれば解のはず!)します。関数行列を実現するには、(関数)リストに次元属性を与え、行列風にアクセスしています。(もっと気の利いたやりかたを御存じなら教えてください)

Fh <- as.list(rep(NA,4))
fh <- as.list(rep(NA,2))

dim(Fh) <- c(2,2) # リスト Fh を行列(風)にする

# Fh の各成分に関数を代入(適当)
Fh[[1,1]] <- function(th, r, Th) { Th[1,1]+th[2]+r[1]-3 }
Fh[[1,2]] <- function(th, r, Th) { Th[1,2]*th[1]-r[2]   }
Fh[[2,1]] <- function(th, r, Th) { Th[2,1]-th[2]*r[1]+1 }
Fh[[2,2]] <- function(th, r, Th) { Th[2,2]*th[2]-r[2]+2 }

# 関数のベクトル(リスト)(適当)
fh[[1]]  <- function(th, r, Th) { Th[1,1]-th[1]+1 }
fh[[2]]  <- function(th, r, Th) { Th[2,2]*th[2]-1 }

# チェック
Fh[[1,1]](runif(2), runif(2), matrix(runif(4),2,2))
fh[[1]](runif(2), runif(2), matrix(runif(4),2,2))

# 最小化すべき目的関数(nlm はベクトル引数しか認めない)
F <- function(x) {
       th <- x[1:2]   # 引数ベクトル x をベクトル th, r と行列 Th に変換
       r <- x[3:4]
       Th <- matrix(x[-(1:4)], 2, 2)
       res <- 0
       for (i in 1:2) { # ループ処理ができることが行列・ベクトル風にアクセスできる最大の利点
         res <- res + fh[[i]](th, r, Th)^2
         for (j in 1:2) {
           res <- res + Fh[[i,j]](th, r, Th)^2
       }}
       return(res)    # 全部で6個の関数の自乗和
     }
  
initparam <- runif(8) # 適当な初期値
F(initparam)          # チェック

# nlm による最小化
Nlm <- nlm(F, initparam)
# 結果
> Nlm
$minimum   # 最小値(ほぼゼロ!?)
[1] 6.021328e-11

$estimate  # 推定パラメータ(つまり連立方程式の解)
[1] 1.6511474 1.1744247 2.9999900 0.4674779 0.6511482 0.3792736 1.8169112
[8] 0.8514759
 
$gradient  # これも(ほぼ)ゼロベクトルになっているべき
[1] -9.835981e-07  3.578133e-07  8.322874e-07  0.000000e+00 -1.153805e-06
[6]  1.404144e-06 -9.657654e-07 -1.673428e-07

$code  # 収束コード(おそらく解、という意味)
[1] 2

$iterations
[1] 20

(70) 引数の数を与えて空の関数を作る (Q&A 中級者コース記事より, 2006.1.31)

f <- G(3) 等とすると空の関数 funtion(x1,x2,x3) {} を作ってくれる関数。関数本体はあとで body(f) <- 何やら として定義する。
(1) テキストファイル経由の方法。定義を一時テキストファイルに書きだし、それを読み込む(野暮ったいが汎用性のある方法)

> G <- function(n) {
         tmpfile <- tempfile()                     # 一時ファイル名を作る 
         .x <- paste("x", 1:n, ",", sep="")        # 引数リスト文字列ベクトル
         .x[n] <- substr(.x[n], 1, nchar(.x[n])-1) # 最後の引数文字列からカンマを取る
         sink(tmpfile)
         cat("function (", .x, ") {}?n")           # 関数定義文字列をファイルに出力
         sink()                                    # ファイル出力を閉じる
         .f <- source(tmpfile)[[1]]                # ファイルから関数定義を読み込む
         unlink(tmpfile)                           # 一時ファイルを削除
         return(.f)
      }
  
> f5 <- G(5)
> f5
 function ( x1, x2, x3, x4, x5 ) {}
> str(f5)
 function (x1, x2, x3, x4, x5)
  - attr(*, "source")= chr "function ( x1, x2, x3, x4, x5 ) {}"

(2) もっとスマートな方法

G <- function(n) eval(parse(text=paste("function(",paste("x",1:n,sep="",collapse=","),"){}")))

(69) 行列の冪乗関数 (from r-help, 2005.8.18)

正方行列 A の任意冪乗 A^n を計算する R 関数二つ。R で単に A^3 等とすると、A*(A*A) と成分毎の冪乗とされ、A%*%(A%*%A) とはならないので注意が必要。

look at function ?mtx.exp() in the Malmig package, e.g.,
  A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)
  n.zero <- c(1,0,0)
  library(Malmig)
  all.equal(A %*% (A %*% n.zero), mtx.exp(A, 2) %*% n.zero)
  all.equal(A %*% (A %*% (A %*% n.zero)), mtx.exp(A, 3) %*% n.zero)
  mtx.exp(A, 15) %*% n.zero
Also, there is expm() in the Matrix package. I'm not sure about the
relative merits. expm() is one of the less dubious implementations of
matrix exponentials, but it does require to and from class "Matrix".
mtx.exp is a bit unfortunately named as it appears to calculate matrix
*powers* (which in this case is what you need).

注意:行列の冪乗を計算する方法にはいろいろあり、怪しげなものも多いらしいとのコメントがついていた。


「怪しげなものも多いらしい」って,どういうことでしょ。 結果がおかしいのなら,怪しいというのはあたらないでしょうし。 遅いかもしれないが,再帰関数で書いてやればいかが。もっとも n がスタックがあふれるほど大きいとかには対応していないけど。こういうものが「怪しい」のかな。。

> A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)
> fun <- function(x, n)
+ {
+ 	if (n==1) x
+ 	else x%*%fun(x, n-1)
+ }
> ans1 <- fun(A, 15)
> ans2 <- mtx.exp(A,15)
> all.equal(ans1, ans2)
[1] TRUE
> ans1
           [,1]      [,2]      [,3]
[1,] 0.46894444 1.5197281 2.7158591
[2,] 0.16295155 0.4689444 0.5859337
[3,] 0.05859337 0.2715859 0.3517577
> ans2
           [,1]      [,2]      [,3]
[1,] 0.46894444 1.5197281 2.7158591
[2,] 0.16295155 0.4689444 0.5859337
[3,] 0.05859337 0.2715859 0.3517577

あそうか。if (n==0) diag(nrow(x)) にすれば,n=0 のときにも対応できる。
n が負の整数のときには mtx.exp も対応できてないし。
なんだあ。ソース見たら,繰り返し計算で実現してるだけじゃん。

> mtx.exp
function (X, n) 
{
    if (n != round(n)) {
        n <- round(n)
        warning("rounding exponent `n' to", n)
    }
    phi <- diag(nrow = nrow(X))
    pot <- X
    while (n > 0) {
        if (n%%2) 
            phi <- phi %*% pot
        n <- n%/%2
        pot <- pot %*% pot
    }
    return(phi)
}

「関心があれば、以下をお読みください」だそうです(元記事投稿者):

Moler C., van Loan C., (2003); _Nineteen dubious ways to compute
    the exponential of a matrix,  twenty-five years later_, SIAM
    Review 45, 3-49.

便利な世の中になりました。インターネットで原論文が得られるようです。

(68) 複数変数への一括代入 (from r-help, 2005.7.9)

次の G. Grothendieck 氏による list[.,.] という関数は二つの変数への同時代入を行なう。 すこしいじれば3変数版もできそう。

list <- structure(NA,class="result")
"[<-.result" <- function(x,...,value) {  # list[.,.] 関数の定義本体
  args <- as.list(match.call())
  args <- args[-c(1:2,length(args))]
  length(value) <- length(args)
  for(i in seq(along=args))
  eval(substitute(x <- v,list(x=args[[i]],v=value[[i]])),env=sys.frame(-1))
  x
} 
 
> a <- 1; b <- 2
> list[a,b] <- list(b,a)  # 中間変数を使わないスワップ
> a; b 
[1] 2
[1] 1

> list[a,b] <- runif(2) # 変数 a,b への同時代入
> a; b
[1] 0.6006348
[1] 0.6300699

# 回帰結果から Coef, Rsid 成分を同名の変数に同時代入   
> list[Coef, Resid] <- lm(rnorm(10) ~ seq(10)) 
> Coef
(Intercept)     seq(10) 
 -0.0375528   0.1387487 
> Resid
          1           2           3           4           5           6 
-0.33104288  0.07878998 -0.62926542  0.88475733  1.15494717  0.72894189 
          7           8           9          10 
-1.69131959 -0.78797801 -0.59520853  1.18737807 
# eigen 関数の返り値から eval, evev 成分を同名の変数に同時代入 
> list[eval, evec] <- eigen(cbind(1,1:3,3:1))
> eval
[1]  5.000000e+00 -1.000000e+00 -2.868799e-16
> evec
          [,1]       [,2]       [,3]
[1,] 0.5773503 -0.8082904 -0.9428090
[2,] 0.5773503 -0.1154701  0.2357023
[3,] 0.5773503  0.5773503  0.2357023

(67) 行列関数の微分 (from r-help, 2005.07.06)

実際は関数リスト行列の成分毎の微分

> e <- lapply(1:4, function(i) bquote(x^.(i)))
> dim(e) <- c(2,2)
> De <- lapply(e, D, "x")
> dim(De) <- c(2,2)
> e
     [,1]       [,2]      
[1,] Expression Expression
[2,] Expression Expression
> De
     [,1]       [,2]      
[1,] 1          Expression
[2,] Expression Expression
> e[[1,2]]
x^3
> De[[1,2]]
3 * x^2

(66) match (2005.06.06)

match については,このページでも記述があるが,そこでは言及されていない使い道について。
factor, as.factor は,アルファベット順に整数値を割り当てるが,任意の数値の割り当てを定義するために match を使うことができる。match の第二引数で定義される順序で,第一引数の文字列を数値化できる。

> x <- c("Male", "Female", "Male", "Female", "Female")
> x
[1] "Male"   "Female" "Male"   "Female" "Female"
> as.factor(x)
[1] Male   Female Male   Female Female
Levels: Female Male
> as.integer(as.factor(x))
[1] 2 1 2 1 1
> match(x, c("Male", "Female"))
[1] 1 2 1 2 2

(65) (整数)ベクトル x に対する x[!!cumsum(!!x)] (from r-help, 2005.05.31)

整数ベクトルの先頭部分にある 0 をあるだけ取り除きたい(ただし途中にある 0 はそのままで)という質問に対する R の魔術師 G. Grothendieck 氏の回答。一読何をしているか首をかしげたというだけの理由で紹介します。以下の例を見ればなるほどとわかります。要するに !! で 0 は FALSE (=0)、ゼロ以外は TRUE (=1) に変えているわけです。

x = c(0,0,1,0,-2,0,0,3,-4,5,-6,0)
!x
[1]  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE
!!x
[1] FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE
cumsum(!!x)
[1] 0 0 1 1 2 2 2 3 4 5 6 6
!cumsum(!!x)
[1]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
!!cumsum(!!x)
[1] FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
x[!!cumsum(!!x)]
[1]  1  0 -2  0  0  3 -4  5 -6  0

なるほど,これはトリッキー過ぎる。私は

x[cumsum(abs(x))>0]

を考えた(後知恵だけど)。文字数が同じ以下の解も。

x[cumsum(x^2)>0]

次の解は,長すぎるし,先頭の0がないときにはだめだな。

x[-c(1:(which(x!=0)[1]-1))]

これで十分かと思うんですが?(トリック大会でなければ,むしろ後で目が点になるほうが恐いような)

x[cumsum(abs(x))!=0]

(64) R からメイルを送る (from r-help, 2005.05.11)

r-help に R からメイルを送る方法が議論されていました。おもな使い方は、時間のかかる作業が終ったらユーザーに通知することのようです。以下は Unix-like なシステムの例ですが、Windows でも適当なソフトを使えばできるようです(2005.05.11 の r-help アーカイブ参照)

For example, the following works on my Linux system:

  send.mail<-function(addr,subject="Mail from R",
                    text="empty text"){
    mail.cmd<-paste("mail ",
                    "-s \"",subject,"\" ",
                    addr,
                    " << EOT &\n",
                    text,"\n",
                    "EOT",
                    sep="",collapse="")
     system(mail.cmd,intern=FALSE)
  }

For example, with: 

  send.mail("ted@localhost",
  subject="Message from R: At Last!",
  text="Good Morning!\nI've been busy all night.\nIt's finished."
  )

after a minute or two for delivery I get the following message
delivered on my machine:

 From ted@compo.my.LAN  Tue May 10 09:50:27 2005
 Date: Tue, 10 May 2005 09:49:50 +0100
 From: Ted Harding <ted@compo.my.LAN>
 To: ted@compo.my.LAN
 Subject: Message from R: At Last!

 Good Morning!
 I've been busy all night.
 It's finished.

I've specified "intern=FALSE" explicitly (though it is the
default), since if you set it TRUE (e.g. in order for R to
verify the sending) then R will hang until it receives the
feedback from the command (the "&" in "EOT &\n" puts the job
in the background so there is no delay).

Hoping this helps! (Of course, if you're not a unixoid, this
is probably no use to you at all, and I have no idea how to
do anything similar in Windows).

Best wishes,
Ted.

(63) 一般関数の数値微分 (from r-help, 2005.05.08)

Richardson 補外による一般関数の数値微分 (r-help 記事引用)

Hi, 

There is a nice article by Fornberg and Sloan (Acta Numerica 1994) on
various order accuracy (Taylor-series based) approximations for different
order derivatives. I had coded a couple of these in R for first and second
order derivatives, with truncation errors of orders 2 and 4.  

Here is the code and a simple example demonstrating its accuracy for a
sharply oscillating function:

############################################
# A function to compute highly accurate first- and second-order derivatives
# From Fornberg and Sloan (Acta Numerica, 1994, p. 203-267; Table 1, page 213)
deriv <- function(x, fun, h=NULL, order=1, accur=4) {
macheps <- .Machine$double.eps

if (order==1) {
if(is.null(h)) h <- macheps^(1/3)* abs(x)
ifelse (accur==2, w <- c(-1/2,1/2), w <- c(1/12,-2/3, 2/3,-1/12))
ifelse (accur==2, xd <- x + h*c(-1,1), xd <- x + h*c(-2,-1,1,2))
return(sum(w*fun(xd))/h)
}
else if (order==2) {
if(is.null(h)) h <- macheps^(1/4)* abs(x)
ifelse (accur==2, w <- c(1,-2,1), w <- c(-1/12,4/3,-5/2,4/3,-1/12))
ifelse (accur==2, xd <- x + h*c(-1,0,1), xd <- x + h*c(-2,-1,0,1,2))
return(sum(w*fun(xd))/h^2)
}
}
############################################
func1 <- function(x){
sin(10*x) - exp(-x)
}
#############################################
> curve(func1,from=0,to=5) 

> x <- 2.04
# first order derivative 
> numd1 <- deriv(x,f=func1)
> exact <- 10*cos(10*x) + exp(-x)
> c(numd1,exact,numd1/exact-1)
[1] 3.335371e-01 3.335371e-01 1.981793e-11
# second order derivative 
> numd2 <- deriv(x,f=func1,order=2)
> exact <- -100*sin(10*x) - exp(-x)
> c(numd2,exact,numd2/exact-1)
[1] -1.001093e+02 -1.001093e+02 -2.300948e-11

Hope this is helpful.

Best, Ravi.
A current version of this code is in the dseplus bundle in the devel
section of CRAN. It does Richardson's extrapolation, which gives an
accurate numerical estimate at the expense of speed. (It does a very
large number of function evaluations.) That may not be what you want.

Paul Gilbert

(62) 行列のランクを求める (from r-help, 2005.5.5)

線形代数では中心的な概念であるランクであるが、一方数値計算的には求めるのがハード(微妙)であることは、知る人ぞ知る。r-help に紹介されていた印象的な例。 hilbert(n) はオーダー n のヒルベルト行列を計算する。これはランク n を持つが、それを数値的に計算するのが困難であることで有名。

> hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
> qr(hilbert(9), tol = 1e-8)$rank  # 一応こうすればランクは計算できる
> [1] 9                            # 正解!
> h9 <- hilbert(9)
> temp <- cbind(h9, h9)
> h9times4 <- rbind(temp, temp)    # hilbert(9) を4つ並べた 18x18 正方行列
> qr(h9times4,tol=1e-7)$rank       # 正解はやはり 9 になるはずが?!  
[1] 7
> qr(h9times4, tol=1e-8)$rank
[1] 10
> qr(h9times4, tol=1e-9)$rank
[1] 11
> qr(h9times4, tol=1e-10)$rank
[1] 12

(61) R の構文を引数として関数に渡す (Q&Aコーナー記事より、2005.4.30)

R の構文を表す文字列を(引用符付き、もしくはなしで)引数として与え、それを実行したい。以下の関数 test は構文的に正しい文字列(というべきなのか)を引数としてとり、それを実行します。

> x <- 1:10
> y <- rnorm(10)
> test <- function(.x) {eval(substitute(.x))} # 混乱を避けるため変数名を .x とした
> test(x^2)
 [1]   1   4   9  16  25  36  49  64  81 100
> test(y*x)
 [1]   0.3400069  -1.0557752   0.9980655  -9.7466238   4.0248083  -2.1227256
 [7]   1.2294612  -0.6428637 -14.4409546  12.6976043
> test(lm(y~x))

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x

同じことを構文を表す(引用符つき)文字列を引数として表すと

> test2 <- function(.x) {eval(parse(text=substitute(.x)))}
> test2("x^2")
 [1]   1   4   9  16  25  36  49  64  81 100
> test2("x*y")
 [1]   0.3400069  -1.0557752   0.9980655  -9.7466238   4.0248083  -2.1227256
 [7]   1.2294612  -0.6428637 -14.4409546  12.6976043
> test2("lm(y~x)")

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x
   -0.38949      0.03299

(60) 行列の一般化内積 (from r-help, 2005.4.30)

r-help の wizard Gabor Grothendieck 氏による次の関数は [i,j] 成分が f(a[i,],b[,j],...) である行列を生成する。内側の apply 関数は i 成分が f(a[i,],x,...) であるベクトルを作り、外側の apply 関数はそれを各 x=b[,j] について繰返し、行列とする。 f <- match.fun(f) は inner の関数引数 f を確実に特定するために使われている。

inner <- function(a,b,f, ...) {
           f <- match.fun(f)
           apply(b,2,function(x)apply(a,1,function(y)f(x,y, ...)))
}

# 例:spearman 相関行列
data(iris); irish <- head(iris)[,-5]  # test data
res <- inner(t(irish), irish, cor, method = "spearman")  # test
res2 <- cor(irish, method = "spearman") # should give same result
identical(res, res2) # TRUE

(59) 文字列ベクトルからその簡略形を作る、関数 abbreviate(), 2005.4.30

使用例はすぐ下の記事を参照。

Usage:

 abbreviate(names.arg, minlength = 4, use.classes = TRUE, dot = FALSE)

Arguments:

names.arg: 省略すべき文字列ベクトル
minlength: 省略文字列の最小長
use.classes: 論理値 (現在 R では無視).
     dot: 論理値;ドットを加えるか

(58) データフレームの変数名を簡略化する (from r-help, 2005.04.29)

データフレームの変数(列)名はしばしば長過ぎて(それはそれで意味がわかりやすいからよいことであるが)後々操作に困る。names 関数を用いて扱いやすい名前に変えることは簡単だが、たくさんあると面倒。以下はそれを自動化する二つの工夫。

> data(iris)   # 例、組み込みデータ iris
> iris
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1            5.1         3.5          1.4         0.2     setosa
2            4.9         3.0          1.4         0.2     setosa
................................................................
> data(iris)
> names(iris) <- gsub("[a-z.]", "", names(iris)) # 列名から小文字とドットを取り去る(重複する可能性あり)
> names(iris)[5] <- "Class"                      # 種名だけは手作業で変更
> 
> iris
     SL  SW  PL  PW      Class
1   5.1 3.5 1.4 0.2     setosa
2   4.9 3.0 1.4 0.2     setosa
..............................
> data(iris)
> names(iris) <- abbreviate(names(iris))  # 文字列からその省略形を自動的に作る(重複しないことが保証される)
> names(iris)[5] <- "Class"
> 
> iris
    Sp.L Sp.W Pt.L Pt.W      Class
1    5.1  3.5  1.4  0.2     setosa
2    4.9  3.0  1.4  0.2     setosa
...................................

(57) 文字列ベクトルから一つを対話的に選ぶ、select.list() 2005.4.27

select.list() 関数は引数の文字列ベクトルの要素をリストアップし、そのうちの一つを対話的に選択する。類似関数に menu() 関数があるがこちらは要素番号を返し、switch(menu(...)) 等と条件付き分岐に使う。

> sort(.packages(all.available = TRUE))
 [1] "base"      "datasets"  "grDevices" "graphics"  "grid"      "methods"
 [7] "mvtnorm"   "splines"   "stats"     "stats4"    "tcltk"     "tools"
[13] "utils"
> x=select.list(sort(.packages(all.available = TRUE)))

 1: base        2: datasets    3: grDevices   4: graphics    5: grid
 6: methods     7: mvtnorm     8: splines     9: stats      10: stats4
11: tcltk      12: tools      13: utils
 
選択:7
> x
[1] "mvtnorm"   #選択した要素文字列を返す
> x=menu(sort(.packages(all.available = TRUE)))

 1: base        2: datasets    3: grDevices   4: graphics    5: grid
 6: methods     7: mvtnorm     8: splines     9: stats      10: stats4
11: tcltk      12: tools      13: utils 

選択:7
> x
[1] 7            # 要素番号を返す

(56) リストはベクトル (r-help, 2005.04.17)

r-help の R 教授のコメントによれば、R のリストはベクトルとして実装されているらしい。昔は pairlist という形式だったらしい。

> a <- list(1,"abc",1:5)
> is.vector(a)             # a はベクトル
[1] TRUE
> is.list(a)               # 当然 a はリストでもある
[1] TRUE
> a[2]                     # ベクトルとしてアクセス(ただし結果はリスト) 
[[1]]
[1] "abc"
> a[[2]]                   # リストとしてアクセス(結果はベクトル)
[1] "abc"
> str(a)
List of 3
 $ : num 1
 $ : chr "abc"
 $ : int [1:5] 1 2 3 4 5

(55) 組み込みデータへの付値に関する注意 (from r-help, 2005.04.16)

組み込みデータ(フレームとその成分変数)を参照するには attach(...) とする。しかし、うっかりその成分変数に付値を行なうと、厄介なことが起きるので注意が必要。

> rm(list=ls())     # オブジェクトをすべて(.GlobalEnv から)消去
> ls()              # 確かに何もない
character(0)
> attach(women)     # 組み込みデータフレーム women を参照できるようにする
> ls()              # しかし .GlobaEnv には変化なし
character(0)
> find(height)      # height 変数を探す
[1] "women"         # women というデータフレーム中にあると回答 
> height[1]
[1] 58
> height[1] <- 48   # height 変数を変更
> ls()
[1] "height"        # height 変数が作られている
> find(height)
[1] ".GlobalEnv" "women"  # height 変数は .GlobalEnv と women の双方にあると回答
> women$height      # women 中の height 変数を参照する確実な方法
 [1] 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
> height            # 一方この height 変数は新たに作られたオブジェクトを指している
 [1] 48 59 60 61 62 63 64 65 66 67 68 69 70 71 72
> detach(women)     # women を取り外す
> ls()              # height 変数は(当然)残る
[1] "height"

オンラインヘルプに書いてある。 Note that by default assignment is not performed in an attached database. Attempting to modify a variable or function in an attached database will actually create a modified version in the user's workspace (the R global environment).

上の匿名氏の(そして r-help の R 教授の)解説のように、attach(women) はデータフレーム women の成分を展開した環境を作り、それを検索パス中に登録するという動作を行なうらしい。次の R 教授による例は更に微妙で、attach、detach 関数が環境を操作する関数であることが良くわかる。

> attach(women)
> height            # この height は women 展開環境中にコピーされた women$height           変数
 [1] 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
> height[1] <<- 0   # この永続付値は women を展開した環境中の height オブジェクトを操作している
> height
 [1]  0 59 60 61 62 63 64 65 66 67 68 69 70 71 72
> women$height      # これは women の成分を直接参照している 
 [1] 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
> height[1] <- 10   # この付値は .GlobalEnv 中に height のコピーを作る
> height            # 単に height とすれば .GlobalEnv 中の height とされる 
 [1] 10 59 60 61 62 63 64 65 66 67 68 69 70 71 72
> rm(height)        # .GlovalEnv 中の height を消す
> height            # 今度は women 展開環境中の height のこととされる
 [1]  0 59 60 61 62 63 64 65 66 67 68 69 70 71 72
> detach(women)     # women 展開環境を消す
> height            # その中の height オブジェクトも道連れで消える 
Error: Object "height" not found
# いずれにしても women データフレーム自身は何ら変更されることはない

(54) NA, NaN, Inf 等がある場合のベクトルの表集計 (from r-help, 2005.04.09)

ベクトル中に NA, NaN, Inf 等がある場合の表集計あれこれ

> x <- c(1,2,3,3,NA, NaN, Inf)
> table(x)       
x
  1   2   3 Inf                          # NaN, NA は除外
  1   1   2   1  
> y <- factor(x)
> table(y)
y
  1   2   3 Inf NaN                 # NA だけ除外
  1   1   2   1   1  
> summary(y)                    #  何でも一応 summary 関数を使ってみると得するかも知れないという例 
   1    2    3  Inf  NaN NA's   # NA, NaN, Inf すべて集計
   1    1    2    1    1    1 

(53) 関数中で未定義のオブジェクトを未定義エラーとする、関数 environment, new.env (fro r-help, 2005.04.09)

関数中で未定義のオブジェクトは既定では自動的にその親(更には先祖)環境の名前空間中で探され、もしあればその値が使われる。しかし、場合によれば、うっかりミスを防ぐため、そうして欲しくない(未定義エラーとなって欲しい)こともある。

> a = 2.0
> mult <- function(x)   return(pi*a*x) 
> mult(1)
[1] 6.283185                       # 関数 mult は変数 a をその親環境の a とした
> environment(mult) <- new.env(parent = NULL)  # 関数 mult の親環境は無いと宣言
> mult(1)            
Error in mult(1) : Object "a" not found  # 親環境がないから変数 a は未定義エラーとされる

(52) リストの成分名を変数で与える (Q&A 初級者記事より、2005.04.07)

> hoge <- function(name1="a", name2="b") 
           {  tmp <- list(rnorm(2), runif(2))
              names(tmp) <- c(name1, name2)
              tmp
           } 
> hoge()
$a
[1] -1.627812797 -0.005876096
$b
[1] 0.8406733 0.6035658
> hoge("A","B")
$A
[1] 1.111073 1.008163
$B
[1] 0.3292604 0.1637242
#次では成分名は "name1, "name2" のままで変わらない
> hoge2 <- function(name1="a", name2="b") 
                        list(name1=rnorm(2), name2=runif(2))

(51) バッチモード実行の R エラーメッセージを標準出力にリダイレクト (from r-help, 2005.04.07)

未確認。linux 以外では? prova が実行スクリプト。

$ R -q --no-save < prova > prova.out 2>> prova.out
$ R -q --no-save < prova > prova.out 2>&1                     # 標準エラー出力

(50) ベクトル x に対する as.vector(x) (2005.04.07)

ベクトル x に対し as.vector(x) とするのは、無意味に見えるが、x に names 属性がついているとそれを除くという重要な副作用がある。 table 関数の結果を頻度ベクトルにするのに使える。

> (x <- sample(1:5, 10, replace=TRUE))
 [1] 4 4 5 5 3 4 2 2 4 5
> names(x) <- letters[1:10]  # 名前ラベルをつける
> x
a b c d e f g h i j 
4 4 5 5 3 4 2 2 4 5 
 > is.vector(x)
[1] TRUE
> as.vector(x)                        # 名前ラベルが消える
 [1] 4 4 5 5 3 4 2 2 4 5
> x <- sample(1:5, 10, replace=TRUE)
> (y <- table(x))    # 集計表
x
1 2 3 5 
2 3 2 3 
> str(y)
 int [, 1:4] 2 3 2 3
 - attr(*, "dimnames")=List of 1   # dimnames 属性というリストがついている
  ..$ x: chr [1:4] "1" "2" "3" "5"
 - attr(*, "class")= chr "table"       # y はクラス ”table" のオブジェクト
> is.matrix(y)                                 # 行列でもないし
[1] FALSE
 > is.vector(y)                                # ベクトルでもない
[1] FALSE
> is.table(y)                                   # table は table 
[1] TRUE
> y[1]+y[2]                                    # ただしベクトル演算も可能
1 
5 
> sum(y)
[1] 10
> as.vector(y)
[1] 2 3 2 3                                     # 頻度ベクトル(これだけが欲しいことも良くある)
> dimnames(y) <- NULL   <- これでも良い
> y
[1] 2 3 2 3

(49) do.call(func, list) は(名前つき成分)リストを引数として、名前が func という関数を実行する。(2005.04.06)

 
# リストの行列化
> L <- list(a=rnorm(5), b=runif(5)) 
> do.call("rbind", L)     # 行列化
        [,1]      [,2]       [,3]       [,4]      [,5]
a -0.2610772 0.5175920 -0.7930714 -0.5240486 -1.527911
b  0.1986414 0.7121043  0.8124055  0.3421709  0.738326
> do.call("cbind", L)   # 行列化
        a                   b
[1,] -0.2610772 0.1986414
[2,]  0.5175920 0.7121043
[3,] -0.7930714 0.8124055
[4,] -0.5240486 0.3421709
[5,] -1.5279110 0.7383260
ためしにリストを行列に強制変換すると、支離滅裂
> as.matrix(L)
     [,1]
[1,] Numeric,5    
[2,] Numeric,5
# リスト成分に名前が無いと...
> L <- list(rnorm(5), runif(5)) 
> do.call("rbind", L)
          [,1]        [,2]       [,3]      [,4]       [,5]
[1,] 0.6996488 1.559108982 -0.8200150 0.1542949 -0.4964565
[2,] 0.1321730 0.007599077  0.4112039 0.7412868  0.4271347

(48) R プロセスの実行優先度を設定する。パッケージ nice 。(2005.04.05)

Unix の nice 命令(Unix 以外では使える?)。バックグラウンドで R を実行し、その間に他の作業を快適に行ないたい時に有効?

nice                  package:nice                  R Documentation

Get or Set UNIX Priority

Description:

     Get or Set UNIX Priority (Niceness) of this R process.

Usage:

     get.my.priority()
     set.my.priority(priority = 15)    # 19 が優先度最低

Arguments:

priority: The UNIX priority, also called niceness, what the UNIX
          commands 'nice' and 'renice' set, (an integer number, usually
          in the range -20 to 19, that we want to set the priority to).
           Generally, ordinary users can only raise the priority, so it
          is an error to try to set the priority lower than that
          returned by 'get.my.priority'.

(47) 関数本体で、関数自体を表す関数 Recall (from r-help, 2005.03.31)

何やら意味不明の説明だが、次の ?Recall の実例を見ると意味がわかる。再帰的関数以外に使い道は?

    ## A trivial (but inefficient!) example:
    fib <- function(n) if(n<=2) {if(n>=0) 1 else 0} else Recall(n-1) + Recall(n-2)
    fibonacci <- fib; rm(fib)
    ## renaming wouldn't work without Recall
    fibonacci(10) # 55

(46) モデル式を関数引数にする as.formula() 関数 (r-help, 2005.03.27)

get.model.fit <- function(my.form)
{
    res <- lm(as.formula(my.form))
    summary(res)$r.squared
    # other stuff happens here...
} 

Then call it with:

get.model.fit("y ~ x1 + x2 + x3")

Internally, the vector will be converted to:

> as.formula("y ~ x1 + x2 + x3")
y ~ x1 + x2 + x3

Marc Schwartz

(45) オブジェクトの強制変換方法を指定する(?)関数 as -- (2005.02.21)

help(as) を読んでも意味、使い方が良くわからん。as.matrix 等の強制変換メソッドを独自に指定するための関数らしいが? とりあえず末端ユーザには不要な関数らしい。methods ライブラリには R 言語の真髄である関数がいっぱいあるが、いずれも意味がなかなか掴めない。R 言語そのものを拡張しようという方には不可欠な内容らしいのでじっくり吟味下さい。

# help(as) 冒頭の説明
These functions manage the relations that allow coercing an object
to a given class.

(44) データとその解析関数の一括リスト化 (Q&Aコーナーの質問から思いついた工夫、2005.02.19)

リストの成分には任意の R オブジェクトが許される。実際、R の関数が返す返り値リストには、結果そのものではなく、結果出力用の関数が成分として入れてある場合があるようだ。ということは、ユーザも一連の解析に登場するデータとその解析関数(更には解析結果)を、すべて一つのリストにまとめておくことができるはず。記憶力と整理整頓能力に自信のないかたは試してみる価値があるかも?私はこれまで

関数中に対象データを含めておく
関連関数とデータを一つの大きなダミー関数中に放り込んでおく
解析毎に作業ディレクトリを変える

などというやぼったい工夫をすることがありました。同様な、もしくはもっと気の効いたやりかたをご存知なら是非教えて下さい。

# つまらない例
> x <- 1:10        # データ
> y <- runif(10)   # データ
> ex1 <- function () lm(ex$y~ex$x) # その解析関数
> ex2 <- function () summary(lm(ex$y~ex$x)) # 同じく
> ex <- list(x=x, y=y, ex1 = ex1, ex2 = ex2) # 一つのリストに一括保存
> ex$ex1()
Call:
lm(formula = ex$y ~ ex$x)
Coefficients:
(Intercept)         ex$x
    0.54553     -0.00734
> ex$ex2()
Call:
lm(formula = ex$y ~ ex$x)
Residuals:
     Min       1Q   Median       3Q      Max
-0.50573 -0.34908  0.09843  0.29182  0.40047
Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.545534   0.253466   2.152   0.0635 .
ex$x        -0.007339   0.040850  -0.180   0.8619
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 0.371 on 8 degrees of freedom
Multiple R-Squared: 0.004019,   Adjusted R-squared: -0.1205
F-statistic: 0.03228 on 1 and 8 DF,  p-value: 0.8619

(43) データフレームを転置すると行列に変わる, from r-help -- 2005.02.11

おまけにデータフレームに sample 関数を適用するとケース毎に抽出が行われる(ヘルプドキュメントにも記載されていない隠し機能?!)

> mat1 <- data.frame(matrix(rnorm(10*6), ncol=10))
> mat2 <- t(mat1)
> is.data.frame(mat1)
[1] TRUE
> is.data.frame(mat2)
[1] FALSE
> is.matrix(mat2)       # 転置(そもそもデータフレームの転置など聞いたことがない)すると行列になっている!
[1] TRUE
>  dim(sample(mat1, 5, replace=T)) # ケース毎の標本抽出からなる行列が得られる
[1] 6 5
# 一方行列である mat2 に対する標本抽出は行列をベクトル化した上での操作だから、結果は長さ 5 のベクトル
>  dim(sample(mat2, 5, replace=T))
NULL
> length(sample(mat2, 5, replace=T))
[1] 5
>  dim(sample(as.data.frame(mat2), 5, replace=T))
[1] 10  5
# もう一つの例
> x = data.frame(V1=c(1,2), V2=c("a","b"))
> x
  V1 V2
1  1  a
2  2  b
> t(x)       # 文字列行列に変わっている!
   1   2  
V1 "1" "2"
V2 "a" "b"
> xx <- data.matrix(x) # データフレームを行列に変える「安全な」方法
> xx
  V1 V2  # 但し V2 列は因子とされ、その内部コーディングである数値に置き換わる
1  1  1
2  2  2
# 更に(少し奇妙な感じだが)。一体どういう場合に使うことを想定しているのか?それともたまたまの単なる副作用か?
> x <- data.frame(foo=1:5, bar=letters[1:5])
> x
  foo bar
1   1   a
2   2   b
3   3   c
4   4   d
5   5   e
> (xx <- sample(x, 6, replace=TRUE) )
  bar foo bar.1 foo.1 foo.2 foo.3
1   a   1     a     1     1     1
2   b   2     b     2     2     2
3   c   3     c     3     3     3
4   d   4     d     4     4     4
5   e   5     e     5     5     5
> is.data.frame(xx)  # 結果はデータフレームであることを注意
[1] TRUE

(コメント) Q モードの分析をしたいというときには有りそうかも。ただ,そのような場合には文字データを含むデータ列があるなんてことは,考えにくい。うっかり,やってしまうと言うことはあるかも。

(コメント) as.matrix(x) でも,数値も文字列になる。
なぜなら,行列は列ごとに別のデータタイプを持てないから。
原理的に考えれば,理の当然。(結果論だろうって,そうですよ)。

> as.matrix(x)
  V1  V2 
1 "1" "a"
2 "2" "b"

本質は,データフレームをマトリックスにしたり, データフレームを転置して(マトリックスにならざるを得ないとき), 結果を表すには一番上位のデータタイプを適用して, それよりも下位のデータタイプのデータを変更せざるを得ないと言うこと。
たとえば,数値(real)と複素数(complex)が共存すると,複素数に合わせるしかないということ。(知ってる人は知っている,知らない人は知らない。これ,あたりまえ。このページのコンテンツはそのようなもの)

> x = data.frame(V1=c(1,2), V2=c(complex(1,1,2), complex(1,3,4)))
> x
  V1   V2
1  1 1+2i
2  2 3+4i
> y <- as.matrix(x)
> y
    V1   V2
1 1+0i 1+2i
2 2+0i 3+4i

別の例を示そう。

>  x = data.frame(V1=c(1,2), V2=c(3,4), V3=c("a","b"))
> apply(x,2,mean)
V1 V2 V3 
NA NA NA 
Warning messages: 
1: argument is not numeric or logical: returning NA in: mean.default(newX[, i], ...) 
2: argument is not numeric or logical: returning NA in: mean.default(newX[, i], ...) 
3: argument is not numeric or logical: returning NA in: mean.default(newX[, i], ...) 

データフレーム全部,各列の最上位のデータタイプとして文字列を仮定されたから,平均値は求まらない。
数値の入っている列だけを指定するとちゃんと計算できる。

> apply(x[,1:2],2,mean)
  V1  V2 
 1.5 3.5 

もっとも,この場合(平均値が必用)に限れば

> summary(x)
       V1             V2       V3   
 Min.   :1.00   Min.   :3.00   a:1  
1st Qu.:1.25   1st Qu.:3.25   b:1  
 Median :1.50   Median :3.50        
 Mean   :1.50   Mean   :3.50        
 3rd Qu.:1.75   3rd Qu.:3.75        
 Max.   :2.00   Max.   :4.00        

とすればいいだけ。


(42) sum 関数の落とし穴 -- 2005.02.11

sum 関数は引数に NA, NaN があると結果はやはり NA,NaN になる。na.rm=TRUE オプションを使えばあらかじめ NA, NaN を取り除いて和を計算してくれる。しかし、空のベクトルの和は 0 を返す(数学ではおなじみの慣習)という落とし穴がある。

> sum(c(1,2,NA))
[1] NA
> sum(c(1,2,NA), na.rm=TRUE)
[1] 3
> sum(c(NA,NA), na.rm=TRUE)
[1] 0
> sum(numeric(0))
[1] 0
> sum(c(NaN))
[1] NaN
> sum(c(1,2,NaN))
[1] NaN
> sum(c(1,2,NaN), na.rm=TRUE)
[1] 3

プログラム中ではうっかり引っかかる恐れもある。面倒臭いが例えば次のような「安全な和関数」を使うと良い(?)

> safesum <- function(x, na.rm=FALSE) {
 ifelse(na.rm==TRUE, y <- na.omit(x), y <- x)
 if(length(y)==0) stop("summing an empty object!")
 else return(sum(y))}
> safesum(c(1,NA), na.rm=TRUE)
[1] 1
> safesum(c(NA,NA), na.rm=TRUE)
Error in safesum(c(NA, NA), na.rm = TRUE) :
       summing an empty object!
> safesum(numeric(0))
Error in safesum(numeric(0)) : summing an empty object!

sum( ) が0になるということが問題になるのは,統計学の場面では少ないかも。というのも,統計学で sum( ) が0になるということを実質的に把握するのは mean( ) がいくつになるかという場面ではないかな?

> x <- c(NA,NA)
> mean(x)
[1] NA
> x <- x[is.na(x)==FALSE]
> mean(x)
[1] NaN
 > x <- c()
> mean(x)
[1] NA
Warning message: 
argument is not numeric or logical: returning NA in: mean.default(x) 

mean をとるときは,有効データ数 N がいくつかを必ず確認しますから(確認しないなら馬鹿),問題は生じないことを期待したものだ(期待は常に裏切られが)。


(41) ユーザーのファイルシステムにあるファイルの情報を与える関数、file.info() -- 2005.02.06

得られる情報は OS 依存である。

> file.info("report.pdf")
            size isdir mode               mtime               ctime
report.pdf 77680 FALSE  644 2005-02-01 20:18:37 2005-02-01 20:18:37
                         atime  uid  gid uname grname
report.pdf 2005-02-01 20:19:42 1000 1000  hoge   hoge
> file.info("/home/hoge")  # ディレクトリ情報
           size isdir mode               mtime               ctime
/home/hoge 4096  TRUE  777 2005-02-06 22:07:24 2005-02-06 22:07:24
                         atime  uid  gid uname grname
/home/hoge 2005-02-06 21:13:43 1000 1000  hoge   hoge

(40) 文字列名で与えた関数をリスト引数の各成分に適用、do.call -- 2005.02.04

> ( mylist <- list(1:3, 4:6, 7:9) )
[[1]]
[1] 1 2 3
[[2]]
[1] 4 5 6
[[3]]
[1] 7 8 9
> rbind(mylist)  # リスト成分からなる行列を作りたかったが...
       [,1]      [,2]      [,3]
mylist Integer,3 Integer,3 Integer,3
> do.call('rbind', mylist)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9

(39) 多倍長整数を操作する関数群、パッケージ gmp -- 2005.02.03

パッケージ gmp は C ライブラリ cmp の R への移植である。はて、どういう使い方があるのやら?操作速度は未検証。

*** gmp/INDEX ***
add.bigz                Basic arithmetic operators for large integers
bigz                    Large sized integer values
==.bigz                 Relational Operators
[.bigz                  Extract or Replace Parts of an Object
gcd.bigz                Great common divisor, Least common multiple
gcdex                   Compute Bezoult coefficient
isprime                 Determine if number is prime
lucnum                  Compute Fibonacci and Lucas numbers
modulus                 Modulus
nextprime               Next prime number
sizeinbase              Compute size of a biginteger in a base
urand.bigz              Generate a random number
>pow.bigz(10,10000) #10の1万乗
[1] "10000・・・(0が全部で10001個)・・・00"

>#同じことを、計算結果を多倍長にする場合、
>#as.bigzを使うと、引数もas.bigzにする必要がある(1つでいい)
>as.bigz(as.bigz(10)^as.bigz(10000))
[1] "10000・・・(0が全部で10001個)・・・00"

>#引数を多倍長にするのを忘れるとInfの多倍長の結果と同じになる
>#(Do not Run!!)
>as.bigz(10^10000)
[1] "17976931348623159・・・(後略)
>as.bigz(Inf)
[1] "17976931348623159・・・(後略)

>#1000の階乗(1000!)
>kaijo<-as.bigz(1)
>for(i in 1:1000)
+	kaijo<-kaijo*i
>kaijo
[1] "40238726007709377・・・(後略)

(38) 隠しオブジェクトの表示 ls(all=TRUE) -- 2005.02.02

現在のワーキングスペースに存在するオブジェクト中、ドットで始まるオブジェクトは ls() では表示されない。場合によると巨大な隠しオブジェクトが知らない間に作られることもあるらしい。ついでにメモすると、r-help 記事によれば plot recording history 機能が有効(どうすると有効、無効にできるか不明)になっていると .SavedPlots ファイルが作られ、更に終了時に .Rdata ファイルの一部として付け加えられるため、首をかしげるほどのサイズになることがあるらしい。

> ls()
[1] "inner" "x"    
> ls(all=TRUE)
[1] ".Random.seed" ".Traceback"   "inner"        "x"    
> rm(list=ls())  # ワーキングスペースに存在するオブジェクトをすべて抹消
> ls()
character(0)
> ls(all=TRUE)
[1] ".Random.seed" ".Traceback"   # しかし隠しオブジェクトは依然として残る
> rm(list=ls(all=TRUE))           # 本当にすべてのオブジェクトを消す(賢いかどうかは別)
> ls(all=TRUE)
character(0)

(37) 引数オブジェクトに指定の属性を与えたオブジェクトを返す関数、structure -- 2005.02.01

> structure(1:6, dim = 2:3)  # ベクトル 1:6 に次元属性 dim を与えたオブジェクト(つまり 2x3 行列)
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> x <- structure(1104508800, class = c("POSIXt","POSIXct")) # 時間クラス属性を与える
> x
[1] "2005-01-01 01:00:00 JST"
> class(x) <- NULL           # class 属性をとれば元に戻る(つまり文字列に変わったわけではない)
> x 
[1] 1104508800

(36) オブジェクト名にマッチする関数を返す関数、match.fun -- 2005.01.30

match.fun 関数は与えられたオブジェクト名を持つ関数を探し、その値を返す。関数名と同じ名前のオブジェクトを「うっかり」作ってしまった場合にも、関数だけを返すので安全になる。プログラム中で変数に関数オブジェクトを付値する安全な方法。

> f <- sin        # 変数 f に関数 sin を付値
> f(0.1)          # 当然 sin(0.1) を計算する
[1] 0.09983342
> sin <- 1:10     # うっかり sin という名前のベクトルを作ってしまった
> f <- sin
> f(0.1)           # いまや sin はベクトルであるから当然エラーとなる
Error: couldn't find function "f"
> f <- match.fun(sin) # match.fun 関数は sin という名前の関数オブジェクトとマッチする
> f(0.1)           # 正しく sin(0.1) を計算
[1] 0.09983342
# 使い方の例 (r-help 記事より)
inner <- function(a,b,f) {
	f <- match.fun(f)    # 引数 f にマッチする関数オブジェクトを変数 f に付値
	apply(b,2,function(x)apply(a,1,function(y)f(x,y)))
}
inner(matrix(1:4,2), matrix(4:1,2), crossprod)

(35) (インストール済みの)パッケージの Description ファイルを見る、packageDescription 関数 -- 2005.01.27

> packageDescription("fields")
Package: fields
Version: 1.7.2
Date: Novemeber 12 2004
Title: Tools for spatial data
Author: Doug Nychka
Maintainer: Doug Nychka <nychka@ucar.edu>
Description: Fields is a collection of programs for curve and function
        fitting with an emphasis on spatial data and spatial
        statistics. The major methods implemented include cubic,
        robust, and thin plate splines, universal Kriging and Kriging
        for large data sets. One main feature is any covariance
        function implemented in R can be used for spatial prediction.
        There are also some useful functions for plotting and working
        with spatial data as images. This package also contains
        implementation of a smooth wavelet basis suitable for spatial
        models.
License: GPL Version 2 or later.
URL: http://www.cgd.ucar.edu/stats/Software/Fields
Packaged: Sat Nov 20 12:32:00 2004; nychka
Built: R 2.0.1; i386-pc-linux-gnu; 2005-01-19 19:22:15; unix

-- File: /usr/local/lib/R/site-library/fields/DESCRIPTION 

(34) 現在の作業環境をファイルにセーブし、それを次回に復元する方法。save.image 関数 -- 2005.01.27

R を q() 関数で終了する際、yes と答えると現在の作業環境(ie. 現在のワークスペース中のすべてのオブジェクト)は起動ディレクトリの XDR というバイナリ書式ファイル .Rdata に保存され、次回 R を起動するとそれが自動的に読み込まれる。もし作業途中の環境を独自に保存したければ save.image 関数を使う。

save.image(file = ".RData", version = NULL, ascii = FALSE,
           compress = FALSE, safe = TRUE)

オプションでアスキーセーブ、圧縮ファイルにすることもできる。 'save.image()' は単に 'save(list = ls(all=TRUE), file = ".RData")' と同値で、'q("yes")' としたのと同じ(但し R を終了しない)となる。もし作業途中の環境を独自に保存したければ、例えば

save.image(file=".Rdata.20050127") 

等とする。次回の R セッションをちょうどこの環境で始めたければ、R をオプション付き(.Rdata を自動読み込みしない)で

R --no-restore

と開始し(詳しくは help(Startup) 参照)、次に命令

load(".Rdata.20050127")

を実行する。大量のファイルを読み込む作業を毎回繰り返したくない時や、試行錯誤で多くの大事なオブジェクトを作った時に、便利かも。


(33) 因子毎の平均を求める ave 関数、2005.01.26

ave(x, ..., FUN=mean) は引数 x を(任意個の)第二引数(通常因子)毎にグループ化し、その平均(既定、任意関数で良い)を求める。つまり x <- ave(y, g1, g2, FUN) に対し、x[i] は

FUN( {y[j]; g1[j]=g1[i], g2[j]=g2[i]} ) 

となる。

> attach(warpbreaks)
> breaks
 [1] 26 30 54 25 70 52 51 26 67 18 21 29 17 12 18 35 30 36 36 21 24 18 10 43 28
[26] 15 26 27 14 29 19 29 31 41 20 44 42 26 19 16 39 28 21 39 29 20 21 24 17 13
[51] 15 15 16 28
> 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
> (x <- ave(breaks, wool))
 [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

# ここで x[1] は wool[1] と同じ因子(つまり A)を持つ breaks 変数のサブセットの平均値、つまり
> mean(breaks[wool==wool[1]])
[1] 31.03704
# 打ち切り平均を使った例
> ave(breaks, wool, FUN = function(x) mean(x, trim = 0.1))
 [1] 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174
 [9] 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174
[17] 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174 29.52174
[25] 29.52174 29.52174 29.52174 24.73913 24.73913 24.73913 24.73913 24.73913
[33] 24.73913 24.73913 24.73913 24.73913 24.73913 24.73913 24.73913 24.73913
[41] 24.73913 24.73913 24.73913 24.73913 24.73913 24.73913 24.73913 24.73913
[49] 24.73913 24.73913 24.73913 24.73913 24.73913 24.73913
# 各 breaks[i] と同じ wool 因子を持つ breaks 変数の個数
> ave(breaks, wool, FUN=length)
 [1] 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27
[26] 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27
[51] 27 27 27 27

(32) プロット点にラベルをつける identify 関数、2005.01.24

散布図のある点をマウスの第一ボタンでクリックすると対応する z 中のラベルをその位置に表示する。一度にいくつでもラベルをつけることができる。終了するには画面の任意位置をマウスの第一ボタン以外のボタンでクリックする。散布図中の特別な点を識別するのに便利。

# データ番号をラベルにつける
# z 引数を省略すると z=seq(along=x) とされる、つまり点の番号自体がラベルになる
x <- runif(26)
y <- runif(26)
z <- 1:26       # 点の番号に対応するラベルをつける(つまり点 (x[1],y[1]) にはラベル z[1])
plot(x,y)
identify(x,y,z) # 終了すると選択された点の番号が返される
[1]  1  2  3  4  5  6  8  9 10 11 12 13 14 16 17 18 19 20 21 22 23 24 25 26
# ラベルは文字列でも良い
x <- runif(26)
y <- runif(26)
z <- letters      # a,b,c,...
plot(x,y)
identify(x,y,z)

(31) 多変量データの thin plate spline 近似 -- 2005.01.23

fields パッケージの Tps 関数は多変量データのスプライン関数近似を計算する。次元は任意で、同じ位置に複数データがあっても良い。より詳しい解説

# Wendelberger's function
> g <- function(x, y) {
         .75*exp(-((9*x-2)^2 + (9*y-2)^2)/4) +
         .75*exp(-((9*x+1)^2/49 + (9*y+1)^2/10)) +
         .50*exp(-((9*x-7)^2 + (9*y-3)^2)/4) -
         .20*exp(-((9*x-4)^2 + (9*y-7)^2)) }
# テスト用三次元関数
> f <- function(x, y, z){1.0 + 0.5*z + 0.2*z*g(x,y)}
# テスト用データ発生領域 [0,200]x[0,100]x[0,10]
> X <- expand.grid( x=c(50, 100, 150), y=c(25, 75), z=c(2, 8))
> X <- as.matrix(X)
> X <- rbind(X, X[1,])   # わざと同じ座標値を入れる

> library(fields)      # パッケージ fields ロード

Y <- numeric(0) 
# 誤差無しの目的変数でチェック
> for ( i in 1:dim(X)[1]) Y[i] <- f(X[i,1],X[i,2],X[i,3]) + rnorm(1, mean=0, sd=0)
# thin plate spline 関数でデータ (X,Y) を近似
> tpsfit <- Tps(X, Y, scale.type="unscaled")
> YY <- predict(tpsfit, X)
> mean(abs(Y-YY))          # 平均絶対誤差
[1] 2.220446e-16           # ほぼデータ値を復元

> Y <- numeric(0)
> for ( i in 1:dim(X)[1]) Y[i] <- f(X[i,1],X[i,2],X[i,3]) + rnorm(1, mean=0, sd=0.5)
> # thin plate spline 関数でデータ (X,Y) を近似
> tpsfit <- Tps(X, Y, scale.type="unscaled")
> YY <- predict(tpsfit, X)
> mean(abs(Y-YY))          #  平均絶対誤差
[1] 0.2446144

> Y <- numeric(0)
> for ( i in 1:dim(X)[1]) Y[i] <- f(X[i,1],X[i,2],X[i,3]) + rnorm(1, mean=0, sd=1)
> # thin plate spline 関数でデータ (X,Y) を近似
> tpsfit <- Tps(X, Y, scale.type="unscaled")
> YY <- predict(tpsfit, X)
> mean(abs(Y-YY))          #  平均絶対誤差
[1] 0.589194

# Tps による近似関数を定義
> F <- function(x,y,z) predict(tpsfit, matrix(c(x,y,z), nr=1))
> F(100,30,25)     
[1] 11.06612

(30) データフレームの既存の欄を変形、新しい欄を加える、transform() -- 2005.01.20

> airquality   # 組み込みデータフレーム airquality 
    Ozone Solar.R Wind Temp Month Day
1      41     190  7.4   67     5   1
2      36     118  8.0   72     5   2
3      12     149 12.6   74     5   3
(以下省略)
> transform(airquality, Ozone = -Ozone) # Ozone 欄の符号を変える
    Ozone Solar.R Wind Temp Month Day
1     -41     190  7.4   67     5   1
2     -36     118  8.0   72     5   2
3     -12     149 12.6   74     5   3
(以下省略)
> transform(airquality, new = -Ozone, Temp = (Temp-32)/1.8) # new 欄を加え、Temp 欄をスケーリング
    Ozone Solar.R Wind     Temp Month Day  new
1      41     190  7.4 19.44444     5   1  -41
2      36     118  8.0 22.22222     5   2  -36
3      12     149 12.6 23.33333     5   3  -12
 (以下省略)
# ベクトル Ozone とその対数値からなるデータフレーム
>      attach(airquality)
>      transform(Ozone, logOzone = log(Ozone)) # marginally interesting ...
      x logOzone
1    41 3.713572
2    36 3.583519
3    12 2.484907
 (以下省略)

(29) ベクトルの成分の系統的な置き換え、car パッケージの recode 関数 --2005.01.20

> (x <- rep(1:3, 3) )
[1] 1 2 3 1 2 3 1 2 3
> recode(x, "c(1,2)='A'; else='B'")      # x 中の 1,2 を "A", その他を "B" で置き換える
[1] "A" "A" "B" "A" "A" "B" "A" "A" "B"
> recode(x, "1:2='A'; 3='B'")            # x 中の 1,2 を "A", 3 を "B" で置き換える
[1] "A" "A" "B" "A" "A" "B" "A" "A" "B"
> recode(x, "c(1,2)='A'")                # 選択肢がなければそのままとされる(が他とあわせて文字化される)
[1] "A" "A" "3" "A" "A" "3" "A" "A" "3"

(28) 関数のソースコードを見る、methods(), getS3method() 関数 -- from r-help, 2005.01.18

R の多くの関数のソースコードは直接参照できない

> kruskal.test              # kruskal.test 関数のソースは?
function (x, ...)
UseMethod("kruskal.test")
<environment: namespace:stats>

# kuruskal.test のメソッド(実体)関数は実は二種類ある(そして隠匿されている)!
> methods(kruskal.test)
[1] kruskal.test.default* kruskal.test.formula*
 Non-visible functions are asterisked
 
> kruskal.test.default
Error: Object "kruskal.test.default" not found
> kruskal.test.formula
Error: Object "kruskal.test.formula" not found

# 隠匿された関数のソースを見るには getS3method() 関数を使う
> getS3method("kruskal.test", "default") # kruscal.test.default 関数のソース
function (x, g, ...)
{
    if (is.list(x)) {
        if (length(x) < 2)
            stop("x must be a list with at least 2 elements")
(途中省略)
    class(RVAL) <- "htest"
    return(RVAL)
}
<environment: namespace:stats>

> getS3method("kruskal.test", "formula")
function (formula, data, subset, na.action, ...)
{
    if (missing(formula) || (length(formula) != 3))
        stop("formula missing or incorrect")
    (途中省略) 
    y$data.name <- DNAME
    y
}
<environment: namespace:stats>

単に三重コロン演算子を使ってもよい

> stats:::kruskal.test.default
function (x, g, ...)
{
   if (is.list(x)) {
       if (length(x) < 2)
           stop("x must be a list with at least 2 elements")
   (途中略)
   class(RVAL) <- "htest"
   return(RVAL)
}
<environment: namespace:stats>

(27) 画面数に応じて適当な mfrow 値を計算する関数、n2mfrow() -- 2005.01.13

> n2mfrow(8)  
[1] 3 3           # プロット数が 8 なら mfrow = c(3,3) が適当(という意味) 
> n <- 5
> x <- seq(-2, 2, len = 51)
> op <- par(mfrow = n2mfrow(n))  # 基本的使い方(5つのプロットを一つの画面にレイアウト)
> for (j in 1:n) plot(x, x^j, main = substitute(x^exp,
                      list(exp = j)), type = "l", col = "blue")
> par(op)

(26) R には未実装(しかし実装が予定・期待されている)の関数、もしくは既存関数の未実装引数を知らせる機能がある、help(NotYetImplemented) を参照 -- 2005.01.13

ユーザーからの contribution を呼びかける意味もあるらしい。

> plot.mlm              # おそらく mixed linear model のプロットメソッド関数?
function (x, ...)
.NotYetImplemented()
<environment: namespace:stats>
> barplot(1:5, inside = TRUE)  # barplot 関数には未実装の inside 引数がある
Warning message:
argument 'inside' is not used (yet) in: .NotYetUsed("inside", error = FALSE)

(25) 与えられたベクトル、因子のすべての組合せから成るデータフレームを作る、expand.grid 関数 -- 2005.01.12

本来何に使うためにあるのか不明だが、ツールとして使う場面はありそう。
→ データフレームを読み込んで集計・統計解析を行う関数を作る際の、デバッグ用ダミーデータを発生させるのに使えそうですね。

## example(expand.grid) より
## 身長、体重、性別のすべての組合せ
> expand.grid(height = seq(60, 80, 5), weight = seq(100,
              300, 50), sex = c("Male", "Female"))
   height weight    sex
1      60    100   Male
2      65    100   Male
3      70    100   Male
4      75    100   Male
5      80    100   Male
6      60    150   Male
7      65    150   Male
8      70    150   Male
9      75    150   Male
10     80    150   Male
(途中省略)
49     75    300 Female
50     80    300 Female

(24) sink 関数よりも簡便なファイルへの出力関数、capture.output 関数 -- 2005.01.11

例えば

attach(some_data)
(some_data を使った作業)
detach(some_data) 

という一連の作業は

with(some_data を使った作業, data=some_data) 

と簡便化できるように、

sink(some_file)
(何かを出力する作業)
sink() 

というルーチンは

capture.output((何かを出力する作業), file = some_file)

と簡便化できる。出力が終るとファイルへの接続は自動的にクローズされる。もしファイル引数を省略すると文字列として R コンソールに出力される。 ?capture.output には次のようなもっと驚くような使い方が紹介されている:

## on Unix with enscript available
ps <- pipe("enscript -o tempout.ps","w") # pipe connection 
## example(glm) の出力を ps ファイルに落す
capture.output(example(glm), file=ps)    # connection ps への出力 (lpr 等でハードコピーできる)
close(ps)

(23) バイトベクトル、raw() -- 2005.01.10

R は raw ベクトル(成分が 1 byte データであるベクトル)を操作する関数群を持つ。

> raw(8)                     # 長さ 8 の raw ベクトル
[1] 00 00 00 00 00 00 00 00
> x <- "A test string"
> (y <- charToRaw(x))        # 文字列を raw ベクトルに変換
 [1] 41 20 74 65 73 74 20 73 74 72 69 6e 67     # 16進表現
> is.vector(y)
[1] TRUE
> rawToChar(y)               # 元の文字列に戻す 
[1] "A test string"
## 関連関数     
raw(length = 0)
as.raw(x)
charToRaw(x)
rawToChar(x, multiple = FALSE)
rawShift(x, n)
rawToBits(x)
intToBits(x)
packBits(x, type = c("raw", "integer"))

(22) ゼロ行列もしくはベクトルを作る関数、mat.or.vec() -- 2005.01.10

> mat.or.vec(3, 1)   # ゼロベクトル
[1] 0 0 0
> mat.or.vec(3, 2)   # 3x2 ゼロ行列
     [,1] [,2]
[1,]    0    0
[2,]    0    0
[3,]    0    0
> mat.or.vec(1,3)    # 1x3 ゼロ行列
     [,1] [,2] [,3]
[1,]    0    0    0

(21) ある値 x が増加数列 v のどこの間にあるかを計算する関数、findInterval() -- 2005.01.09

> v <- c(1,3,7,10,19)        # テスト用増加列
> x <- 4
> (i <- findInterval(x, v))
[1] 2
> c(v[i],x,v[i+1])           # x=4 は 区間 [v[i], v[i+1]) の間にある
[1] 3 4 7
> x <- 3
> (i <- findInterval(x, v))
[1] 2
> c(v[i],x,v[i+1])           # x=3 は区間 [v[i], v[i+1]) の間にある
[1] 3 3 7
> x <- 20
> (i <- findInterval(x, v))
[1] 5
> c(v[i],x,v[i+1])           # x が大き過ぎる場合
[1] 19 20 NA
> x <- 0
> (i <- findInterval(x, v))
[1] 0
> c(v[i],x,v[i+1])           # x が小さ過ぎる場合
[1] 0 1

注意:関数は引数 x についてベクトル化されている。 注意:データ X = c(X[1].X[2],...,X[n]) に対し、findInterval(t, sort(X))/n は X の経験分布関数の t における値となる。


(20) ベクトル(従って行列、配列) x 中の指定添字位置を他の値で置き換える、repalce 関数 (from r-help, 2005.01.05)

ただし x 自身を書き換えるわけではないので注意。必要なら x 自身に再付値する。

> x=1:10
> replace(x, which(x<4), NA)
 [1] NA NA NA  4  5  6  7  8  9 10
> x=matrix(sample(c(NA,1:4), 16, replace=TRUE, p=rep(1,5)/5), 4,4)
> x
     [,1] [,2] [,3] [,4]
[1,]   NA    4    1    4
[2,]    3   NA    1    2
[3,]    1   NA   NA    3
[4,]    4    4    1    3
## 行列中の NA 値を 0 で置き換える
> replace(x, is.na(x), 0)
     [,1] [,2] [,3] [,4]
[1,]    0    4    1    4
[2,]    3    0    1    2
[3,]    1    0    0    3
[4,]    4    4    1    3
## 行列中のすべての NA 値を順に 100,200 で置き換える
> x=matrix(sample(c(NA,1:4), 16, replace=TRUE, p=rep(1,5)/5), 4,4)
> replace(x, is.na(x), c(100,200))
     [,1] [,2] [,3] [,4]
[1,]    2    4    4    1
[2,]    3    1    2  200
[3,]  100    1    4  100
[4,]    2    4    3  200

(19) 二つの因子のすべての組合せを作る関数 interaction() -- (from r-help, 2005.01.05)

> a <- sample(c(0,1,2), 100, replace=TRUE, p=c(1,1,1)/3)
> b <- sample(c(0,1,2), 100, replace=TRUE, p=c(1,1,1)/3)
> x <- data.frame(A=a, B=b)
> attach(x)
> interaction(A,B)  ## すべてのケースに対する A, B の組合せ(因子)を計算
  [1] 1.2 0.0 1.0 0.2 1.2 2.0 0.0 0.0 2.2 2.2 0.1 2.1 1.2 0.2 1.0 0.0 2.1 2.0
 [19] 0.2 0.1 1.2 2.0 1.0 0.2 2.0 1.2 2.0 0.0 1.0 0.2 0.2 2.2 1.1 2.2 1.1 0.1
 [37] 1.0 2.0 2.0 2.2 1.0 1.1 0.0 1.2 0.2 0.0 0.2 1.0 1.2 2.1 0.0 1.1 1.0 0.1
 [55] 2.0 2.1 0.1 0.1 0.2 1.0 2.1 1.0 1.1 0.2 1.0 0.0 1.0 0.2 2.2 0.2 0.2 0.2
 [73] 2.0 2.2 1.1 0.0 0.2 1.0 1.0 0.1 2.2 0.1 0.0 0.1 2.0 0.0 0.0 2.2 2.2 1.2
 [91] 2.1 1.0 0.2 2.0 1.0 0.1 1.1 1.0 0.1 1.2
Levels: 0.0 1.0 2.0 0.1 1.1 2.1 0.2 1.2 2.2

## 例えば table 関数と組合せて、A, B の各組合せの度数を計算
> table(interaction(A,B))
0.0 1.0 2.0 0.1 1.1 2.1 0.2 1.2 2.2  # A=2,B=2 のケースは 10 組あった
 13  17  11  11   7   6  16   9  10

(18) R オブジェクトが占める(近似的)メモリサイズ, object.size() 関数 -- (2005.01.04)

よほど大規模な作業を行なわない限り、気にする必要はないでしょうが。

> object.size(letters)
[1] 1068
> object.size(ls)
[1] 9284
## base パッケージ中のオブジェクトのメモリサイズ10傑
> z <- sapply(ls("package:base"), function(x) object.size(get(x, 
             envir = NULL)))
> as.matrix(rev(sort(z))[1:10])
                   [,1]
library          108088
[<-.data.frame    72144
loadNamespace     60076
rbind.data.frame  54896
read.table        54792
merge.data.frame  42324
data.frame        42212
seq.POSIXt        35316
eigen             34956
write.table       30332
> x <- rnorm(1000000)
> object.size(x)
[1] 8000028
> x <- rnorm(10000000)
> object.size(x)
[1] 80000028

(17) 先頭文字を与えて一時ファイル名を作る関数、tempfile() -- 2005.01.02

具体的使用例は ファイルを読み込む tips 集(暫定版)参照

> tempdir()           # 一時ファイルを置いておく既定のディレクトリ(Linux の場合)
[1] "/tmp/RtmpFVxBmG"
> tempfile("abc")                  # 先頭文字 "abc" で一時ファイル名を作る
[1] "/tmp/RtmpFVxBmG/abc2ae8944a"
> tempfile("abc")                  # 同じセッションでは毎回異なるファイル名になる
[1] "/tmp/RtmpFVxBmG/abc238e1f29"
> tempfile("abc")
[1] "/tmp/RtmpFVxBmG/abc46e87ccd"
> tempfile("abc")
[1] "/tmp/RtmpFVxBmG/abc3d1b58ba"
> tempfile("abc", "~/temp")        # ディレクトリ名も指定する例
[1] "~/temp/abc625558ec"

(16) 二つのオブジェクトが(誤差範囲内で)等しいかどうか確認する「正しい」方法 -- 2004.12.30

数値の比較には計算機による実数処理にかかわる気づきにくい落し穴がある

> 0.3 == 0.7 - 0.4       # 一見 TRUE になりそうだが
[1] FALSE
> options(digits=22)
> 0.7 - 0.4         
[1] 0.2999999999999999333866  # 確かに 0.3 とは異なる
> 0.3
[1] 0.2999999999999999888978  # 一方(0.3 は 0.3 ではない!)

★ 0.3 を 16進数で表すと 0.4CCCCCCCC… となり,R が 0.3 と表示するのもある意味理解できないことではあるでしょうが。

数値オブジェクトの比較には all.equal() 関数を使う
all.equal() 関数は(既定では) .Machine$double.eps^0.5 以内の差は無視する
二つの数値ベクトル(したがって数値行列・配列の)の比較にも使える

>  .Machine$double.eps^0.5
[1] 1.490116119384765625000e-08
> all.equal(0.3, 0.7 - 0.4)
[1] TRUE

all.equal 関数は TRUE でないときは FALSE ではなく食い違い情報を返すので要注意

> all.equal(pi, 355/113)
[1] "Mean relative  difference: 8.491368e-08"

identical() 関数も数値オブジェクトの比較に使えるが all.equal() とは微妙な差がある。
identical() 関数は必ず論理値を返す

> 1 == as.integer(1)
[1] TRUE 
> all.equal(1, as.integer(1)) 
[1] TRUE
> identical(1, as.integer(1))      # 1 は実数、 as.integer(1) は整数でモードが異なるから
[1] FALSE
> all.equal(0.3, 0.7-0.4)
[1] TRUE
> identical(0.3, 0.7-0.4)    # 厳密な比較
[1] FALSE

いったいどうしたら良いのという時は(結構面倒だが)

# 二つの数値オブジェクトの同等性を(誤差範囲内)で確認する正しい構文
> x <- 1
> y <- 0.99999999999
> identical(all.equal(x, y), TRUE) 
[1] TRUE

(15) 「:」演算子の罠 -- from r-help. 2004.12.30

等差数列を作るのに便利な「:」演算子だが、うっかりすると思わぬことが。身に覚えのある人も多いはず。単項・二項演算子間の優先順位に気をつけよう。

> sequence <- 1:5
> 1:length(sequence) 
[1] 1 2 3 4 5
> 1:length(sequence)-1   # 1:(length(sequence)-1) のつもりでいるとやばい
[1] 0 1 2 3 4
> 1:length(sequence)+1
[1] 2 3 4 5 6
> 1:length(sequence)*2
[1]  2  4  6  8 10
> 1:length(sequence)/2
[1] 0.5 1.0 1.5 2.0 2.5
> 1:length(sequence)^2  # ところが冪乗演算子は 「:」よりも優先順位が高いから厄介
 [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

P. Dalgaard 氏はプログラム中では「:」演算子の代わりに seq() 関数を使用することを勧めている。なぜなら

> seq(length=length(sequence))       # 当然 1:length(sequence) と同じ
[1] 1 2 3 4 5
> seq(length=length(sequence)-1)     # 先ほどの思い込みミスが避けられる
[1] 1 2 3 4
> seq(l=length(sequence))            # 当然省略記法が使える
[1] 1 2 3 4 5 
> seq(along=sequence)                # 1:length(sequence) と同じ
[1] 1 2 3 4 5
> seq(a=sequence)                    # こうすれば 1:length(sequence) より簡潔だし
[1] 1 2 3 4 5
> seq(2, length=length(sequence)-1)  # 2:(length(sequence)-1) 
[1] 2 3 4 5
> seq(3)
[1] 1 2 3
> seq(0)     #  但しこれは危険(S との互換性のためにこうなっている)
[1] 1 0

汎用的関数のプログラムではベクトルがたまたま空になるケースも考えておく必要が起きる

> sequence <- numeric(0)             # もしベクトルが空だったら?
> 1:length(sequence)                 # length(numeric(0)) は 0 とされるから!
[1] 1 0
> seq(length=length(sequence))       # これらなら間違いにすぐ気づく
numeric(0)
> seq(along=sequence)                # 同じく
numeric(0)
> seq(length=length(sequence)-1)     # この場合もエラーとなるから安全
Error in seq.default(length = length(sequence) - 1) :
        Length must be non-negative number

(14) 重複する要素の添字を返す、duplicated() --- 2004.12.26

ただし重複する要素の二度目以降の位置の添字だけを返す

> x <- c(9:20, 1:5, 3:7, 0:8)
> x
 [1]  9 10 11 12 13 14 15 16 17 18 19 20  1  2  3  4  5  3  4  5  6  7  0  1   2
[26]  3  4  5  6  7  8
> duplicated(x)   # 同じものが二度目以降に現れた位置が TRUE とされる
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE FALSE FALSE FALSE  TRUE
[25]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
> x[duplicated(x)]
 [1] 3 4 5 1 2 3 4 5 6 7
> x[!duplicated(x)]        # unique(x) と同じだが、 unique(x) の方が効率的
 [1]  9 10 11 12 13 14 15 16 17 18 19 20  1  2  3  4  5  6  7  0  8
> unique(x)                
 [1]  9 10 11 12 13 14 15 16 17 18 19 20  1  2  3  4  5  6  7  0  8
> x[duplicated(x)]         # 重複する要素(二度目以降)
 [1] 3 4 5 1 2 3 4 5 6 7

(13) unstack() 関数は因子の値に基づきベクトルを分割する。stack() 関数はその逆操作 -- 200412.26

# テスト用データフレーム(因子毎のデータ数が等しい時)
> (x <- data.frame(x=6:1, ind=c(rep("a",2),rep("b",2),rep("c",2))))
  x ind
1 6   a
2 5   a
3 4   b
4 3   b
5 2   c
6 1   c
> (y <- unstack(x))   # 因子毎に分類 
  a b c
1 6 4 2
2 5 3 1
> stack(y)             # 再構成、x に戻る
  values ind
1      6   a
2      5   a
3      4   b
4      3   b
5      2   c
6      1   c
> stack(y, select=-c)    # c グループを除いて再構成 
  values ind
1      6   a
2      5   a
3      4   b
4      3   b
# 因子毎のデータ数が異なる場合 
> (x <- data.frame(x=6:1, ind=c(rep("a",3),rep("b",2),rep("c",1))))
  x ind
1 6   a
2 5   a
3 4   a
4 3   b
5 2   b
6 1   c
> y <- unstack(x)   # 因子毎のデータ数が異なる時はリストに分類
$a
[1] 6 5 4
$b
[1] 3 2
$c
[1] 1
> stack(unstack(x))   # この際も stack() は unstack() の逆操作
  values ind
1      6   a
2      5   a
3      4   a
4      3   b
5      2   b
6      1   c

(12) 数列の並行生成、少し面白い -- from r-help, 2004.12.23

> x = 1:4
> y = 2*(2:5)
> mapply(seq, x, y)
[[1]]
[1] 1 2 3 4                 # seq(1,4)

[[2]]
[1] 2 3 4 5 6               # seq(2,6)

[[3]]
[1] 3 4 5 6 7 8             # seq(3,8)

[[4]]
[1]  4  5  6  7  8  9 10    # seq(4,10)

(11) y のどれかに一致する x の要素の添字を返す -- match(x, y) 関数、2004.12.23

> x <- c(1,1,2,4,3,2,2,1)
> match(x, 1:4)                # x[3] は y[2] と一致するので 2 が返される等々
[1] 1 1 2 4 3 2 2 1
> match(x, 4:1)
[1] 4 4 3 1 2 3 3 4
> match(x, 1:2)                # 一致しない x[i] に対しては NA が返される 
[1]  1  1  2 NA NA  2  2  1
> z <- c("one", "tow", "three")
> z[match(x, 1:3)]
[1] "one"   "one"   "tow"   NA      "three" "tow"   "tow"   "one"
> match(letters[1:3], letters[1:5])
[1] 1 2 3
> c("one","tow","three")[match(x, 1:3)]
[1] "one"   "one"   "tow"   NA      "three" "tow"   "tow"   "one"
# R の内部二項演算子  %in% はつぎのように定義されている
"%in%" <- function(x, table) match(x, table, nomatch = 0) > 0
# y 中に無い x 要素の添字を返す二項演算子(?match より)
"%w/o%" <- function(x,y) x[!x %in% y]
# 使用例 
> (1:10) %w/o% c(3, 7, 12)
[1]  1  2  4  5  6  8  9 10

(10) show 関数 (2004.12.18)

初級者用 Q&A の記事で「発見した」関数。通常知る(使う)必要の無い関数で、どうも、print メソッドの実体らしい(?) print.default() は showDefault() を呼び出すと解説がある。

> x <- 1:10
> x            
# (恐らく)内部的には print(x) (print.default(x)) が呼び出され、
# 更に showDefalt(x) が呼び出される、というややこしい仕組み(?)
 [1]  1  2  3  4  5  6  7  8  9 10
> print(x) 
 [1]  1  2  3  4  5  6  7  8  9 10
> print.default(x)
 [1]  1  2  3  4  5  6  7  8  9 10
> show(x)
 [1]  1  2  3  4  5  6  7  8  9 10
> showDefault(x)
 [1]  1  2  3  4  5  6  7  8  9 10

(9) if 文のびっくりする使い方 (2004.12.17)

関数ではありませんが、今日 if 文のびっくりする使い方の例を目にしたのでメモしておきます。 言われてみれば当たり前ですが、こんな使い方は普通思いつかない(はず?)

> x=c(3,4)
> y <- "a"; x * if(y=="a") 1 else 2  # x*1 になる
[1] 3 4
> y <- "b"; x * if(y=="a") 1 else 2  # x*2 になる
[1] 6 8
> y <- "a"; x[if(y=="a") 1 else 2]   # x[1] になる
[1] 3
> y <- "b"; x[if(y=="a") 1 else 2]   # x[2] になる
[1] 4

追加コメントします
この,マイルドな使い方は

a <- if (i %% 2) "odd" else "even"

普通の言語では

if (i %% 2) a <- "odd" else a <- "even"

というようになるのにという趣旨でしょう。(ALGOL は R と同じような構文を許していたかな?C では a = 論理式 ? 式1 : 式 2 で,これは R での a <- ifelse(論理式,式1, 式2) と同じですね。ただ,if と ifelse は,論理式にベクトルを許すかどうかの違いがあるし,C では当然論理式にベクトルは使えないが)
if (論理式) 式1 else 式2 自体が式ですから,その結果を付値して何が不都合かということでしょうね。

もう少しひねくれているのは

> x = 1
> a <- (if (x %% 2) sin else cos)(3) 
> a
[1] 0.14112
> x=2
> a <- (if (x %% 2) sin else cos)(3) 
> a
[1] -0.9899925
> sin(3)
[1] 0.14112
> cos(3)
[1] -0.9899925

if ~ else で関数名を割り当てている。びっくりもんですか?(使いましたよ。。)

↑ これはもう「こんなこと許しても本当にいいの?!」というレベルですね。以下から分かるように確かにこの if 文は関数オブジェクト(正確には内部関数 sin へのエントリポイント)を返していますね。驚きです。

> x=1; if (x %% 2) sin else cos
.Primitive("sin")
> .Primitive("sin")(2)
[1] 0.9092974

(8) 指定したモードと長さのベクトルを作る関数 vector() (2004.12.14)


> (out <- vector("list", 4)) # リストからなる長さ4のベクトルを作る
[[1]]
NULL  # NULL で初期化されている
[[2]]
NULL
[[3]]
NULL
[[4]]
NULL
> out[2]  # 第二成分はリスト
[[1]]
NULL

> (mout <- matrix(vector("list", 4), 2,2)) # リストの行列!
     [,1] [,2]
[1,] NULL NULL
[2,] NULL NULL
> mout[1,2]       # [1,2] 成分は確かにリスト
[[1]]
NULL

(7) 文字列を指定した長さに分割した文字列ベクトルを作る、関数 strwrap()、(2004.12.05)

但し、分割は空白文字位置で行なわれ、単語を分割することはない。また余分な空白は取り除かれる。

> x="Each character string in the input is first split into paragraphs" # テスト用文字ベクトル
> strwrap(x,40)  # 40 文字で分割
[1] "Each character string in the input is"
[2] "first split into paragraphs"
> strwrap(x,20)  # 20 文字で分割
[1] "Each character" "string in the"  "input is first" "split into"
[5] "paragraphs"
> strwrap(x,10)  # 10 文字で分割   
[1] "Each"       "character"  "string"     "in the"     "input is"
[6] "first"      "split"      "into"       "paragraphs"
# 一続きの単語が分割長より長ければ、そのまま分割されない。
> strwrap(x,7)
 [1] "Each"       "character"  "string"     "in"         "the"
 [6] "input"      "is"         "first"      "split"      "into"
[11] "paragraphs"
# 特に長さ 1 を指定すれば、単語に分割される。
> strwrap(x,1)
[1] "Each"       "character"  "string"     "in"         "the"
[6] "input"      "is"         "first"      "split"      "into"
# 次の構文はグラフィックスタイトルを複数行に折り返す時便利。指定長で分割の上、改行文字 "?n" を挟んで再結合している。
> sapply(lapply(x, strwrap, 15), paste, collapse = "?n")
[1] "Each?ncharacter?nstring in the?ninput is?nfirst split?ninto?nparagraphs"
> sapply(lapply(x, strwrap, 20), paste, collapse = "?n")
[1] "Each character?nstring in the?ninput is first?nsplit into?nparagraphs"
> sapply(lapply(x, strwrap, 40), paste, collapse = "?n")
[1] "Each character string in the input is?nfirst split into paragraphs"

(6) ベクトル、行列、データフレームをある性質で分割する、split(x,f)

> x <- matrix(rnorm(16), 4,4)  # テスト用行列
> x
          [,1]       [,2]       [,3]       [,4]
[1,]  1.261265 -1.1751028  1.0451261  0.2706234
[2,] -1.269447  0.7282239 -0.7784738 -1.1093583
[3,]  2.649739  1.2154682  1.1069235 -1.1391058
[4,] -1.024763  0.2041062  0.1663222  0.5010538
> split(x, col(x))  # 列番号で分割
$"1"
[1]  1.261265 -1.269447  2.649739 -1.024763 # x[,1] と同じ
$"2"
[1] -1.1751028  0.7282239  1.2154682  0.2041062
$"3"
[1]  1.0451261 -0.7784738  1.1069235  0.1663222
$"4"
[1]  0.2706234 -1.1093583 -1.1391058  0.5010538
> split(x, row(x))  # 行番号で分割            
$"1"
[1]  1.2612648 -1.1751028  1.0451261  0.2706234 # x[1,] と同じ
$"2"
[1] -1.2694475  0.7282239 -0.7784738 -1.1093583
$"3"
[1]  2.649739  1.215468  1.106924 -1.139106
$"4"
[1] -1.0247627  0.2041062  0.1663222  0.5010538
> L3 <- LETTERS[1:3]
> (d <- data.frame(cbind(x=1, y=1:10), fac=sample(L3, 10, repl=TRUE)))
   x  y fac  # テスト用データフレーム
1  1  1   A
2  1  2   A
3  1  3   A
4  1  4   A
5  1  5   B
6  1  6   C
7  1  7   A
8  1  8   A
9  1  9   B
10 1 10   B
> split(d, d[,3]) # fac に関し分割、結果はデータフレームのリスト
$A
  x y fac
1 1 1   A
2 1 2   A
3 1 3   A
4 1 4   A
7 1 7   A
8 1 8   A
$B
   x  y fac
5  1  5   B
9  1  9   B
10 1 10   B
$C
  x y fac
6 1 6   C 

(5) もう一つのソート関数、sort.list()

ソートといえば sort() 関数だが、sort.list() 関数は文字ベクトルでもソートできる。 100000 以下の整数だけからなるベクトルなら sort.list(x, method="radix") は sort(x) よりも早くて正確らしい。

> x <- sample(letters, 1e6, replace=TRUE) # 小文字アルファベット10万個のベクトル
> gc(); system.time(o <- sort.list(x))
[1] 11.82  0.02 11.87  0.00  0.00
> xx <- sample(1:100000, 1e7, replace=TRUE) # 長さ一千万個整数ベクトル
> gc(); system.time(o <- sort.list(xx, method="radix")) 
[1] 5.27 0.16 5.53 0.00 0.00
> gc(); system.time(o <- sort(xx)) 
[1] 11.03  0.22 11.28  0.00  0.00

(4) namespace を操作する関数達

R のパッケージ(例:すべてのベースパッケージと推奨パッケージ)には namespace と呼ばれる、ユーザー定義オブジェクトや他のパッケージのオブジェクトとオブジェクト名の衝突を避けるためメカニズムがある(らしい、知ったかぶり)。例えば関数 printCoefmat のコードを表示すると最後に <environment: namespace:stats> と表示されるが、これはこの関数がパッケージ stats の namespace に登録されている(と言うべきなのか?)ことを示している。こうした関数を通常の仕方で書き換えても、実際には依然 stats の namespace に登録されたものが呼び出されるので、無効になる。こうした「保護された」関数を操作するための関数には次のようなものがある。普通のユーザーが使うようなものではなく、R の基本関数の仕様をどうしても変えたい時にのみ慎重に使う。具体的な例が Q&A コーナーの記事「経済学における統計的有意の記号」で議論されている。

以下は M. Maechler 氏による namespace の解説(r-help 記事より)

Utility functions to access and replace the non-exported functions
in a namespace, for use in developing packages with namespaces.
    getFromNamespace(x, ns, pos = -1, envir = as.environment(pos))
    assignInNamespace(x, value, ns, pos = -1, envir = as.environment(pos))
    fixInNamespace(x, ns, pos = -1, envir = as.environment(pos), ...)
please do learn about namespaces;  many good R packages (incl.
all standard and recommended packages) nowadays make use of
namespaces, and we expect even more in the future.

One big advantage of 'namespacing a package' is that
you can make sure that your functions call your other functions
even if they have the same name as someone's functions in
another package.  Another big "pro" is that you can use as many
``package-internal'' objects {"helper-functions" in my useR talk}
as you want: They won't be visible to the user and not clutter
the global namespace.
One side effect of namespaces: We usually do not export (S3)
methods such as silhouette.default() : 

A user should call silhouette(...) and the method dispatching
chooses the corresponding method.

Now, to see "hidden" objects, you can 
1) use  getAnywhere()	        { getAnywhere("silhouette.default") }
2) use  <namespace>:::<object>  { cluster:::silhouette.default } 
and for S3 methods,
3) getS3method("<generic>", "<class>")
   i.e., in this case getS3method("silhouette", "default")

(3) ベクトルから同じ数が引き続く回数とその数を取り出す、rle()

いわゆる run length を計算する。inverse.rle() はその逆関数で、元のベクトルを再現する。

> (x <- rep(6:10, 1:5))  # テスト用ベクトル
 [1]  6  7  7  8  8  8  9  9  9  9 10 10 10 10 10
> (y <- rle(x))          # x の run length を計算
Run Length Encoding
  lengths: int [1:5] 1 2 3 4 5   # 同じ数が繰り返される回数
  values : int [1:5] 6 7 8 9 10  # 繰り返されている数 (例えば最後に 10 が 5 回続く)
> inverse.rle(y)                 # y から x を復元する 
 [1]  6  7  7  8  8  8  9  9  9  9 10 10 10 10 10

(2) 行列・配列の特定の次元を軸にあるベクトルを順に操作する、関数 sweep()

> x = matrix(1:16, 4,4)
> x
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
> y= 1:4                  # sweep されるベクトル
> sweep(x,1,y)            # x の各列から y を引き去る(既定動作) 
     [,1] [,2] [,3] [,4]
[1,]    0    4    8   12
[2,]    0    4    8   12
[3,]    0    4    8   12
[4,]    0    4    8   12  
> sweep(x,2,y)            # x の各行から y を引き去る(既定動作) 
     [,1] [,2] [,3] [,4]
[1,]    0    3    6    9
[2,]    1    4    7   10
[3,]    2    5    8   11
[4,]    3    6    9   12
> sweep(x,1,y, "*")       # x の各列に y を掛ける 
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    4   12   20   28
[3,]    9   21   33   45
[4,]   16   32   48   64

(1) 文字グラフィックス、関数 symnum()

関数 symnum() は数値・論理値ベクトル、配列を適当な区間毎に対応する文字記号に置き換え表示する。有意水準記号 ***,**,** の表示の際使われている。区間と対応記号をオプションで指定できる。

> x <- cor(matrix(rexp(30, 1), 5, 18)) # 指数乱数の行列の相関行列 
> symnum(x)

 [1,] 1
 [2,] * 1
 [3,]     1
 [4,] . . , 1
 [5,] . . , . 1
 [6,] . . ,   . 1
 [7,] 1 *   . . . 1
 [8,] * 1   . . . * 1
 [9,]     1 , , ,     1
[10,] . . , 1 .   . . , 1
[11,] . . , . 1 . . . , . 1
[12,] . . ,   . 1 . . ,   . 1
[13,] 1 *   . . . 1 *   . . . 1
[14,] * 1   . . . * 1   . . . * 1
[15,]     1 , , ,     1 , , ,     1
[16,] . . , 1 .   . . , 1 .   . . , 1
[17,] . . , . 1 . . . , . 1 . . . , . 1
[18,] . . ,   . 1 . . ,   . 1 . . ,   . 1
attr(,"legend")  # コーディングに使われた区間と対応文字の説明
[1] 0 ` ' 0.3 `.' 0.6 `,' 0.8 `+' 0.9 `*' 0.95 `B' 1

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-03-25 (土) 11:19:16