Rは統計解析のために作られたソフトだが、普通のプログラミング言語と同等な機能を備え、科学計算の強力な道具として利用可能である。また、基本的なグラフィックスを簡単な命令で実現することができるので、結果を視覚的に確認する道具としても利用価値が十分にある。統計解析者だけに占有させるのはもったいない。ここでは、統計解析に特化しない便利な使い方を解説する。
インストールは簡単、インストールのページを検索してダウンロードするだけ。
Rconsole 画面は次のようなもので、左下の「>」が入力待ちのプロンプトである。数式を入力すれば直ちに結果が表示される。
大量のデータ入力や、グラフィックス画面の自動保存などで扱いが簡単になる RStudio という統合環境があるので、それを使うと良い、というよりも、効率を考えるとほとんど必須の道具である。
R をインストールした後に、「rstudio install」で検索すると、rstudio.com というページを見いだすことができるので、自分のコンピュータのOSに合った Installer をダンロードする。日本語化されていないので、検索は英単語で。
次が RStudio の標準的画面である。左下が Rconsole で結果が文字で表示される。右下の「Plots」タブにはグラフィックスが表示され、「Help」には関数のオンラインマニュアルが表示される。右上は Rconsole に入力された変数や関数の値が記録される。左上はエディタで、プログラムを編集することができる。
RStudio の便利な機能（一部）
RStudio を立ち上げると、幾つかの画面（ペインと呼ばれる）が同時に表示される。そのうちの左下（デフォルト）がconsole ウィンドウになっている。（RStudio を使わない場合）R を立ち上げると、最初に「RGui」ウィンドウが表示され、 その中に「R console」ウィンドウが表示される。
console ウィンドウは、計算結果やエラーメッセージを文字で表示するための必須ウィンドウである。開始時にいろいろなメッセージが表示され、最後に「>」記号が表示されて入力待ちになる。この状態で、計算指示の命令文を入力しEnterキーを押すと、計算結果が表示される。たとえば次のように。
> 1+2*3
[1] 7
> 355/113 - pi [1] 2.667642e-07
> sin(1)^2 + cos(1)^2
[1] 1
> exp(800) [1] Inf > exp(-700) [1] 9.859677e-305
計算式はExcelと同じように、１列に入力する。かけ算は「*」、割り算は「/」、べき乗は「^」、通常の優先順位に従って左から右に計算される。優先順位を乱すためのカッコは「( )」しか使えない。pi は円周率を表す定数である。
計算結果は、小数点以下がなければ整数として、さもなければ小数点付きで、数値の範囲によっては有効桁に10のべき乗を掛けた指数形式で表示される。有効桁の７桁はデフォルト値で、options(digits=16) などで指定できる。あまり大きな数は Inf という記号で表示される（無限大 INFinity のつもり）。
関数電卓にあるような数学関数は全て使える。たとえば、三角関数は「sin(), cos()」などと書く、というような約束事を覚えれば、高性能の電卓のように、あるいはExcelのように便利な道具として使うことが出来る。
スペースは 無視されるので、数式が見やすくなるように、適宜挿入してよい。
完全な数式を入力したつもりで結果が表示されない場合、エラーメッセージが表示されるか、「+」が表示されて入力待ちになる。後者の場合は、カッコが閉じられていない、などの入力ミスが考えられるので、直前の入力行をチェックして（直前の行の続きを）追加入力し、命令を完成させてEnterキーを押す。最初から入力し直したい場合はEscキーを押すと、直前の入力はキャンセルされて「>」表示に変わる。
エラーメッセージが表示された場合は、入力した命令が間違っているので、エラーメッセージに従って修正してから入力し直す。「↑」キーを押すたびに、すでに入力した式が行単位で次から次へと表示されるので、もし一度入力したものを小修正すればよいのであれば、↑キーで目的の行を表示させ、「←」「→」キーを使って修正箇所までカーソルを移動させ、修正してからEnterキーを押すことで入力が完了する。
単純な数式ならば、コンソール画面でリアルタイムの作業しても構わないが、複雑な作業はエディタを使って数式を完成させ、１行ごとに console に送って実行させるというバッチ作業方式の方が効率的である。
RStudio 画面の左上がエディタ（Rソース入力）である。（RStudio を使わない場合）R の場合は、R Gui 画面で「ファイル」メニューから「新しいスクリプト」サブメニューを選択すると、Rエディタの「無題」という名前のウィンドウが開き、入力可能になる。
エディタで数式、あるいは命令文を入力し、入力行にカーソルを置いて F5 キーを押す（Ctrlキーを押しながら r キーを押しても同じ、以後 Ctrl+r キーと記す）と、その行が R console 画面に送られて実行される。mac の場合は ?+return キーを押す。
複数行を入力し、必要な行をすべて選択し、Ctrl+rキーを押せば、入力したものが1行ずつ R console 画面に送られて順次実行される。
エディタ画面は「file/save」メニューで保存できる。「file/new」メニューで別のエディタ画面を生成できる。 RStudio の場合はタブで画面を切り替えることができる。また、「file/close」メニューを実行しなければ、再度立ち上げた時、前回の作業文書も表示され、直ちに作業を再開できる状態になっている。
R の場合は、「ファイル」メニューから「保存」サブメニューを選び、適当なファイル名を付けて保存すると、拡張子「.R」が付加されて保存される。再度立ち上げた時は、「ファイル」メニューの「ファイルを開く」サブメニューを選ぶことによって、呼び出さなければいけない。
エディタなので、計算内容やアルゴリズムなどを日本語で書いておけば、文書管理としても有効である。
R は統計データを処理するために作られたソフトなので、データの集合としてのベクトルを対象とした計算が効率的に実行できるように工夫されている。統計処理でなくても、この処理方式を取り入れれば、計算を効率化できる。例えば、利率 r=0.05 で n=100 年間運用した場合の元利合計（複利計算）(1+r)n は 1.05^100 と入力すれば良いが、いろいろな年数、たとえば n=1,2,3,5,10,20,50,100 について計算したい場合、８通りの同じような数式を入力する代わりに、べき乗の n をベクトルにして、次のように書けば８通りの計算を一つの式で代用することができる。
> 1.05^100
[1] 131.501258 > 1.05^c(1,2,3,5,10,20,50,100) [1] 1.050000 1.102500 1.157625 1.276282 1.628895 2.653298 11.467400 [8] 131.501258
c(... ) は、sin, cos などの数学関数と同じように、c という名前の関数であることを表し、その引数を一つのベクトルにまとめて返す、と決められている。べき乗の肩にベクトルが書かれると、ベクトルの各要素にべき乗の計算を施し、結果をベクトルで返す、と決められている。結果の各行の先頭にはこのように鉤括弧で囲まれた整数が表示されるが、これは結果のベクトルの要素番号を示す。ただし、改行された先頭の要素についてだけ表示される。この場合、c(...) で定義されたベクトルの要素数は８なので、結果のベクトルも８個の要素からなり、２行目の先頭にある [8] は、８番目の要素は 131.5,,, であることを示す。n=100 の場合だけを計算する場合、計算結果はスカラー量であるが、左端にベクトルの要素番号が表示されている。これは、R がスカラーを長さ１のベクトルとみなして、要素番号を付けて表示しているのである。
いろいろな利率について計算したい場合も同様に、r をベクトルにして実行させると、年数を変えた場合と同様に、利率を変えた様々な元利合計がベクトルとして表示される。
> (1 + c(0.01, 0.02, 0.05, 0.1))^10 [1] 1.104622 1.218994 1.628895 2.593742
いろいろな利率を色々な年数について計算したい場合は、２次元の行列になるので、外積関数 outer を利用する。
> (r = c(0.01, 0.02, 0.05, 0.1))
[1] 0.01 0.02 0.05 0.10
> (n = c(1,2,3,5,10,20,50,100))
[1] 1 2 3 5 10 20 50 100
> outer(r, n, function(r,n) (1+r)^n)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1.01 1.0201 1.030301 1.051010 1.104622 1.220190 1.644632 2.704814
[2,] 1.02 1.0404 1.061208 1.104081 1.218994 1.485947 2.691588 7.244646
[3,] 1.05 1.1025 1.157625 1.276282 1.628895 2.653298 11.467400 131.501258
[4,] 1.10 1.2100 1.331000 1.610510 2.593742 6.727500 117.390853 13780.612340
結果の数値は４行８列の行列の形をしており、左端の鉤括弧内数字は行番号、上端の鉤括弧内数字は列番号を表している。r 行 n 列に (1+r)n の値が表示されている。
自分で行列を定義する場合は matrix 関数を使う。ベクトルと行数を引数とすると、ベクトルの要素を行優先で、与えられた行数を持つ行列の形に順番に流し込む。
> matrix(c(11,12,21,22,31,32), 2)
[,1] [,2] [,3]
[1,] 11 21 31
[2,] 12 22 32
複雑な計算を実行するために、計算途中の結果を、名前を付けた変数に一時的に記憶させておくことができる。変数に記憶させる命令は代入文と呼ばれ、「変数 = 計算式」あるいは「変数 <- 計算式」と書く。たとえば、
> x = 3 > y = 4 > (z = sqrt(x^2 + y^2)) [1] 5 > s <- (x + y + z)/2 > sqrt(s * (s-x) * (s-y) * (s-z)) [1] 6
等号記号は左右両辺が等しいという意味ではない。計算結果を変数に送る、というつもりで矢印に見える記号「<-」を代入記号とする流儀もあるが、慣れてしまえば、キーストロークが少ない方が便利。代入文の左辺は必ず一つの変数でなければいけない。
>」が表示され、結果が表示されない。代入文全体を「()」で 囲むと代入された変数の値が表示されるようになる。
計算式の結果がベクトルあるいは行列になる場合は、代入された変数はそのベクトルあるいは行列全体の名前になり、それに鉤括弧を添えたものが要素の名前になる。
> (r = c(0.01, 0.02, 0.05, 0.1)) [1] 0.01 0.02 0.05 0.10
> r[4]
[1] 0.1
> (A = matrix(c(11,21,12,22,13,23), 2)) [,1] [,2] [,3] [1,] 11 12 13 [2,] 21 22 23 > A[2,1]
[1] 21
最後の例の表示は奇妙に見えるかもしれない。A[2,1] はベクトルではなく、一つの数だが、R はそれを要素数１のベクトルとみなして、先頭に鉤括弧付き数字を表示しているのである。
R の計算処理は、四則演算などの初等的な数式計算は別として、普通は「関数」を呼び出すことによって実行される。関数と言っても、R における関数は、数学用語の関数という役割を持つものもあるが、プログラミング言語における「手続き」と同じで、一般に「まとまった処理（計算、作業）を実行する単位」を表す。
sqrt(x), abs(x), sin(x), exp(x) などは、数学関数の例である。ただし、独立変数 x はベクトルでもよく、その場合は、関数記号がベクトルの各要素に適用された結果をベクトルとして返す。次の例で、2:10 は公差１の等差数列を表す。
> sqrt(2) [1] 1.414214
> sqrt(2:10) [1] 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751 2.828427 3.000000 [9] 3.162278
max(x), mean(x) は、ベクトルを引数として一つの値を返す関数の例である。
sort(x) は引数のベクトルを大きさの順に並べ替えたベクトルを返す関数の例である。
matrix(...) は引数のベクトルを行列に変える関数、solve は引数の行列をその逆行列という別の行列に変える関数、eigen は引数の行列の固有値と固有ベクトルを返す関数である。
lm は、引数のデータを用いて線形回帰モデル分析を実行し、その結果を返す関数である。同様に、aov は分散分析を実行する。
optimize は最適化計算を実行する関数である。
plot は平面曲線を描く関数である。
R を賢く使うためには、分析目的に合う適切な関数を見つけて利用することが肝心である。
curve 関数を使って、関数のグラフを描くことができる。定義域を均等に分割した n 点で関数値を計算し、折れ線で結んだもので近似される。デフォルトは n=101 であるが、パラメータで指定できる。次の例で、dnorm は正規分布の密度関数の名称である。abline 関数を使って、現在表示されている図形のスクリーンに軸や通常の直線を追加することができる。h=0 は原点を通る垂直線、v=0 は原点を通る水平線を描く。
> curve(dnorm, -3, 3) > abline(h=0) > abline(v=0)
例えば次のように描かれる。
add=TRUE オプションを使うと、複数のグラフを同じ画面に描画して比較することができる。lty= は線の種類を指定するオプションである。legend 関数は凡例を作成して表示する。
> curve(dt(x, 2), add=TRUE, lty="dotted")
> curve(dt(x, 10), add=TRUE, lty="dashed", col=2)
> legend(2, 0.35, c("norm", "t(2)", "t(10)"), lty=c(1,3,2), cex=0.8, y.intersp = 1.5)
例えば次のように描かれる。
任意の平面曲線は、分点を十分に細かく取ってその座標を計算し、それらの点を折れ線で結ぶことで、近似的に描くことができる。次の例で、type="l" オプションは折れ線を描くための指定である。type= オプションのデフォルトは計算された座標点をプロットするだけである。points 関数は、現在描かれているスクリーンに点列を追加表示する。rnorm(n) は標準正規乱数を n 個生成する。
# 曲線を描く
> t = seq(0, 2*pi+0.1, by=0.1) > plot(4 * cos(t), y = 2 * sin(t), type="l")
# 散布図を描く
> n = 1000 > points(rnorm(n), rnorm(n)/2)
例えば次のように描かれる。
persp 関数は２変数関数を透視図法で描画する。次の例は２変量正規分布の３Dグラフを描いたものである。outer(x, y, f) は、外積を計算する関数で、x, y のすべての組み合わせつについて２変数関数 f の値を計算して行列として返す。persp は、計算された曲面上の点をワイヤーフレームで表示する関数である。theta=, phi= オプションは視点を指定する。
# ２変量正規分布の密度関数の3D画像
> f = function(x,y,r=0.6) exp(-(x^2-2*r*x*y+y^2)/2/(1-r^2)) / sqrt(1-r^2) / 2 / pi
> x = y = seq(-2, 2, length=50) > z = outer(x, y, f) > persp(x, y, z, theta=120, phi=20, expand=0.5, ticktype="detailed")
例えば次のように描かれる。
contour 関数を使うと、等高線表示ができる。asp=1 は縦軸と横軸の単位長さを等しくするオプションである。image 関数は等高線の代わりに違う高さには違う色を使った塗り絵を表示する。
# ２変量正規分布の密度関数の等高線と色分け図 > contour(x, y, z, asp=1)
> image(x, y, z, asp=1)
２変量正規分布の密度関数に対して、左が contour 関数、右が image 関数を使った描画例である。
Rは統計解析のためのパッケージなので、データ集計や、統計解析の結果を視覚的に表現する様々なグラフが利用可能である。barplot は棒グラフ、hist はヒストグラム、piechart はエングラフ、boxplot は箱ヒゲ図、などなど。
Rは様々な関数を組み合わせて目的の分析を実行する。パブリックドメインのソフトなので、探せば誰かが目的にあった関数を作って公開している可能性がある。もしなければ自分で作れば良い。
簡単な数式評価だけではなく、条件分岐や繰り返しを含むような複雑な計算も、if else 構文や、while, repeat, for 構文のようなプログラミング機能を使えば実現することができる。
例えば、次はニュートン法による方程式の解を計算するものである。function(x) { ... } は関数を定義する。f, df は目的の関数と導関数を定義している。repeat { ... } は、かっこに挟まれた命令文群を繰り返すという命令、if( ... ) break は条件が満たされたら repeat ループを脱出するという命令である。print 関数は変数を表示する。
> f = function(x) x * (x^2 - 1)
> df = function(x) 3*x^2 - 1
> x = -0.2
> epsilon = 1e-10
> repeat {
+ xa = x
+ x = xa - f(x) / df(x)
+ print(x)
+ if(abs(x - xa) < epsilon) break
+ }
[1] 0.01818182
[1] -1.203297e-05
[1] 3.484565e-15
[1] 0
１次元のベクトルを決められた個数ごとに区切って２次元の表に並べかえたものを行列という。matrix 関数は、ベクトルを行列に置き換える（といっても新たな実態を作るわけではなく、２つの添え字で参照できるようにするだけ）。区切ったものを列として並べるか（これがデフォルト）、行として並べるか、指定することができる。同じサイズのベクトルを並べて行列を構成することもできる。rbind 関数は与えられたベクトルたちを行ベクトルとして行列化し、cbind 関数は列ベクトルとして行列化する。rbind 関数は、行列に新たな行を追加する場合にも使用される。diag 関数は単位行列や対角行列を生成する。
> (A = matrix(1:4, 2))
[,1] [,2]
[1,] 1 3
[2,] 2 4
> rbind(c(1,2,3,4), c(5,6,7,8))
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
> rbind(A, c(3,5))
[,1] [,2]
[1,] 1 3
[2,] 2 4
[3,] 3 5
> cbind(c(-1,0), A)
[,1] [,2] [,3]
[1,] -1 1 3
[2,] 0 2 4
> diag(2)
[,1] [,2]
[1,] 1 0
[2,] 0 1
二つの行列 A, B の行列としての積は A %*% B とする。A * B とした場合は要素ごとに掛け算が実行される。したがって、後者の場合は、行列のサイズが同じでなければいけない。同じように、* の代わりに +, -, / とした場合も要素ごとに演算が実行される。正方行列 A の逆行列は solve(A) で計算できる。連立一次方程式 A x = b は solve(A) にベクトル b を掛けたものが解ベクトルになるが、solve(A, b) で求めることができる。
> A = matrix(1:4, 2)
> solve(A)
[,1] [,2]
[1,] -2 1.5
[2,] 1 -0.5
> A %*% solve(A)
[,1] [,2]
[1,] 1 0
[2,] 0 1
> b = c(5:6)
> solve(A, b)
[1] -1 2
> solve(A) %*% b
[,1]
[1,] -1
[2,] 2
> A %*% solve(A, b)
[,1]
[1,] 5
[2,] 6
det(A) は行列 A の行列式、eigen(A) は行列 A の固有値と固有ベクトルを計算する。eigen(A) の出力は$value によって固有値を、$vector によって固有ベクトルを取り出すことができる。
> det(matrix(1:4, 2))
[1] -2
>
> (v = eigen(A)) $values [1] 5.3722813 -0.3722813 $vectors [,1] [,2] [1,] -0.5657675 -0.9093767 [2,] -0.8245648 0.4159736
> A %*% v$vector[,1] - v$value[1] * v$vector[,1] [,1] [1,] 0.000000e+00 [2,] -8.881784e-16
例（統計数値表） qnorm 関数は正規分布の累積分布関数の逆関数である。qt(x, df) は自由度 df の t 分布の累積分布関数の逆関数である。これらを使って、正規分布、t分布の統計数値表を作ることができる。colnames 関数は列に名前を付け、rownames 関数は行に名前を付ける。
> a = c(0.8, 0.9, 0.95, 0.975, 0.99, 0.995)
> tt = rbind(qt(a, 3), qt(a, 10), qt(a, 20), qt(a, 60), qnorm(a))
> colnames(tt) = a
> rownames(tt) = c("t_3", "t_10", "t_20", "t_60", "norm")
> tt
0.8 0.9 0.95 0.975 0.99 0.995
t_3 0.9784723 1.637744 2.353363 3.182446 4.540703 5.840909
t_10 0.8790578 1.372184 1.812461 2.228139 2.763769 3.169273
t_20 0.8599644 1.325341 1.724718 2.085963 2.527977 2.845340
t_60 0.8476530 1.295821 1.670649 2.000298 2.390119 2.660283
norm 0.8416212 1.281552 1.644854 1.959964 2.326348 2.575829
１と自分自身以外に約数を持たない正整数を素数という。ある数が素数かどうかをチェックするには、それより小さい数で割った余りがすべて０でないことをチェックする必要がある。チェックする数の最大値はその数の平方根である。順番にチェックする代わりに、ベクトル計算を使えば、一つの式で実現できる。次の例で、%% は余りを計算する演算子だが、チェックすべき数すべてについてあまり計算を実施し、それが０でなければ TRUE さもなければ FALSE という論理値に置き換え、その論理値がすべて TRUE ならばその数は素数であると言える。all 関数は論理値ベクトルを引数とし、それらがすべて TRUE であるとき TRUE を返す。
ある範囲の数を同時にチェックして素数表を作ることができるアルゴリズムとして、エラトステネスの篩（ふるい）がある。数の集合に対して、２で割り切れるものをふるい落とし、残ったものの中から３で割り切れるものをふるい落とし、残ったものの中から５で割ったものをふるい落とし、... という作業を続け、最後に残されたものが素数である、という方法である。これはプログラミングの問題としてよく取り上げられるが、普通に思いつくのは for 構文のプログラムである。for 構文と同じような振る舞いをする関数に sapply 関数がある。これを使えば、関数だけで篩を実行することができる。例えば次のように。
# 素数の判定
> m = 2^31-1 > if(all(m %% 2:sqrt(m) != 0)) "素数" else "合成数" [1] "素数"
# エラトステネスの篩
> n = 100 > prime = 1:n > prime[seq(4,n,2)] = 0 > prime[1] = 0 > z = sapply(seq(3,sqrt(n),2), function(p) prime[seq(2*p,n,p)] = 0) > prime[prime != 0] [1] 2 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 [27] 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99
## sapply 関数は次と同じ
> for(p in seq(3,sqrt(n),2)) prime[seq(2*p,n,p)] = 0
整数の範囲で解を見つける問題には、見かけは単純で解決のつかない問題が多数ある。コラッツの問題もその一つである。
正整数を一つ決め、偶数ならば２で割り、奇数ならば３倍して１を足す、という計算を繰り返すと、そのうちに「必ず」１に到達することを証明せよ、という問題。
どのようにして１に到達するのか確かめようとすれば、地道に repeat 構造を使ってプログラムを書くしかない。まとめて計算する場合は、上と同じように sapply 関数や lapply 関数を使うことができる。
> collatz = function(n = 100) {
+ w = n
+ repeat{
+ if(n%%2 == 0) {
+ n = n/2
+ } else n = 3 * n + 1
+ w = c(w, n)
+ if(n == 1) break
+ }
+ return(w)
+ }
> collatz()
[1] 100 50 25 76 38 19 58 29 88 44 22 11 34 17 52 26 13 40 20
[20] 10 5 16 8 4 2 1
> collatz(201)
[1] 201 604 302 151 454 227 682 341 1024 512 256 128 64 32 16
[16] 8 4 2 1
# 大きな数で（まとめて）確かめてみる
> x = 100000000
> z = sapply(x:(x+10000), collatz)
# １に到達するまでの長さに特徴があるか調べてみる
> table(unlist(lapply(z, length)))
82 108 139 157 162 170 188 214 276 307 320 325 338 351 369 413 444 488
771 1715 1608 1682 865 226 1254 661 282 592 117 23 8 33 74 2 86 2
連続関数で区間の端点で関数値の符号が異なるとき、uniroot 関数はその関数のゼロ点を返す。
> uniroot(cos, c(1,3)) $root [1] 1.570813 $f.root [1] -1.688029e-05 $iter [1] 4
関数名を引数とするが、関数定義を書き込んでもよい。区間の端点で同符号となる場合はエラー表示される。収束の途中経過を表示させることもできる。
> uniroot(function(x) x^3 - 3*x + 1, c(1,10)) $root [1] 1.532084 $f.root [1] -2.134805e-05 $iter [1] 13
> uniroot(function(x) x^3 - 3*x + 1, c(2,10)) 以下にエラー uniroot(function(x) x^3 - 3 * x + 1, c(2, 10)) : f() の端点での値が異なった符号を持ちません
> uniroot(function(x) print(x)^3 - 3*x + 1, c(-10,10)) [1] -10 [1] 10 [1] -0.01030928 [1] -0.02092609 [1] -5.010463 [1] -0.06877696 以下略
n 次多項式関数のゼロ点を求めたい場合は、polyroot 関数を使えば、複素数解を含めて n 個の解を同時に求めることができる。
# 実数解だけの場合： 1 - 3x + x^3 = 0
> polyroot(c(1,-3,0,1)) [1] 0.3472964-0i -1.8793852+0i 1.5320889+0i
# 複素数解を含む場合： 1 - x + x^3 = 0 > polyroot(c(1,-1,0,1)) [1] 0.662359+0.5622795i -1.324718-0.0000000i 0.662359-0.5622795i
線形連立方程式は、solve 関数を使って解くことができる。
> (A = matrix(1:4, 2))
[,1] [,2]
[1,] 1 3
[2,] 2 4
> b = c(5:6)
> solve(A, b)
[1] -1 2
D 関数は、導関数を式として出力する。関数は expression 関数の引数として与える。本格的な数式処理は期待しない方が良い。
> D(expression(x^2), "x") 2 * x > D(expression(1 / sqrt(1 - x^2)), "x") 0.5 * (2 * x * (1 - x^2)^-0.5)/sqrt(1 - x^2)^2
deriv 関数は数値微分を計算する関数を定義する。定義された関数は、関数値と数値微分の値を返す。数値微分の値だけを取り出す場合は attr(, "gradient") とする。
> (df = deriv(~ (exp(-x^2)), "x", func = T))
function (x)
{
.expr3 <- exp(-x^2)
.value <- .expr3
.grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
.grad[, "x"] <- -(.expr3 * (2 * x))
attr(.value, "gradient") <- .grad
.value
}
> df(0:3)
[1] 1.0000000000 0.3678794412 0.0183156389 0.0001234098
attr(,"gradient")
x
[1,] 0.0000000000
[2,] -0.7357588823
[3,] -0.0732625556
[4,] -0.0007404588
# 導関数のグラフを描く
> xx = seq(-3,3,0.1)
> plot(xx, attr(df(xx),"gradient"), xlim=c(-3,3), type="l")
integrate 関数は定積分を計算する。被積分関数は、その場で定義しても良いし、関数定義を書いてその関数名を引用しても良い。
Inf を無限大を表す記号として、無限区間の広義積分を計算することもできる。また、関数値が端点で発散する場合は、開区間での積分を計算する。
> integrate(function(x) 1-x^2, -1, 1) 1.333333 with absolute error < 1.5e-14
> integrate(function(x) 1/sqrt(x), 0, 1) 2 with absolute error < 5.8e-15
> integrate(function(x) 1/x^2, 1, Inf) 1 with absolute error < 1.1e-14
> integrate(function(x) 1/sqrt(x), 0, Inf) 以下にエラー integrate(function(x) 1/sqrt(x), 0, Inf) : the integral is probably divergent
シミュレーションのための、各種確率分布に従う乱数を生成することができる。一様乱数は runif 関数、正規乱数は rnorm 関数など、分布の名前の前に r をつけたものが乱数生成用の関数である。
> rnorm(5) [1] 1.1828169 0.4154934 -0.2823300 -0.9743875 -1.6157993
sample 関数は、一般の離散分布に従う乱数を生成する。sample 関数のデフォルト値は非復元抽出である。抽出個数を指定しないとランダム置換になる。replace=TRUE オプションで、復元抽出が実行できる。次の例で、sort 関数は配列を昇順に並べ替える。apply 関数は行列の行、あるいは列に同じ操作を適用 apply する。
# さいころ投げのつもり
> sample(1:6, 20, replace=TRUE) [1] 4 4 4 4 3 3 5 5 5 3 6 1 2 1 2 2 5 1 5 1 > table(sample(1:6, 6000, replace=TRUE)) 1 2 3 4 5 6 999 980 1114 932 975 1000
# ビンゴボードのつもり
> (bingo = matrix(sort(sample(1:75, 25)), 5)) [,1] [,2] [,3] [,4] [,5] [1,] 1 12 30 39 51 [2,] 2 17 31 41 56 [3,] 6 23 32 42 71 [4,] 7 24 36 45 74 [5,] 8 29 38 50 75
# あるいは、もうひとひねり > apply(bingo, 2, function(w) w[sample(1:5)]) [,1] [,2] [,3] [,4] [,5] [1,] 6 12 31 39 51 [2,] 7 23 30 45 56 [3,] 8 29 32 50 75 [4,] 2 24 36 42 71 [5,] 1 17 38 41 74
モンテカルロ法は、決定論的問題を確率変数の期待値として定式化し、乱数実験でその期待値を推定する方法である。複雑な問題で簡単には解けない問題でも、推定誤差という問題はあってもある程度の精度で解を求められることから、いろいろな問題で用いられている強力な問題解決の道具である。
例えば、円周率πは、区間 [0, 1] で定義された独立な確率変数 X, Y に対して、X2 + Y2 < 1 ならば１さもなければ０という値をとる確率変数の期待値の４倍として定式化できる。一様乱数は runif 関数で生成し、条件判定の結果はそのまま 0-1 確率変数の値となる（TRUE は 1、 FALSE は０）ことを利用して、mean 関数で標本平均を計算する。統計理論に従って 95% 信頼区間を計算した。
> n = 10000
> N = 1000
> z = sapply(rep(n,N), function(n) 4 * mean(runif(n)^2 + runif(n)^2 < 1))
> cat("[", mean(z)-2*sd(z)/sqrt(N), ",", mean(z)+2*sd(z)/sqrt(N), "]\n")
[ 3.140048 , 3.142117 ]
もう少し複雑な例で、同じく円周率を推定する「ビュフォンの針の実験」について、説明する。床に等間隔に平行線を引き、たくさんの針を投げて、それらの平行線のどれかと交わった針の本数を数えると円周率が推定できる、というのがビュフォンの針の実験である。
L を平行線の間隔、b を針の長さとする。平行線と交わるかどうかだけを考えるのであれば、２本の平行線とそれに垂直に交わる直線を考え、針の中心が垂直線の平行線に挟まれる部分上にランダムに落ちると考えれば良い。
針の向きは 0からπまでランダムに選ばれるとすると、針が平行線と交わるということは、中心より上の部分が上の平行線と交わるか、中心より下の部分が下の平行線と交わるかのいずれか、と考えると次ようなプログラムができる。
runif は一様乱数を生成する関数なので、t がランダムな針の向き、y がランダムに置かれた針の中心座標を与える。その結果、y+yn が針の上端、y-yn が針の下端になり、針が平行線と交わるためには、上端が L を超えるか、下端が 0 より小さくなるかすれば良い。mean の中身は針が平行線と交われば TRUE、さもなければ FALSE という論理値になるが、mean は TRUE を１、FALSE を０として平均を計算するので、条件を満たす針の比率が計算される。
plot 関数の type="n" オプションは、領域だけ確保して、プロット領域には何も描かない、というオプションである。lines で n 本の針を描き、points で針の中心を描いている。
> L = 0.25 > b = 0.2 > n = 200 > t = runif(n, 0, pi) > y = runif(n)
> yn = b/2*sin(t) > (est = 2 * b / L / mean(y%%L + yn > L | y%%L - yn < 0)) [1] 3.106796
# グラフィックデモ用
> x = runif(n); xn = b/2*cos(t) > plot(c(0,1), c(0,1), type="n") > abline(h=L*(0:4)) > for(i in 1:n) lines(c(x[i]+xn[i],x[i]-xn[i]),c(y[i]+yn[i],y[i]-yn[i])) > points(x,y)
次は実行例である。
統計データの基本統計量を計算する関数には次のようなものがある。
| 関数名 | 説明 | 関数名 | 説明 | 関数名 | 説明 |
|---|---|---|---|---|---|
mean |
平均 | var |
不偏分散 | sd |
標準偏差 |
max |
最大値 | min |
最小値 | sum |
総和 |
range |
範囲 | median |
中央値 | quantile |
分位点 |
cor |
相関係数 | cov |
共分散 | summary |
四分位数、平均値 |
大量のデータを入力する場合は、Excel やエディタのような道具を使ってデータを整え、それらを一挙に取り込むようにすると効率が良い。コピーペーストで扱える程度の量であれば、scan 関数を利用する。キーボードから入力する場合、データの区切りはスペースが標準（sep="," でカンマ区切りでも読める）、Excel の行、あるいは列の一部をコピーして、入力プロンプトの後にペーストしても良い。ペーストの後、次の入力プロンプトが表示されるので Enter キーを押してスキャンモードを終了させる。下の例では「4:」の後に空行が入力されている。
> z = scan()
1: 1 3 0.01 4: Read 3 items > z [1] 1.00 3.00 0.01
統計データのような２次元の表形式のデータは、Excelシートで整形してからRに取り込むとよい。表の部分をコピーし、read.table 関数を使ってクリップボード経由で取り込むか、カンマ区切りのCSVファイル形式で保存して、read.csv 関数を利用する。
ファイルは MyDocument ファイルに保存してあればそのまま読めるが、そうでない場合は、完全なファイルパス名を入力するか、Rconsole 画面で「ファイル／ディレクトリの変更...」メニューを選び、保存したデータファイル の入っているディレクトリを指定する必要がある。現在のワーキングディレクトリは getwd 関数で知ることができる。
# クリップボード経由（Excelデータをコピーした後実行する）
df = read.table("clipboard", header=F) # windows の場合
df = read.delim(pipe("pbpaste")) # mac の場合
# csv 形式で保存されたファイルから読み込む場合
df = read.csv("stocktestdata.csv", header=T, skip=2)
計算結果をファイルに保存したい場合は write 関数を使う。データはタブ区切りになる。ファイルはパス名を書かない場合、現在のワーキングディレクトリに保存される。「.txt」拡張子をつけておけば、内容をエディタで見ることができる。
write(t(A),"data.txt",ncol=ncol(A)) # 行列Aの転置（行優先）を書き出す
sink 関数は、Rconsoleへの計算結果をファイルに保存することを指示する。
sink("result.txt") # ここに書かれた命令の計算結果がすべてファイル保存される
sink() # 何も指定しなければ「マイドキュメント」に保存される
もどる