東京工業大学の間瀬茂です.拙著「Rプログラミングマニュアル」の出版から7年がたちました.当時Rのバージョンは2.6,2.7でした.この四月にRはついにバージョン3の大台に到達しました(Rがささやかにネットに公開されてから丁度20年の節目なのは偶然でしょうか?).

拙著の基本方針は「Rのプログラミング機能を網羅的に紹介する」でした.今回のバージョンアップにあたって,果たしてRがどのくらい変わっているのか気になり調べています.結果は想像以上で,拙著は今でも有効であるが,新規に加わった機能を紹介しないのは,「網羅的」という趣旨に反するということで,現在出版社と相談の上改訂版をだす作業中です(私の現役 最後の仕事になります).出版時期はまだ未定です.

改訂作業中に特に「!!!」と思った点が幾つかありました.これらは単にこれまで私が気づかなかっただけのもありますし,新たに付け加わったものもあるようです.この頁でその内の幾つかを改訂版原稿からピックアップして紹介したいと思います.

もし,私の誤解,理解不足がありましたら,コメントいただけると幸いです.但し,参考になったコメントは改訂版に反映させて頂きますので,ご了解お願いします.

注意:頁タイトルが誤解を与えるかもしれませんが,ここに書いていることは私の本の初版と予定改訂版の差分(つまりこの7年間分,一部は前回私が見落としたものが混ざっているかもしれません)のうち私が面白いと感じたものです.バージョン3ではじめて登場したオブジェクトという意味ではありません.

(独り言:Rバージョン4が出る頃,私はまだ生きているかな...と考えました.)



「apply関数族は簡潔な既述ができるが,実体は for ループを隠蔽しているだけで,速度向上は望めない」とこれまでは理解していましたが,最近(?)のRでは実際に速くなるようです.

lapply, vapply, sapply, rapply 関数は内部関数を用いて計算するため,処理内容によっても異なる可能性があるが,単なるループ処理よりは計算速度自体が向上すると思われる.それ以外もバイトコンパイルされている分計算速度は向上すると思われる.

以下は sapply 関数と for ループを使った同値な処理をそれぞれ100回繰り返した際の時間計測結果である.僅かであるが sapply 関数を使ったほうが速いことが分かる.

foo <- function(x) x # つまらない関数
X <- numeric(1e6)    # 念のため,予め変数を用意
# sapply関数を用いた処理を100回繰り返す
Y <- sapply(1:100, function(i) system.time(X <- sapply(1:1e6,foo)))
# そのユーザ時間の要約
> c(fivenum(Y[1,]),mean(Y[1,]),sd(Y[1,]))
[1] 2.048 2.240 2.280 2.388 2.604 2.303 0.09771
# forループを用いた処理を100回繰り返す
Z <- sapply(1:100, function(x) system.time(for(i in 1:1e6) X[i] <- foo(i)))
# そのユーザ時間要約
> c(fivenum(Z[1,]),mean(Z[1,]),sd(Z[1,]))
[1] 2.548 2.592 2.606 2.624 2.684 2.609 0.02449
  • OS やメモリーなどの環境によるのかもしれませんが,以下のような結果になりますので,参考までに。 示されたプログラムでは system.time(X <- sapply(1:1e6,foo)) にされていますが,私の環境では gc が起動するようで,sapply, for の真の所要時間が測れているかどうかよくわからないので {gc(); gc(); system.time(X <- sapply(1:m, foo))} のようにして計測しました。また,対象要素数が 1e6 の場合と,例えば 1e4 の場合では様相が違うようなので,以下のようなシミュレーションプログラムで速度を見ました。
    foo <- function(x) x # つまらない関数
    sim <- function(n, m) {
    	X <- numeric(m)    # 念のため,予め変数を用意
    	Y <- replicate(n, {gc(); gc(); system.time(X <- sapply(1:m, foo))})[1, ]
    	Z <- replicate(n, {gc(); gc(); system.time(for(i in 1:m) X[i] <- foo(i))})[1, ]
    	return(list(s=Y, f=Z))
    }
    得られた測定データの最小値,最大値,平均値,標準偏差を以下に示します(n は対象要素数)。
              n=1e4               n=1e6
         sapply    for       sapply    for      
    min  0.00700   0.01600   2.11900   1.73000
    max  0.01900   0.02300   2.91500   2.76100
    mean 0.00868 < 0.01838   2.45032 > 1.84305
    sd   0.00140   0.00183   0.12012   0.22043
    n=1e4 のとき(左図)は for は sapply の倍くらい時間がかかる。しかし,n=1e6 のとき(右図)では形勢逆転で sapply のほうがかなりおそい。ヒストグラムを見ると for では時間がかなりかかる場合が若干あるために平均値が引き上げられているので,実行時間の差は平均値で見る以上に大きいようだ。 -- 2013-10-23 (水) 18:49:18
    fig1.png fig2.png
  • 早速のコメントありがとうございます.実行環境により様子が随分と違うようですね.ちなみに私の場合1e4,1e5回の繰り返しでもsapplyの方が速いという結果になります.それと相対的に速い場合と遅い場合の二つの山ができる傾向があるようですね. -- 間瀬? 2013-10-23 (水) 20:48:34
    n=1e5 sapply
    0.132 0.140 0.144 0.148 0.180 0.145 0.00770
    n=1e5 for ループ
    0.248 0.256 0.260 0.264 0.292 0.260 0.00780
    n=1e4 sapply
    0.0120 0.0120 0.0120 0.0140 0.0160 0.0130 0.00174
    n=1e4 forループ
    0.0240 0.0240 0.0240 0.0280 0.0320 0.0253 0.00205
  • あ, sapplyはlapplyのただのラッパーです. 性能のバラツキはキャッシュミスによるものが殆どでしょうからsapply(内側:計測対象)がやってる余計な事(!)を省いてlapplyを使って計測すればもうちょっと安定すると思います. Rのforやwhile等のループは以前はもっと遅くって問題視()されていましたが現在はかなり(でも当時と比べて)高速化されています. まあlist処理なのでapplyはキャッシュミスが出やすいですね. 一方forの方は付値主体(あまりリスト処理が無い)なので純粋にコードの為の実行時間が見えているのだと思います. -- 2013-10-24 (木) 11:23:35

「Rシステム」はあるが「R言語」は無い.

現在「R言語」という呼び方がしばしばみられますが,私は二重の意味でこれは好ましくないと考えています.何よりもRは統計システム・環境であり,言語機能は即統計解析機能ではなく,単に計算機言語と捉えるのは矮小化していると思うからです.更にRの言語機能はS言語のフリーな移植(エンジン,方言)であり,開発者達自身がRの言語機能はS言語であると明記しています(以前は無かったので,恐らく世界的にR言語と呼ぶ風潮が広まっているためと思われます).

例えば,公式マニュアル An Introduction to R には

”a well developed, simple and effective programming language (called ‘S’)  
which includes conditionals, loops, user defined recursive functions and input
and output facilities. (Indeed most of the system supplied functions are 
themselves written in the S language.)”

公式FAQ集の記述では

"We can regard S as a language with three current implementations 
or“engines”, the “old S engine” (S version 3; S-PLUS 3.x and 4.x), 
the “new S engine”(S version 4; S-PLUS 5.x and above), and R. 
Given this understanding, asking for “the differences between R and S”
really amounts to asking for the specifics of the R implementation of
the S language, i.e., the difference between the R and S engines.”

Tierny氏の compiler パッケージが推奨パッケージになった

R開発メンバーのTierny氏の compiler パッケージのおかげで,Rの実行時間が知らぬ間に数割早くなっていたようです.現在のR(及びほとんどの貢献パッケージ)の内部関数を利用する関数以外は予めバイトコンパイルされているか,読み込み時にバイトコンパイルされるようです.

関数 enableJIT は実行時(Just-In-Time)コンパイルを指示します.JIT機能を許可すれば,それ以降全ての関数は最初の実行時に黙ってバイトコンパイルされます(最適化レベルを上げれば必ず速くなるわけでもない).level 引数が0なら許可しない.1ならクロージャは最初の使用時にコンパイルされる.もし2なら,クロージャ自体もコンパイルされる.もし3なら,全てのループが実行前にコンパイルされる. 関数 compilePKGS はパッケージを読み込む際に,含まれる関数をバイトコンパイルするかどうかを指示する.

# lapply関数の初期の定義
> la1 <- function(X,FUN,...) {  
    FUN <- match.fun(FUN)
    if (!is.list(X)) 
      X <- as.list(X)
    rval <- vector("list",length(X))
    for(i in seq(along=X)) 
      rval[i] <- list(FUN(X[[i]],...))
    names(rval) <- names(X) 
    return(rval) }
# 時間のかかるまずいコード
> tmp <- function(n) {
    x <- numeric(0)
    for(i in 1:n) 
      x <- c(x,runif(1))
    return(x) }
# コンパイルする
> la1c <- cmpfun(la1)
> tmpc <- cmpfun(tmp)
> la1c
function(X,FUN, ...) {
  FUN <- match.fun(FUN)
  if (!is.list(X)) X <- as.list(X)
  rval <- vector("list",length(X))
  for(i in seq(along=X)) rval[i] <- list(FUN(X[[i]],...))
     names(rval) <- names(X) 
  return(rval) }
<bytecode: 0xa9a3388> # バイトコンパイルされている証拠
> tmpc
function(n) {x <- numeric(0); for(i in 1:n) x <- c(x,runif(1)); x}
<bytecode: 0xa64e254>
# 以下は la1(1:100,tmp) を1万回実行したユーザ時間の要約
   Min. 1st Qu.  Median  Mean 3rd Qu.  Max. 
  0.016   0.020   0.020 0.021  0.020  0.024 # JITコンパイル無し
  0.004   0.012   0.012 0.011  0.012  0.024 # enableJIT(1)
  0.004   0.012   0.012 0.011  0.012  0.016 # enableJIT(2)
  0.008   0.012   0.012 0.011  0.012  0.016 # enableJIT(3)

バイト型データと関連関数の導入

Rの原子オブジェクトに新しく(何時から?)バイト(raw)型データが付け加えられました.「生型」では落ち着きが悪いですね.そのままRAW型と書いたり,ロー型と訳されている例がありますが意訳してバイト型と訳そうと考えていますが,すでに良い訳語があるでしょうか.一つのバイト型データは16進数(0,1,2,...,9,a,b,c,d,e,f,A,B,C,D,E,Fでも良い)の対で表現されます.一方でバイト型データのビット列(2進法)表現があり,00, 01 (それぞれ2進法の0,1を表す)の列で表現されます.一つのバイト型データは長さ32のビット列で表現されます.とりあえず私には使う場面は思いつきませんが.

# 255は16進法ではff
> x <- as.raw(255); x
[1] ff  
> str(x)
raw ff
> x <- "A test string"
# xのバイト列表現
> y <- charToRaw(x); y
[1] 41 20 74 65 73 74 20 73 74 72 69 6e 67  
# 文字列に戻す
> rawToChar(y) 
[1] "A test string"
# 整数8のビット列表現
> intToBits(8L)
 [1] 00 00 00 01 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
[26] 00 00 00 00 00 00 00

隠し変数 .Last.value

Rのトップレベル(つまりRコンソールのプロンプト > に対して直接入力された)の命令を内部的に評価した値を隠し変数 .Last.value に記録していることはご存じの方も多いかと思います(おそらく結果を print 関数でコンソール出力するため).従って巨大なオブジェクトを操作すると,それと同じサイズの .Last.value も作られることになるので注意がいります.一見メモリにまだ十分余裕があるはずなのにメモリ不足エラーになる理由の一つがこれと思われます(その他に変形途中に一時的ファイルが作られることが挙げられます).ガベージコレクション gc() を実行すると,結果として .Last.value も gc() の出力に書き換えられ,直近のそれの分メモリに大きな隙間が生じる可能性があります.このためガベージコレクションは 続けて二回実行することが勧められています.

> x <- log(10); x
[1] 2.302585
> .Last.value
[1] 2.302585
# library関数の返り値は不可視
> library(splines)
# .Last.valueには記録
> .Last.value
[1] "splines" "stats"
--以下省略--
# 巨大なオブジェクトx
> x <- rnorm(1e6)
# .Last.valueにも記録
> identical(x,.Last.value)
[1] TRUE
# この時点で.Last.valueはすでに更新されている
> str(.Last.value)
TRUE

ユーザレベルでも .Last.value 変数の便利な使い方があります.ある時間のかかる計算をうっかり結果を変数に付値するのを忘れて実行してしまったことありませんか(私はよくあります).この時は慌てて実行中断せず,直後に res <- .Last.value 等とすれば結果を保存できます.

> sin(1)
[1] 0.841471
> res <- .Last.value; res
[1] 0.841471

マルチコア並列処理パッケージ parallel が推奨パッケージになった

Rの推奨パッケージに新しく付け加わった並列処理パッケージ parallel はマルチコアCPUを持つ一台の計算機で並列処理を行うことを想定しており,二つの貢献パッケージ multicore と snow の内容を受け継ぐ(両者は実装の仕方が異なる).最近の計算機のCPUはマルチコア化されており,一台の計算機で複数のプロセスを同時に動作させられる.一つのコアには複数のプロセスを割り当てることもできるため,実際の子プロセス数は物理的コア数を超えても良い(が結局無理をすることになる).

# クラスタを作り並列処理する例
> require(parallel)
# 4コアのクラスタを生成
> cl <- makeCluster(4)
# 変数xxと関数fooを定義
> xx <- 1; foo <- function(y) xx + y
# 小プロセスにそれを移出,オブジェクト名は文字列で与える
> clusterExport(cl,c("xx","foo"))
# sapply関数の並列版parSapplyを用いfoo(1:10)をクラスタに割り当てて計算
> parSapply(cl,1:10,foo)
 [1]  2  3  4  5  6  7  8  9 10 11
# 2変数関数を定義し,小プロセスにそれを移出
> bar <- function(x,y) x + y
> clusterExport(cl,"bar")
# 多変数apply関数であるmapply関数の並列版clusterMap.以下のリスト成分はクラスタが分担計算
> clusterMap(cl,bar,c(a=1,b=2,c=3),
                    c(10,0,-10))
$a     # つまりbar(a=1,10)
[1] 11
$b
[1] 2
$c
[1] -7
> stopCluster(cl) # 最後にクラスタを閉じる

関数 pvec はクラスタの生成,apply 関数風の並列計算と結果の統合,そしてクラスタの終了、という一連の処理を自動で行う便利な関数である.オプション mc.cores で子プロセス数を指定できる(既定では2).

次の結果は日付文字列一万個を as.POSIXct 関数で日付オブジェクトに変換することを1,000回繰り返した際の経過時間の要約と標準偏差である.pvec を使った並列計算が高速であること,実行毎に所要時間は相当異なることがわかる.子プロセスが増えればそれだけ早くなるわけでもない(並列化のオーバーヘッドのためと思われる).この例では子プロセス数8(使用計算機の総コア数)の場合が一番高速であった.なお,並列処理の時間計測では,経過時間を使うのが好ましい.ユーザ時間は計算中子プロセスの管理だけをしている親プロセスのそれを計測しており,かなり小さくなる.平均計算時間比は329:330:185:133:117:130:113:102:100:112:109:349.

# 日付文字列一万個のベクトルdatesをas.POSIXctで日付オブジェクトに変換する時間を計測
# すべて経過時間を示す
> dates <- sprintf('%04d-%02d-%02d',as.integer(2000+rnorm(1e4)),
                   as.integer(runif(1e4,1,12)),as.integer(runif(1e4,1,28)))
> N <- 1:1000
# 非並列化版
> x <- sapply(N, function(i) system.time(as.POSIXct(dates))[3])
> summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1750  0.1770  0.1780  0.1795  0.1810  0.2530 
# 以下子プロセス数n=1,2,...,10,100で並列計算の経過時間を計測
> x <- sapply(N, function(i) system.time(pvec(dates,as.POSIXct,mc.cores=n))[3])
# summary(x)の結果
       最小値  第1四分偏差 中央値  平均    第3四分偏差  最大値
       0.175   0.177      0.178   0.180   0.181       0.253 # 非並列版 
 n=1   0.176   0.177      0.178   0.180   0.182       0.195 # 親プロセスを含めれば並列版
 n=2   0.094   0.099      0.100   0.101   0.102       0.130 
 n=3   0.066   0.071      0.072   0.072   0.073       0.113 
 n=4   0.052   0.057      0.059   0.064   0.068       0.091 
 n=5   0.059   0.070      0.072   0.071   0.073       0.079 
 n=6   0.051   0.060      0.061   0.062   0.063       0.095 
 n=7   0.051   0.053      0.054   0.056   0.056       0.090 
 n=8   0.048   0.050      0.052   0.055   0.059       0.072 # 最も高速
 n=9   0.053   0.058      0.061   0.061   0.063       0.086 
 n=10  0.052   0.057      0.059   0.060   0.061       0.079 
 n=100 0.148   0.181      0.192   0.190   0.201       0.260
  • 推奨パッケージ内のbootはsnow及びmulticoreを使えるようになってます -- 2013-10-24 (木) 11:36:40

long vectors 機構の登場

大規模データ処理の必要性から,バージョン3.0のRよりlong vectors 機構が採用され,実数換算で32ペタバイト相当にもなる最大 2^{52}=4.503610^{15} (約4,504兆)のサイズのオブジェクトが可能になった.Rの公式マニュアル「R Internals」の最新版の12章には long vectors の実装に関する解説がある.現在全ての原子的 (raw, logical, integer, numeric, complex,character) ベクトル,リスト そして 表現式が対象になっている.

しかし,これは十分なメモリ(long vectors オブジェクトは結局使用出来るRAMのサイズを超えることはできない)を搭載した64ビットマシンで,使用OSがそれをサポートしている場合に限られる.行列,配列の各々の次元は依然 2^{32}-1 以下の必要があり, 一次元ベクトルの長さは依然として 2^{32}-1 以下である. そうしたベクトルに対する数値演算には,対応関数が long vectors をサポートしている 必要がある.詳細は help(”long vectors”) を参照せよ.

long vectors に対する添字は(整数を表す)倍精度実数であり,length 関数は整数でなく,倍精度実数値を返す.

実際のところ,私はそんなにハイスペックの計算機を持っていませんので,自分では試せないでいます.どなたか使用経験のある方がいらっしゃったら感想をお聞かせください.


国際化の進展

中間さんのRの日本語化に端を発し,Rの国際化が始まったことはご存知かと思います. 現在のRでは特に UTF エンコーディングのサポートが充実してきているようです.

\unnnn, \u{nnnn}, \Unnnnnnnn, \U{nnnnnnnn} は 4桁もしくは8桁の16進数を用いたUnicodeによる文字表現を与える.

# キリル(ロシア)文字の例
> c("\u0414","\u{0411}")
[1] "Д" "Б"

Rの推奨パッケージ tools 中には,非アスキー文字や多バイト文字の各種エンコーディング(特にUTF-8)における文字とその整数表現を扱う関数がある.utf8ToInt? は長さ1のUTF-8文字列を整数に変換する.intToUtf8 はその逆(既定では単一文字列に変換)である.

# intToUtf8(946)でもよい
> intToUtf8(0x03B2L) 
[1] "β"
# 元に戻る
> utf8ToInt(intToUtf8(0x03B2L))
[1] 946
> intToUtf8(945:950)
[1] "αβγδεζ"
> intToUtf8(945:950,TRUE)
[1] "α" "β" "γ" "δ" "ε" "ζ"
# エキゾチックな半角文字がいっぱい
> x <- intToUtf8(as.integer(
    c(160:383,0x0192,0x02C6,0x02C7,0x02CA,
      0x02D8,0x02D9,0x02DD,0x200C,0x2018,
      0x2019,0x201C,0x201D,0x2020,0x2022,
      0x2026,0x20AC)),multiple=TRUE)
# 見やすくする
> matrix(x,ncol=16,byrow=TRUE)
---出力省略---

群盲象を撫でる

R3.0の基本パッケージ中には総計1,264(内部関数は488)個のオブジェクトがある.CRANから入手できるバージョン3.0.0の基本,標準そして推奨パッケージのヘルプ文章を網羅した3,604頁に及ぶPDFファイル R-fullrefman.pdf がある.こうした目録が必要なほどRシステムは巨大になっている.Rに関しては「群盲象を撫でる」,「車輪を二度発明する」という故事成語が決して冗談では無い.CRANに登録されているアドオンパッケージ数も4,700余りであり,CRAN未登録分を含めれば5,000を超えることは確実と思われる.


noquote関数と文字列出力

cat関数とprint関数は文字列の出力形式(引用符の有無)が異なることはご存知かと思います.dataframeの出力でも文字変数は引用符無しで出力されます(意外に言われないと気づかない?).私も長い間気づきませんでしたが,この動作の背後にあるのがnoquote関数のようです.

noquote関数は文字列ベクトルにクラス属性 "noquote" を加え,print関数による出力の際に文字列引用符無しで出力(cat関数では既定動作)されるようにする.鈎括弧演算子 [ ] による要素抽出は文字列ベクトルと同様にできる.

> cat(letters[1:3],"\n") # cat関数は文字列を引用符なしで出力
a b c 
> print(letters[1:3]) # print関数では引用符付き
[1] "a" "b" "c"
> print(noquote(letters[1:3])) # noquote関数を使えばprint関数でも引用符がつかない
[1] a b c
# クラス'noquote'を持つ文字列ベクトルになる
> str(noquote(letters[1:3]))
Class 'noquote'  chr [1:3] "a" "b" "c"
> noquote(letters)[1:3]
[1] a b c

cat関数による出力に制御文字を使う

cat関数では制御文字 "\n" が改行の意味を持つことはよく知られていますが,その他の制御文字も使えます.

cat関数に使用でき,次の出力位置を制御できる制御文字(エスケープシーケンス)には以下のようなものがある:

 \a  ベルを鳴らす
 \b  一文字分戻る
 \f  現在の位置から一行進む
 \n  改行し行先頭に
 \r  現在の行の先頭に戻る
 \v  垂直タブ,一定行分進む
 \t  水平タブ,次のタブ位置から行う

これらの制御文字は cat 関数によるコンソール出力で有効であるが,print 関数では無視される.

> cat("abc\n")
abc  
# 途中で二回連続改行
> cat("abc\n\ndef\nijk\n")
abc

def
ijk
# 行頭に戻って次を出力し最後に改行
> cat("abc\rdef\n")
def  
# abcの後,一文字戻ってdefを出力後改行
> cat("abc\bdef\n")
abdef 
# print関数では無視される
> print("abc\bdef")
[1] "abc\bdef"
# 2文字戻って次を出力し改行
> cat("abc\b\b\n)
adef 
# 3文字戻って次を出力
> cat("abc\b\b\b\n")
def 
# abc出力後,行先頭に戻りdefを出力
> cat("abc\rdef\n")
def
# 水平タブを二回使用
> cat("abc\tdef\tijk\n")
abc     def     ijk
# 垂直タブを二回使用
> cat("abc\vdef\vijk\n")
abc
   def
      ijk
# 途中で二回行下げ
> cat("abc\fdef\fhij\n")
abc
   def
      hij
abdehij
> cat("abc\vdef\rijk\n")
abc
ijkdef
> cat("abc\vdef\v\rijk\n")
abc
   def
ijk
# 垂直タブを二回連続使用
> cat("abc\vdef\v\v\rijk\n")
abc
   def

ijk

画面への出力で行頭に戻って出力する制御文字 "\r" は,何回もの計算結果を繰り返しコンソールに出力する際,大量の出力で画面が流れ去れないようにするのに有用である.以下の例では,2秒毎に数字1,2,...,100を同じ行に更新出力し,最後の出力100だけが画面に残る.

# 途中出力は上書きされ,最後の出力だけが残る
> for(i in 1:100) 
    {Sys.sleep(2) 
     cat("\r",i,"番目の出力")}; cat("\n")
 100 番目の出力
# 更に空白を除く
> for(i in 1:100) 
    {Sys.sleep(1)
      cat("\r",i,"\b番目の出力")}
    cat("\n")
 100番目の出力

行列・配列のコンパクトな表示

Rオブジェクトのコンソールへの出力は多くの場合好ましい形式が選ばれるが,場合によると煩わしいこともある.

行列や配列のコンソールへの(print関数による)表示の際は暗黙のうちに次元からなる次元名が付く.以下の関数 no.dimanames は次元表示無しにコンパクトに表示(空文字ベクトルからなるリストで次元名を与えている)する.

> no.dimnames <- function(a) {
    d <- list(); l <- 0
    for(i in dim(a)) 
      d[[l <- l + 1]] <- rep("", i)
    dimnames(a) <- d
    a }
# 既定の表示
> 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
# 次元名なしに表示
> no.dimnames(X)
 1 5  9 13
 2 6 10 14
 3 7 11 15
 4 8 12 16
# 配列にも使える
> X <- array(1:8,c(2,2,2)); X
, , 1
     [,1] [,2]
[1,]    1    3
[2,]    2    4
, , 2
     [,1] [,2]
[1,]    5    7
[2,]    6    8
> no.dimnames(X)
, , 
 1 3
 2 4
, , 
 5 7
 6 8
  • X <- matrix(1:16,4,4); X は (X <- matrix(1:16,4,4)) のように例示される方が良いと思います -- 2013-10-25 (金) 21:33:15
  • 念のため,このような機能は,一番最初 Version 0.49 Beta (April 23, 1997) からありました。 -- 2013-10-25 (金) 21:35:58

Lisp(Reduce)風の構文を持つ関数 Reduce, Filter, Find, Position, Map, Negate

Lisp(Reduce)にある機能を真似た関数が幾つか登場しています.Lisp の達人 Tierny 氏の作でしょう.以下で紹介するように面白い使い方ができます.なお,いまだにRはLisp(Scheme)で書かれていると信じている人がいて驚いたことがあります.RはSchemeのアイデアを取り入れているのは事実ですが,Lispで書かれているわけではありません.

関数 Reduce は二項演算子をベクトル x 中の要素に逐次適用した結果をかえす.Filter は(論理値)関数 f の値が TRUE になるベクトル x 中の要素を抜き出す.Find はそうした要素の先頭と末尾を与える.Position はそうした要素の位置を与える.Map は mapply 関数のラッパー関数であり,関数をベクトルの要素に順に適用した結果をリストで返す.Negate は与えられた(論理値)関数の否定を作る.

以下で定義する関数 add, cadd はそれぞれ汎用加算関数と汎用累積和関数であり,和が定義できるものなら何でも可.

 > add <- function(x) Reduce("+",x)
 > add(list(1,2,3))
 [1] 6
 > cadd <- function(x) Reduce("+",x,accumulate=TRUE)
 > cadd(seq_len(7))
 [1] 1 3 6 10 15 21 28
 > cadd(list(1:3,2:4))  # ベクトルの累積和
 [[1]]
 [1] 1 2 3
 [[2]]
 [1] 3 5 7
 > cadd(list(matrix(1,2,2), matrix(2,2,2)))  # 行列の累積和
[[1]]
     [,1] [,2]
[1,]    1    1
[2,]    1    1
[[2]]
     [,1] [,2]
[1,]    3    3
[2,]    3    3

連分数を計算する.

# 連分数を計算する関数
> cfrac <- function(x) Reduce(function(u,v) u+1/v,x,right=TRUE)
> cfrac(c(3,7,15,1,292)) # 円周率を近似する連分数
[1] 3.141593
> cfrac(c(2,1,2,1,1,4,1,1,6,1,1,8)) # exp(1)を近似する連分数
[1] 2.718282

関数の重複適用.

> Funcall <- function(f,...) f(...)
# log(exp(acos(cos(0))を計算
> Reduce(Funcall,list(log,exp,acos,cos),0,right=TRUE)
[1] 0
# 黄金比の連分数近似
> cfrac(rep.int(1,31))
[1] 1.618034
# 関数 t |-> (t+x/t)/2 の不動点としてsqrt(x)を近似計算
> asqrt <- function(x,n) Iterate(function(t) (t+x/t)/2,n)
# 初期値を正の数とする
> asqrt(2,30)(10)
[1] 1.414214
# 初期値を負の数とする
> asqrt(2,30)(-1)
[1] -1.414214

マニュアル,ヘルプ文章も進化している

Rの魅力の一つが充実したマニュアルとオブジェクト毎のヘルプ文章ですが,これらも折に触れて改訂がされているようです.使い慣れた関数のヘルプ文章も一度読みなおすと新しい発見があるかもしれません.関数自体も仕様が変更されていることもあります.


Rの商業的利用に関するR Foundationの見解

Rの商業的利用もかっての,上司に内緒で会社のパソコンにインストールして使っている,という笑い話のような状況から,現在ではかなり公に使われているようです.実際,企業のメイン業務でRを使っている,検討しているという話も聞きます.ビッグデータブームを背景にグーグル等の世界的企業がRの使用を後押ししていることも追い風になっているようです.ご存知のように,RはGPL条項の下で公開されているオープンソースシステムです.Rの商業的利用に付いてR Foundation自体がどのような見解を持っているのかはっきりしませんでしたが,公式 FAQ に次のような解説がありました.要するにRコードを含むプログラム自体で儲けるには制限があるが,それ以外の商業的利用はご勝手にという事のようです.なお,R Foundation自体は昔からこういう法律的問題にはノータッチの姿勢です.気になるのはすでに複数存在するらしい,Rに独自の機能を付け加えたシステムを販売している組織です.普及版は無償提供,高機能版は有償というスタンスのように見えますが.これは Linux の世界でも同様の問題ですが.

2.11 Can I use R for commercial purposes?

R is released under the GNU General Public License (GPL), version 2. 
If you have any questions regarding the legality of using R 
in any particular situation you should bring it up with your legal counsel. 
We are in no position to offer legal advice.

It is the opinion of the R Core Team that one can use R for commercial
purposes  (e.g., in business or in consulting). The GPL, like all Open Source 
licenses, permits all and any use of the package. It only restricts 
distribution of R or of other programs containing code from R. This is made 
clear in clause 6 (“No Discrimination Against Fields of Endeavor”) of 
the Open Source Definition:

    The license must not restrict anyone from making use of the program 
in a specific field of endeavor. For example, it may not restrict the program
from  being used in a business, or from being used for genetic research. 
 
It is also explicitly stated in clause 0 of the GPL, which says in part

    Activities other than copying, distribution and modification are not
covered by this License; they are outside its scope. The act of running 
the Program is not restricted, and the output from the Program is covered 
only if its contents constitute a work based on the Program. 

Most add-on packages, including all recommended ones, also explicitly allow
commercial use in this way. A few packages are restricted to “non-commercial 
use”; you should contact the author to clarify whether these may be used or 
seek the advice of your legal counsel.

None of the discussion in this section constitutes legal advice. The R Core
Team does not provide legal advice under any circumstances. 

パッケージ data.table

データフレームはRのデータ形式の中心であり,多くの統計処理関数は引数としてデータフレームを想定しています.一方で,データフレームは行列風に見せかけたリストで,計算効率の面では行列操作等に比べれば格段に遅くなることが悩みの種でした.コミュニティーで開発されているパッケージ data.table はいずれデータフレームに代わる可能性を持つデータテーブルオブジェクトを提供します.但し操作方法が独自で慣れるのに手間がかかります.同じような趣旨で開発されているパッケージも複数あるようです.

パッケージ data.table (Rの推奨パッケージではない)はデータフレームの代替オブジェクトであるデータテーブル(data.table) オブジェクトを提供する.データフレームは便利なデータ形式であるが,本質はリストで操作時間は行列等と比べれば格段に遅い.データテーブルはより速い操作を与える.data.table クラスは data.frame クラスを継承しており,データフレームが要求されるシステム関数でもそのまま使える.データテーブルオブジェクトは操作が速いだけではなく,部分抽出,変形,データテーブル同士の結合等に対する簡潔で多彩な操作方法を与える.また変形に際しても一時的コピーを作らず,元オブジェクトそのものを変形する(参照渡し)工夫がされており,データフレームに比べればメモリへの負担が少ないことも利点である.

以下は同じ内容の大きなデータフレームとデータテーブルの一部を取り出す例である.ユーザ時間比は 15:5:1,経過時間比は 24.7:7.7:1 であった.キーを設定したデータテーブルに対する操作は高速なことがわかる.

> x <- sample(1:10, 1e6,rep=TRUE)
> y <- sample(letters[1:10], 1e6,rep=TRUE)
> DF <- data.frame(x=x,y=y)
> DT <- as.data.table(DF)
# データフレーム版
> z <- subset(DF,y=="a", select=x)
# データテーブル版
> z <- DT[y=="a"]
# キー指定したデータテーブル版
> setkey(DT,y)
> z <- DT["a"]

次の例はデータテーブルの操作がデータフレームのそれとは異なり,内部的にオブジェクトのコピーを利用せず.既存オブジェクトを「その場で(参照)」加工(せいぜい列単位のコピーのみ)していることを示す.tracemem 関数でマークされたオブジェクトは,内部的にコピーされるとそれを報告する.

> x <- sample(0:9,1e5,replace=TRUE)
> y <- sample(letters[0:9],1e5,replace=TRUE))
> z <- runif(1e5)
> DF <- data.frame(x=x,y=y)
> tracemem(DT)
[1] "<0x9d195e0>"
# データフレームの連結操作.2回の内部コピー
> DF <- cbind(DF,z=z)
tracemem[0x9d195e0 -> 0xcc2e068]: data.frame cbind cbind 
tracemem[0xcc2e068 -> 0xcc2c958]: data.frame cbind cbind 
> DF <- data.frame(x=x,y=y)
> tracemem(DF)
[1] "<0xad5ad58>"
# 同じ事を別の操作で.4回の内部コピー
> DF <- transform(DF,z=z)
tracemem[0xb36b928 -> 0xb367928]: do.call transform.data.frame transform 
tracemem[0xb367928 -> 0xb367a88]: do.call transform.data.frame transform 
tracemem[0xb367a88 -> 0xb367d68]: data.frame do.call transform.data.frame transform 
tracemem[0xb367d68 -> 0xb364658]: data.frame do.call transform.data.frame transform 
> DT <- data.table(x=x,y=y)
> tracemem(DT)
[1] "<0x8354770>"
# データテーブルに新しい列を加える.コピーされていない
> DT[,z:=z]

print関数のオプション zero.print

Rのコンソールへの出力は,実はprint関数の無数のメソッド関数が,オブジェクトの形式に合わせて好ましい書式に整形して出力しています.コンソールに単にオブジェクト名をタイプしたり,命令を実行すると結果が出力されるのも影でprint関数が出動しているわけですが,初心者はしばしば混乱するようですね.オプション zero.print="."を用いるとゼロをドットで表示します(それ以外の文字列も可).

> t1 <- round(abs(rt(200,df=1.8)))
> t2 <- round(abs(rt(200,df=1.4)))
# メソッドprint.table使用
> table(t1,t2) 
    t2
t1  0 1 2 3 4 5 6 7 8 10 17 21 30
  0 21 22 14 4 1 0 1 1 1 0 1 0 0
  1 25 21  7 3 4 2 1 1 1 1 0 0 0
  --途中省略--
  12 1 0 0 0 0 0 0 0 0 0 0 0 0
# 値0をドットで表現
> print(table(t1,t2),zero.print=".")
    t2
t1   0  1  2  3  4  5  6  7  8 10 17 21 30
  0 21 22 14  4  1  .  1  1  1  .  1  .  .
  1 25 21  7  3  4  2  1  1  1  1  .  .  .
  --途中省略--
  12 1 .  .  .  .  .  .  .  .  .  .  .  .
  • Version 0.60 Alpha (December 2, 1997) から使えるようになっていました。 -- 2013-10-25 (金) 21:57:14

こんなのあったけ? 行列・配列編

今回 R-fullrefman.pdf をあれこれ見ているうちに「こんなのあったけ?」と思う関数が結構ありました.少なくとも7年前には無かったはず...(?).

関数 rowsum(x,group,reorder=TRUE,...) (rowSumsではありません)は group 変数でグルーピングされた行の和を計算する.group 変数は行数と同じ長さのベクトルである.数値ベクトル,データフレームに対しても使える.ベクトル x は列数1の行列と見なされる. 結果は行数 length(unique(group)) で x と同じ列数を持つ行列である.他にも,プログラム内部で使うことを想定したという .colSums, .rowSums, .colMeans, .rowMeans という関数もあった.

> rowsum(1:5,c(1,2,2,2,2))  # 各列をグループ1,2に分けてグループ毎の総和
  [,1]
1    1   # グループ1の総和   
2   14   # グループ2の総和   
> rowsum(1:5,c(1,2,2,3,3))  # グループ数3
  [,1]
1    1   # 1
2    5   # 2+3
3    9   # 4+5
> 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
> rowsum(x,c(1,2,2,3))    
  [,1] [,2] [,3] [,4]
1    1    5    9   13   # グループ番号1の行の列和 (第1列そのまま)
2    5   13   21   29   # グループ番号2の行の列和 (第2,3列の和)
3    4    8   12   16   # グループ番号3の行の列和 (第4列そのまま)
> rowsum(x,c(1,2,2,2))
  [,1] [,2] [,3] [,4]
1    1    5    9   13  # グループ番号1の行の列総和 (第1列の和)
2    9   21   33   45  # グループ番号2の行の列総和(第2,3,4列の和)
> rowsum(x, c("a","b","a","b"))  # グループ変数は何でも良い
  [,1] [,2] [,3] [,4]
a    4   12   20   28
b    6   14   22   30

ベクトルは行列と異なり次元属性を持たないので nrow, ncol は NULL を返す.一方,関数 NROW, NCOL は行列・配列に対しては nrow, ncol と同じ値を返すが,ベクトルは as.matrix を用いて縦ベクトルと見なす.

> y <- 1:(3*4*5)
> c(NROW(y),NCOL(y))  # ベクトルに対するNROW,NCOL
[1] 60 1 

配列 x とその次元番号 MARGIN に対する slice.index(x,MARGIN) は行列に対する row, col 関数の配列版である.返り値は x と同じサイズの配列になる.

 > x <- array(1:24,c(2,3,4))
 > x1 <- slice.index(x,1)
 # x1[n,i,j]はn
 > x1[1,,]
      [,1] [,2] [,3] [,4]
 [1,]    1    1    1    1
 [2,]    1    1    1    1
 [3,]    1    1    1    1
 > x1[2,,]
      [,1] [,2] [,3] [,4]
 [1,]    2    2    2    2
 [2,]    2    2    2    2
 [3,]    2    2    2    2
> x <- array(1:8,c(2,2,2))
> x1 <- slice.index(x,1)  
> x2 <- slice.index(x,2)
> x3 <- slice.index(x,3)
> x[x2 == x1 & x3 == x1]   # 配列の一般化対角成分 x[1,1,1]とx[2,2,2]
[1] 1 8
  • Version 1.8.1 (2003-11-21) には既にありました。
    rowsum は Version 0.63.1 (December 5, 1998) からです。
    slice.index は Version 1.7.0 (2003-04-16) には既にありました。 -- 2013-10-24 (木) 21:15:42

こんなのあったけ? apply関数族編

vapply 関数は出力書式指定の sapply 関数だそうです.

関数 vapply(X,FUN,FUN.VALUE,...,USE.NAMES=TRUE) は sapply に似るが、引数 FUN.VALUE で出力書式をあらかじめ指定できる.安全であり、実効速度も速くなることがある.もし length(FUN.VALUE)==1 ならば X と同じ長さのベクトルが、さもなければ配列が返される.もし FUN.VALUE が array でなければ、結果は行数 length(FUN.VALUE)、列数 length(X) の行列になる.さもなければ次元が c(dim(FUN.VALUE),length(X)) の配列になる.もし FUN.VALUE が名前付きならばそれから配列の(次元)名前ラベルが、さもなければ適用関数から名前ラベルが作られる.行列の列名や、配列の最終次元の名前は sapply 関数に於けるように X から設定される.

> i39 <- sapply(3:9,seq); i39 # ベクトルのリスト
[[1]]
[1] 1 2 3
[[2]]
[1] 1 2 3 4
--途中省略--
[[7]]
[1] 1 2 3 4 5 6 7 8 9
> vapply(i39,fivenum,c(0,0,0,0,0))   # 一つの結果は長さ5のベクトルであることを明示
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]  1.0  1.0    1  1.0  1.0  1.0    1
[2,]  1.5  1.5    2  2.0  2.5  2.5    3
[3,]  2.0  2.5    3  3.5  4.0  4.5    5
[4,]  2.5  3.5    4  5.0  5.5  6.5    7
[5,]  3.0  4.0    5  6.0  7.0  8.0    9

こんなのあったけ? 作表関数編

関数 prop.table(x,margin=NULL) は分割表の各項目を margin に対する周辺和との比率で表す.margin は周辺を指定する添字(ベクトル)である.これは実際は sweep(x,margin,margin.table(x, margin),"/") である.もし margin が空であれば x/sum(x) を計算する.

> m <- matrix(1:4,2)
> margin.table(m,1)
[1] 4 6
> margin.table(m,2)
[1] 3 7
> prop.table(m,1)
          [,1]      [,2]
[1,] 0.2500000 0.7500000
[2,] 0.3333333 0.6666667
> prop.table(m,2)
          [,1]      [,2]
[1,] 0.3333333 0.4285714
[2,] 0.6666667 0.5714286
  • Version 1.8.1 (2003-11-21) には既にありました。
    Version 1.0.1 (April 14, 2000) からの関数です。当初は単に
    function (x, margin) 
    sweep(x, margin, margin.table(x, margin), "/")
    でした(今も大差ない) -- 2013-10-24 (木) 21:10:32
    > example(prop.table)
    
    prp.tb> m <- matrix(1:4, 2)
    
    prp.tb> m
         [,1] [,2]
    [1,]    1    3
    [2,]    2    4
    
    prp.tb> prop.table(m, 1)
              [,1]      [,2]
    [1,] 0.2500000 0.7500000
    [2,] 0.3333333 0.6666667

こんなのあったけ? タイマー編

setTimeLimit? はトップレベルの各計算(つまり,コマンド行入力命令による計算)に上限時間を設ける.もしある計算中で使われると,既定ではその計算及びそれ以降の計算に適用される.transient=TRUE ではその計算だけに適用される.上限に達するとエラーメッセージが出力され計算は中断される.setSessionTimeLimit? の実行はそれ以降のRセッションの継続時間に上限を設ける.上限に達するとエラーメッセージが出力され,上限設定はリセットされる.経過時間の測定が行われるのは,ユーザが実行を中断したければできるようなタイミング(及び Sys,sleep 関数実行中)に限られ,外部言語サブルーチン実行中は行われない.

> setSessionTimeLimit(10,10)
> for(i in seq(1e4)) {
    cat(i); Sys.sleep(1)}
12345678910 以下にエラー Sys.sleep(1) : セッション時間が経過して上限に達しました 
> {setTimeLimit(10,10); 
   for(i in seq(1e4)) 
     {cat(i); Sys.sleep(1)}}
12345678910 以下にエラー Sys.sleep(1) : 時間が経過して上限に達しました 

gc.time はガベージコレクションを行い,必要だった時間を返す.

# 巨大なベクトルを作り,消してからガベージコレクション時間を計測.少しずつ必要時間が大きくなる?
> sapply(1:10, function(i) {x <- rnorm(1e8); rm(x); gc.time()}
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,] 1.574 1.584 1.602 1.612 1.622 1.644 1.654 1.664 1.686 1.696  # ユーザ時間
[2,] 1.026 1.072 1.122 1.168 1.218 1.268 1.318 1.364 1.414 1.464  # システム時間
[3,] 1.097 1.138 1.190 1.231 1.272 1.326 1.367 1.407 1.459 1.501  # 経過時間
[4,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[5,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
  • R version 2.9.1 (2009-06-26) に既にありました.
    なかったんじゃないかなとか,知らなかった,あったんだっけ?など曖昧な基準で判断しないようにしたほうが良いと思います。 -- 2013-10-24 (木) 22:36:41
  • ご指摘ありがとうございます.ここで「あったけ?」と言っているのはRプログラミングマニュアル執筆当時には,私は気づかなかったという意味です.逆に教えていただけると参考になりますが,あるオブジェクトが何時からRに導入されたのか調べる便利な方法があったら教えてください.私はそうした方法を思いつかず,悩んでおります. -- 間瀬? 2013-10-24 (木) 23:11:08
  • 全てのバージョンのRをインストールしておけばよい。あるバージョンでのエラーがいつ直ったかなども実際にやってみれば明らか。 -- 2013-10-24 (木) 23:50:58

遅延コピー(とでもいうんだろうか?)

help(tracemem)を見ていていまさらだが気づいたこと(当たり前のことでしたら御免).b <- a とすると a のコピー b が作られるわけだが,必ずコピー(つまり a,b が異なるメモリ領域に対応)が作られるわけでは無く,実際に a もしくは b の内容が変更されるまでメモリ上でのコピーは実行されない.確かにそれはその方が賢いと勝手に感動しました.

> a <- 1:10
> tracemem(a)      # aが指すメモリ領域がメモリ上でコピーされたら通知するようにする
[1] "<0x990d438>"
> b <- a
> b[1:3]           # この段階ではa,bは同じメモリ領域を共有していれば済むから本当のコピーは不要
[1] 1 2 3
> sum(b)           # 同じく
[1] 55
> a[1] <- 10       # a が変更されたので実際のコピーが必要になった 
tracemem[0x990d438 -> 0x990d558]:  # 二度コピーされているのは一時変数を作成しているから?
tracemem[0x990d558 -> 0x8ea59e8]: 
> untracemem(a)
> a <- 1:10
> tracemem(a)
[1] "<0xa3e9cb0>"
> b <- a
> b[1] <- 10       # bが変更されればやはり実際のコピーが必要になる
tracemem[0xa3e9cb0 -> 0xa3e9c68]: 
tracemem[0xa3e9c68 -> 0xa2c5e08]: 
> untracemem(a)
> tracemem(a)
[1] "<0xa3e9cb0>"
> b <- a
> rm(a)            # コピー動作は無い,考えてみれば当然

S3, S4クラス,メソッド

S言語本家の Chambers 氏が開発メンバーに加わったせいがあると思うが,S4(そして関連するS3)関連の関数が充実してきた,というより混乱してきた,という印象がある.少なくともヘルプ文章を読んでも,何をどうすれば良いのかわかりにくい.一般ユーザが気楽に使えるものではないとの印象を受けるし,苦労して使うメリットもあまり感じられない.7年前にはいずれS4オブジェクトがS3オブジェクトを駆逐するのかと予想したが,現実はそうはなっていないようだ.Google のRコーディングスタイルでもS3オブジェクトが推奨されている.結局のところRは計算機言語の専門家用のシステムではなく,統計ユーザが使うものだと考えると,計算機言語としての完成度よりは,素人にも使いやすい「変態言語」のままが賢明なのかもしれない.


豆知識 #line directive

Rでは,#記号から行末まではコメントとされ実行時には無視される,ことになっていましたが,今は例外があるようです.ソースコードに含まれる

#line nn "ファイル名"

(#line は行頭から5文字にあるべき)はC言語のそれと同じく,parser にソースファイルとその行番号 nn を指示する directive とみなされるそうです.今後他の directive も導入される可能性があるようです.とりあえず,うっかり #line で始まるコメント行を入れても実害はない?

R news の記事の一部を引用

The #line directive

In some cases, R source code is written by a program,
not by a human being. For example, Sweave() extracts 
lines of code from Sweave documents before sending 
the lines to R for parsing and evaluation. To support 
such preprocessors, the R 2.10.0 parser recognizes 
a new directive of the form
#line nn "filename"
where nn is an integer. As with the same-named
directive in the C language, this tells the parser to
assume that the next line of source is line nn
from the given filename for the purpose of constructing
source references. The Sweave() function doesn’t
currently make use of this, but in the future, it (and
other preprocessors) could output #line directives
so that source references and syntax errors refer to
the original source location rather than to an inter-
mediate file.
The #line directive was a late addition to R
2.10.0. Support for this in Sweave() appeared in R 2.12.0

ローマ数字の話

Rには正整数をローマ数字に変換する関数があります.いまさらですが,ローマ数字では4,000以上は表現できないことを知り驚きました.あの古代ローマの壮大な建築物はどうやって建設したんでしょう.事情は中国,日本でも同じですが.

関数 as.roman(n) は正整数 n をローマ数字表記に変換する.ローマ数字には0や負の数は無い.3,899(MMMDCCCXCIX)までが一意的に表現可能であり,4,000以上の数は表現できない.四則演算ができ,結果はローマ数字である.

# ローマ数字表記
> x <- as.roman(c(1,5,10,100)); x
[1] I   V   X   C  
# ローマ数字には0や負の数は無い
> as.roman(c(-1,0))
[1] <NA> <NA>  
# クラス'roman'を持つ整数である
> str(x)
Class 'roman'  int [1:6] 1 5 10 100
# ローマ数字に対する四則演算例
> as.roman(10)+1
[1] XI    
> as.roman(10)-1
[1] IX
# ローマ数字では3,899までが一意的に表現可能
# 4,000以上の数は表現できない
> as.roman(3899); as.roman(3900)
[1] MMMDCCCXCIX
[1] <NA> 
# ただの整数に戻る 
> c(as.roman(3),as.roman(10))
[1]  3 10
  • http://www.unicode.org/charts/nameslist/n_2150.html によると, Archaic Roman numerals の中に Roman Numeral One Hundred Thousand (U2188) と言うのがありますから、古代は使っていたのではないでしょうか? -- 2013-10-28 (月) 12:19:44
  • as.roman のヘルプが,http://en.wikipedia.org/w/index.php?title=Roman_numerals&oldid=78252134 をひいていますが,そのページでは,それは「古いバージョンである」として,最新の http://en.wikipedia.org/wiki/Roman_numerals が挙げられています。そこに,Alternate forms という項目があります(旧版にも該当部あり)。
    なお,as.roman では,3899 までが表現可能ということになっていますが,その上限の由来は書かれていません。しかし,google で「3999 はどのようにあらあわすか」という質疑応答 3999 = MMMCMXCIX があったりはします。日本語版 wikipedia 参照。
    ローマ数値記法については,Project Euler での説明 も役に立ちそうです。 -- 2013-10-28 (月) 14:12:45


添付ファイル: filefig1.png 438件 [詳細] filefig2.png 466件 [詳細]

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