コンソール画面で「>」が表示されている状態から入力を始める。> はプロンプト文字と呼ばれる。
数式、もしくは命令文を入力して Enter キーを押すと、Ｒは入力内容の処理を始める。Enter キーが入力されるまでは、←、→キーを使って、その行の入力した内容を修正することができる。
入力内容が完全ならば、直ちに結果が表示される。もし、左端に「+」が表示されて入力待ちになった場合は、直前の入力が完全な数式、もしくは命令文になっていない、ということを意味する。続きがあれば入力して Enter キーを押すか、さもなければ Esc キーを押してそれまでの入力をキャンセルする。その結果、プロンプト文字「>」の表示状態に戻る。
> x <- 3; y <- 4 # 「;」で区切って複数の命令文を入力できる > sqrt(1/x^2 + y^2 + ) # 最後に付け加えれば数式としては完成する
[1] 4.013865 > sqrt(1/(x^2 + y^2)) # ↑キーを1回押して、直前の入力行を修正する（本来はこちらが正しい） [1] 0.2
過去に入力したものを再入力したい場合は、上向き矢印キーを何回か押して、目的の入力行が表示されたところで Enter キーを押す。
過去に入力したものの一部を書き換えて再入力したい場合も同様に、目的の入力行を表示させ、左向き矢印キーを使って修正箇所までカーソルを移動させ、修正して から Enterキーを押す。
コンソールでの入力は、簡単な数式ならば良いが、かっこが何重にもなっていたり、少しずつ異なる数式を何行も入力するような場合は、煩わしく感じられる。エディタを使って入力する式をあらかじめ整えておいて、１行ずつまとめて入力できるようになれば、デバッグなどの作業が楽になる。
RStudio ならば左上にエディタが表示されているが、R Gui 画面ならば「ファイル」メニューから「新しいスクリプト」サブメニューを選択することで、Rエディタの「無題」という名前のエディタが開き、入力可能になる。
エディタで R で入力する関数や命令文を（何行でも）入力したあと、ある行にカーソルを置き、F5 キー（mac では?キーを押しながら Enterキー）を押すと、その行が R console 画面に送られて実行される。複数行を選択して F5 キー（?+Enter）を押すと入力したものが1行ずつ R console 画面に送られて順次実行される。（windows では F5 キーの代わりに Ctrlキーを押しながら r キーを押しても同じ、以後 Ctrl+r キーと記す）
「ファイル」メニューから「保存」サブメニューを選び、適当なファイル名を付けて保存すると、拡張子「.R」が付加されて保存 され る。次回、R Gui 画面で、「ファイル」メニューの「ファイルを開く」サブメニューを選ぶことによって、呼び出すことが出来るようになる。RStudio ならば、一々ファイルを閉じる、開くの操作をしなくても、次回立ち上げた時に前回の作業状態が復元される。エディタなので、計算内容やアルゴリズムなどを日本語で書いておけば、文書管理としても有効である。
Ｒの扱うデータには、数、文字列、論理定数、などがあり、それらは別々の型(type)を持つ。 代表的な型には numeric(integer, double), complex, character, logical がある。
integer, double はコンピュータの内部表現の違いで、実際の計算では両方を合わせて numeric 型として扱われる。
計算では、それぞれの場合に適切なデータの型を選ばないと、不適切な結果となるか、あるいは場合によってはエラーになるので注意が必要である。
異なるデータ型が一つの式に混在している場合は次の順番に従って、型を自動変換する。
logical → numeric → numeric → character
ある変数がどの型として扱われているかをチェックするために is.xxx という関数が使われる。また、ある変数を強制的にある型の変数とするために as.xxx という関数が使われる。上の例で何が起きたのかをチェックする。rm 関数は引数の変数を未定義状態にする（それまでのユーザ定義を無効にする）。ここでは rm(T) は T を TRUE の省略形として使えるようにするという意味を持つ。
> T = 0 > is.logical(T) [1] FALSE > rm(T) > is.logical(T) [1] TRUE
> as.complex(1) [1] 1+0i
四則演算は「+, -, * /」を使う。「-」は符号を変える単項演算子としても用いられる。
べき乗計算は「^」を使う。
整数としての割り算（整数同士でなくてもよい）の商は「%/%」、あまりは「%%」 を使う。
計算の順番を変えるために使うカッコは「( )」のみ許される。
式の途中にある空白文字は無視される。
> 2.5 %/% 1.2 [1] 2 > 2.5 %% 1.2 [1] 0.1
pi は円周率定数として定義されている。options(digits=xx) は表示桁数を指定する。大きな数字を設定しても、実際の有効桁は 16 桁ていど。
大きな数は 10 のべき乗に１以上10未満の数をかけた指数形式で表示される。演算結果が表現しきれない大きな数（21024以上、10進数で307桁程度）になった場合は Inf という記号が返される。
> options(digits=30) 以下にエラー options(digits = 30) : 'digits' のパラメータが不正です: 0...22 が認められています > 2^83 > pi # 本当は pi = 3.141592653589793238... [1] 3.141592653589793116
> for(i in 1:18) print(10^i + 10^(-i) + i) [1] 11.099999999999999645 [1] 102.01000000000000512 [1] 1003.0009999999999764 [1] 10004.000099999999293 [1] 100005.00001000000339 [1] 1000006.0000010000076 [1] 10000007.000000100583 [1] 100000008.0000000149 [1] 1000000009 [1] 10000000010 [1] 100000000011 [1] 1000000000012 [1] 10000000000013 [1] 100000000000014 [1] 1000000000000015 [1] 10000000000000016 [1] 100000000000000016 [1] 1e+18
計算結果を一時的に記憶し、その結果を別の数式に利用することが出来る。記憶する場所を変数という。代入文は、変数に記憶させるための命令で、代入演算子「<-」あるいは「=」の左辺に変数名、右辺に記憶させる数（を計算する為の数式など）を書く。等号記号は左辺と右辺が等しいことを表すものではない。例えば、x = x+1 は x に１を足したものを新たな x とする、という意味である。
代入文を実行した結果は左辺の変数に代入される（記憶される）だけで、コンソール画面には表示されない。代入と同時に画面にも表示させたい場合は、代入文全体を「( )」で囲む。
異なる数式、もしくは命令文を「;」区切りで一行にして入力しても良い。
データの基本はベクトルである。スカラーも長さ１のベクトルとして扱われる。行列はそれらを２次元構造に見立てたもので、データフレームは列毎にデータ型の違うデータを行列に仕立てたものである。さらに、３次元、４次元などの構造に見立てたものは配列と呼ばれる。
c 関数は、同じ型のデータを一つのベクトルにまとめる。ベクトルの要素は、先頭から１番目、２番目のように順番が付けられ、[1], [2], のようにかぎ括弧を付けた添え字を変数に添えた「添え字付き変数」で参照できる。添え字は０からではなく、１から始まる。
ベクトルを表示する際、各要素にその添え字を添えて表示してもらえれば分かりやすいが、一方で煩わしくもある。そこで、各行の先頭の要素だけその添え字を表示することにする。計算結果が表示されると、１行目の先頭に必ず [1] と表示されるのは、そのあと最初に表示されるものがベクトルの先頭要素（添え字１）であることを意味する。
次の例は 501 から 550 までの 50 個の整数をベクトルとして w という名前に記憶させている。501:550 は 501,502,...,550 を意味する省略記号である。先頭の [1], [24], [47] はそれぞれその行の先頭のデータ 501, 524, 547 の添え字がそれぞれ 1, 24, 47 であること、すなわち w[1] = 501, w[24] = 524, w[47] = 547 であることを示している。
例えば添え字付き変数 w[30] と入力するとその中身が表示されるが、その先頭に [1] と表示されるのは、入力した変数が新たな長さ１のベクトルを定義したことになり、その先頭要素（要素は一つしかないが）の添え字は [1] となることを示している（あまり意味はない）。
> (w = 501:550) [1] 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 [24] 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 [47] 547 548 549 550
> w[30] [1] 530
matrix 関数は同じ型のデータを行列にまとめる。ベクトル z と整数 n に対して matrix(z, n) とすると、z を n 個ずつ区切って列ベクトルとし、それらを左から順に並べて行列を作る。ベクトルの要素数を N とすると、出来上がる行列の大きさは n x (N/n) である。byrow=TRUE オプションを使うと、z を N/n 個ずつ区切って行ベクトルとし、それらを上から順に並べて行列を作る。
各要素は２次元の添え字付き変数（カンマ区切りで添え字２つを並べ、カギ括弧で囲んだものを変数の後に書いたもの）で参照できる。添え字は０からではなく、１から始まる。
cbind 関数はベクトルあるいは行列を横方向に並べて新たな行列を作る。rbind 関数は cbind 関数同様だが、縦方向に並べることが違う。
> z = 1:8
> matrix(z, 2)
[,1] [,2] [,3] [,4]
[1,] 1 3 5 7
[2,] 2 4 6 8
> (zz = matrix(z, 2, byrow=TRUE))
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
> zz[2,3]
[1] 7
> cbind(zz, c(0,1))
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 0
[2,] 5 6 7 8 1
> rbind(c(-3,-2,-1,0), zz)
[,1] [,2] [,3] [,4]
[1,] -3 -2 -1 0
[2,] 1 2 3 4
[3,] 5 6 7 8
> cbind(c(1:3), c(11:13), c(21:23), c(31:33))
[,1] [,2] [,3] [,4]
[1,] 1 11 21 31
[2,] 2 12 22 32
[3,] 3 13 23 33
data.frame 関数は、同じ型とは限らないベクトルを一つの行列（データフレームと呼ぶ）にまとめる。cbind 関数を使うと、行列は同じデータ型でなければいけないので、全て文字型になってしまう。> ### 異なるデータ型ベクトルを行列（表）にまとめる
> name = c("相田", "池田", "太田", "蒲田", "島田", "高田", "富田", "古田", "松田") > height = c(177.1, 171.7, 169.2, 173.4, 181.9, 171.9, 165.0, 175.0, 175.4) > weight = c(66.8, 72.6, 77.8, 75.7, 81.2, 70.3, 75.4, 81.5, 72.2) > personalData = data.frame(name, height, weight)
> personalData name height weight 1 相田 177.1 66.8 2 池田 171.7 72.6 3 太田 169.2 77.8 4 蒲田 173.4 75.7 5 島田 181.9 81.2 6 高田 171.9 70.3 7 富田 165.0 75.4 8 古田 175.0 81.5 9 松田 175.4 72.2
# 行列を作ると身長も体重も文字型データに変換される
> cbind(name, height, weight) name height weight [1,] "相田" "177.1" "66.8" [2,] "池田" "171.7" "72.6"
。。。以下省略
seq(x, y, by=a) は初項 x 公差 a の等差数列を y を超えない範囲まで生成する。
seq(x, y, length=n) は x と y の間を n-1 等分した等差数列を生成する。
公差が１あるいはマイナス１ならば x:y と書いても良い。
> seq(1, 3, by=0.3) [1] 1.0 1.3 1.6 1.9 2.2 2.5 2.8 > (0:6) * 0.3 + 1 [1] 1.0 1.3 1.6 1.9 2.2 2.5 2.8 > seq(1,3,length=9) [1] 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00
rep(a, n) は n 個の a を並べたベクトルを生成する。
numeric(n) は rep(0, n) と同じで、n 個の０を並べたベクトルを生成する。
数に関する条件判定式の中で、数を配列に置き換えると、配列の各要素ごとに条件判定を行われ、結果の論理値（TRUE, FALSE）を配列の要素数分並べたベクトルが生成される。
which 関数は、上のような配列に関する条件判定式を引数として、TRUE を持つ添え字集合を返す。
> 1:4 < pi [1] TRUE TRUE TRUE FALSE > which(11:21 %% 2 != 0) [1] 1 3 5 7 9 11
２番目の例は、11を１番目とする配列の中で奇数の位置（添え字）を取り出しているのでこのようになる。
diag(n)（n は整数）は n 次元単位行列を生成する。
diag(a)（a は要素数 n のベクトル）は対角要素が a の n 次元正方行列を生成する。
ベクトルの一部を取り出すには添え字付き変数の添え字の代わりに添え字集合を書けば良い。
which 関数を使って、ある条件を満たす要素だけを抽出することもできる。
添え字集合としてマイナスの添え字を書くと、その添え字を持つ要素を除外したベクトルを生成する。
ベクトルの一部の要素を書き換えるには、上の方法で部分抽出したものを代入文の左辺において代入操作を実行すれば良い。
同じ値で置き換えるためには replace 関数を使う。replace(z, A, a) は、配列 z の添え字集合 A に含まれる添え字を持つ要素をすべて a で置き換える。
c 関数の引数としてベクトルを指定すると、そのベクトルを含む新たなベクトルが生成される。
append 関数はあるベクトルの任意の場所に別のベクトルを挿入して新たなベクトルを生成する。
> x = 1:5 > y = 10:15 > append(x, y, after=3) [1] 1 2 3 10 11 12 13 14 15 4 5
> append(x, y) [1] 1 2 3 4 5 10 11 12 13 14 15 > c(x, y) [1] 1 2 3 4 5 10 11 12 13 14 15
sort 関数は、引数のベクトルを昇順に並べる。decreasing=TRUE オプションを使うと降順に並べることもできる。
rev 関数は、引数のベクトルの順番を逆順にする。
> (z = 37^(1:10) %% 100) [1] 37 69 53 61 57 9 33 21 77 49 > sort(z) [1] 9 21 33 37 49 53 57 61 69 77 > sort(z, decreasing=TRUE) [1] 77 69 61 57 53 49 37 33 21 9 > rev(sort(z)) [1] 77 69 61 57 53 49 37 33 21 9
平方根、指数、対数、三角関数、絶対値、四捨五入などの計算には関数を使う。Excelの関数とほぼ同じ感覚で使うことが出来る。
関数の引数（独立変数）としてベクトルを入力すると、ベクトルの各要素ごとに数学関数が適用され、結果は引数と同じ次元のベクトルになる。
> sqrt(2) [1] 1.414214 > sqrt(2:8) [1] 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751 2.828427
2:8 は 2 から 8 までの整数の集合を表す。したがって、この場合は √2、√3、...、√8 を計算している。
|
|
|
注意：trunc 関数は絶対値を取って切り捨ててから符号を戻すので、floor 関数とは負数の扱いが異なる。
> (x <- seq(-1.64, 2.3, by=0.41)) [1] -1.64 -1.23 -0.82 -0.41 0.00 0.41 0.82 1.23 1.64 2.05 > round(x) [1] -2 -1 -1 0 0 0 1 1 2 2 > ceiling(x) [1] -1 -1 0 0 0 1 1 2 2 3 > floor(x) [1] -2 -2 -1 -1 0 0 0 1 1 2 > trunc(x) [1] -1 -1 0 0 0 0 0 1 1 2 > sign(x) [1] -1 -1 -1 -1 0 1 1 1 1 1
R には、複雑な計算アルゴリズムにとって必須の条件分岐や繰り返しを実現する構文が用意されている。これらを組み合わせて与えれた問題を解く R の命令文群（プログラムという）を作ることをプログラミングという。
R はインタープリタ型なので、プログラムを実行させると、毎回１行ずつ機械語に翻訳して実行するので実行速度は遅い。そのため、プログラムは関数として定義し、あらかじめ機械語に翻訳させておき、必要に応じて関数を呼び出して計算させる、というのが R 流の使い方である。
ここではプログラミングと関数定義について説明する。
条件 A が満たされる場合だけ命令文群 P を実行したい場合は と書く。命令文群がただ一つの命令文からなる場合は
if(A) {P} {} を省くことができる。
条件 A が満たされる場合は命令文群 P を実行し、満たされない場合は命令文群 Q を実行したい場合は と書く。命令文群 Q が if から始まっても良いので、別の条件 B で分岐させることもできる：
if(A) {P} else {Q} if(A) {P} else {if(B) {Q} else {R}}
if ... else 構文で、複数行に分けて入力したい場合、else の直前で改行すると、それまでの入力文群が if 構文と解釈されてしまい、次は新規の命令文入力待ちとなる。else は命令文の先頭に来ることはないのでエラーとみなされるので注意が必要である。例えば、if(a) {P で改行し、直前の入力文がまだ入力途中であることをアピールする。
> a = 0 > if(a > 0) print(paste(a, "は正")) else print(paste(a, "は非正")) [1] "0 は非正"
> ## 正しくない例
> if(a > 0) print(paste(a, "は正")) > else print(paste(a, "は非正")) エラー: 予想外の 'else' です in "else"
> ## 正しい例
> if(a > 0) {print(paste(a, "は正")) + } else print(paste(a, "は非正")) [1] "0 は非正"
命令文群 P を繰り返し実行したい場合は と書く。
repeat{P}
P の中で break 文を実行するとこの繰り返しから抜けて次の実行文に移動する。break 文は通常 if 文に伴われて使用される。
条件 A が満たされている間命令文群 P を実行するという場合は と書く。P を実行している間に A が変わらなければ永久に P を実行し続ける。
while(A) {P}
変数 i の値の取りうる集合 T を与え、T の全てに対して変数 i の値を変えながら命令文群 P を実行したい場合は と書く。
for(i in T) {P}
命令文群の中で next 文を実行すると、残りをスキップし、再び命令文群の先頭から実行を開始する。
fname という名前の関数を定義するには
fname = function(X) {P}
P として実行可能な全ての R の命令が使える。関数定義に固有のものとして return 関数がある。return(z) を実行すると関数の実行が終了し、z を関数値として返す。return 文がない場合は、 P を全部実行し、最後に計算された値が関数値として返される。
関数の内部で代入文の左辺に現れる変数はローカル変数と呼ばれ、その代入結果は関数内部でのみ有効である。それ以外の変数はグローバル変数と呼ばれ、関数呼び出し時に記憶されている値が使われる。ローカル変数でも、代入文より前に使用されている場合はグローバル変数として扱われ、関数呼び出し時の値が使われる。グローバル変数に値を代入したい場合は永続代入演算子 <<- を使う。
次の例で、a は代入文の前はグローバル変数、その後はローカル変数として扱われる。関数実行が終了した時、関数内部での代入文はなかったものとして扱われる。
関数定義の中の a = 2 の代わりに永続代入演算子を使って a <<- 2 とした場合、最後の print(a) で表示されるのは 2 である。
> test = function() {
+ print(1:a)
+ a = 2
+ print(1:a)
+ }
> a = 3
> test()
[1] 1 2 3
[1] 1 2
> print(a)
[1] 3
R は入力された命令の文字列をその都度コンピュータ用のプログラムに置き換えて実行するインタープリタ型の処理系なので、 while 構文のような繰り返し命令をそのまま実行させるのはあまり効率的でない。apply 関数系と呼ばれる関数群など、同じような計算を繰り返し実行するために定義された関数を利用した方が計算時間が早い。例えば、データの平均値を計算する場合、for 文を使ってデータの累積値を計算してデータ数で割る、という方法の代わりに mean() 関数を利用するようなものである。
apply(df,rc,f) は、2次元のデータフレーム df の各行、あるいは各列に関数 f を適用する関数である。rc=1 ならば行ベクトル、rc=2 ならば列ベクトルが適用の対象になる。rc=c(1,2) とすると、全要素に適用される。
例
# 集計作業の例
> (ss <- matrix(runif(15),3,5)) [,1] [,2] [,3] [,4] [,5] [1,] 0.3165592 0.7425624 0.3938019 0.5726331 0.3156723 [2,] 0.7175305 0.8096806 0.9768243 0.9216174 0.5916798 [3,] 0.7836617 0.3334118 0.9215368 0.2805691 0.9762343
# 行和 > apply(ss,1,sum) [1] 2.341229 4.017333 3.295414
# 列平均 > apply(ss,2,mean) [1] 0.6059172 0.6285516 0.7640543 0.5916065 0.6278621
# 行毎の最大、最小値、四分位、... > apply(ss, 1, summary) [,1] [,2] [,3] Min. 0.3157 0.5917 0.2806 1st Qu. 0.3166 0.7175 0.3334 Median 0.3938 0.8097 0.7837 Mean 0.4682 0.8035 0.6591 3rd Qu. 0.5726 0.9216 0.9215 Max. 0.7426 0.9768 0.9762
# 各列の度数分布 > ss <- matrix(rpois(300,1),100,2) > apply(ss,2,table) [[1]] 0 1 2 3 4 45 36 15 3 1 [[2]] 0 1 2 3 4 37 34 20 7 2
# 全要素に同じ計算
> z = matrix(1:6,2) > apply(z, c(1,2), sqrt) [,1] [,2] [,3] [1,] 1.000000 1.732051 2.236068 [2,] 1.414214 2.000000 2.449490
sapply, lapply
sapply(v, f) は、ベクトル v の値を使って関数 f を評価した結果を返すベクトル値関数である。replicate() 関数も参照のこと。v が list の場合は lapply 関数を使う。
例
> piMonteCarlo = function(n=100) {
+ z = 4 * sqrt(1-runif(n)^2)
+ w = c(mean(z), sd(z)/sqrt(n))
+ names(w) = c("平均 ", "標準偏差")
+ return(w)
+ }
>
> sapply(c(100,1000,10000,100000), piMonteCarlo)
[,1] [,2] [,3] [,4]
平均 3.20056590 3.17135502 3.146771176 3.142442417
標準偏差 0.08570012 0.02734449 0.008861248 0.002815803
tapply各データごとに属性値が与えられている時、同じ属性値ごとに同じ計算を実行する。例えば、男女別に集計するなど。
replicatereplicate 関数は、シミュレーションのように、同じ計算を複数回繰り返してその結果をベクトルで返す。
> sapply(rep(1000,4), function(n) mean(rnorm(n))) # 正規乱数1000個の平均値を4通り計算する [1] 0.02675481 -0.03500802 -0.01315011 0.01796544
# 次のようにしても同じ > replicate(4, mean(rnorm(1000))) # 正規乱数1000個の平均値を4通り計算する
例題 表の確率が p のコインを n 回投げて相対度数を計算する、という実験を m回繰り返して、平均値と誤差を見積もり、ヒストグラムを描く
> p = 0.3 > n = 100 > m = 1000 > xa = replicate(m, mean(sample(1:0,prob=c(p,1-p),n,replace=T))) > c(mean(xa),sd(xa)/sqrt(m)) > hist(xa, freq=F)
print 関数は、配列データを編集して表示する。digits= オプションは表示有効桁数を指定する（臨時の options(digits=..) と同じ）。結果は配列として表示される。
> z = rnorm(8) > z [1] 1.5493767 -0.2196001 -1.0388382 -0.2145552 -0.7014241 0.8423489 -0.5512281 -0.3178271
# オプションを使わなければ、ただの配列表示
> print(z) [1] 1.5493767 -0.2196001 -1.0388382 -0.2145552 -0.7014241 0.8423489 -0.5512281 -0.3178271
# 有効桁を設定
> print(z, digits=2) [1] 1.55 -0.22 -1.04 -0.21 -0.70 0.84 -0.55 -0.32
cat 関数は、引数を文字列に変換し半角スペースでつなげた一つの文字列として表示する。自動的に改行はしないので、改行させたい場合は、最後に改行の制御コード "\n" を文字列として指定する。sep= オプションを使うと、半角スペースを挟まない文字列を生成して表示する。
> a <- 1; b <- 2 > cat(a, "+", b, "=", a+b, "\n") 1 + 2 = 3
# print との違い
> cat(1:5) 1 2 3 4 5
> cat(1:5, sep="") 12345 > print(1:5) [1] 1 2 3 4 5
sprintf() は、与えられた書式にしたがって、引数を表示する関数である。仕様は C 言語の printf() 関数と同じ。結果は一つの文字列データとして保存されるので、「""」で囲まれて表示される。「""」を取り除きたければ、結果を cat() 関数で変換する。
> a = 1; b = 2:3;
> sprintf("%d + %d = %d\n", a, b, a+b)
[1] "1 + 2 = 3\n" "1 + 3 = 4\n"
> cat(sprintf("%d + %d = %d\n", a, b, a+b))
1 + 2 = 3
1 + 3 = 4
sink() 関数は、Rconsoleの表示結果と同じものをファイルにコピーする。
> sink("result.txt") # ここに書かれた命令の計算結果がすべてファイル保存される
> sink() # 何も指定しなければ「マイドキュメント」に保存される
write 関数は、計算結果を R Console に表示することなく、直接ファイルに保存する。引数としてファイル名を指定すると MyDocument にそのファイルが作成され、データが保存される。
MyDocument 以外の場所に保存したい場合はディレクトリの変更が必要である。R Console をアクティブにして「ファイル／ディレクトリの変更」メニューを選択、「フォルダの参照」ウィンドウで指定する。
write(t(A), "data.txt", ncol=ncol(A)) # 行列Aの転置（行優先）を書き出す
write.table 関数は、表形式の統計データ（データフレーム）を保存する。
> name = c("相田", "池田", "太田", "蒲田", "島田", "高田", "富田", "古田", "松田")
> height = c(177.1, 171.7, 169.2, 173.4, 181.9, 171.9, 165.0, 175.0, 175.4)
> weight = c(66.8, 72.6, 77.8, 75.7, 81.2, 70.3, 75.4, 81.5, 72.2)
> (df = as.data.frame(list(name=name, height=height, weight=weight)))
name height weight
1 相田 177.1 66.8
2 池田 171.7 72.6
3 太田 169.2 77.8
4 蒲田 173.4 75.7
5 島田 181.9 81.2
6 高田 171.9 70.3
7 富田 165.0 75.4
8 古田 175.0 81.5
9 松田 175.4 72.2
> write.table(df, "datafr.txt")
scan() は、ベクトルデータを読み込む関数である。オプションを何も指定せず空の引数とした場合は、数値入力が期待される。複数のデータをスペース区切りで１行にして入力することができる。改行すると、次のデータ入力を促す数値（次の数が入力される配列の番号）が表示されるが、そこで直ちに Enter キーを押すと入力を終了する。
例 キーボードからの入力
> x <- scan() 1: 1.2 3 4 99 5: Read 4 items > x [1] 1.2 3.0 4.0 99.0
> (x = scan()) 1: # 何も入力せず Read 0 items numeric(0)
what="" オプションを指定すると（半角空白文字を含まない）文字列が入力される。
> z = scan(what="") 1: a b あ d 5: Read 4 items > z [1] "a" "b" "あ" "d"
関数の実行中に、キーボードからデータを入力させたい場合は、readline() 関数を使う。入力された行を文字列として受け取るので、strsplit(z, "[ ,]") 関数で文字列に分解する。"[ ,]" は正規表現で、データはスペース区切り「または」カンマ区切りとする、という記号である。文字列は list 構造になっているので、unlist 関数でベクトル化し、as.numeric 関数を使って数値に置き換えるという作業が必要になる。
# 入力データの平均と標準偏差を計算する
> test <- function() {
+ z <- readline()
+ w <- as.numeric(unlist(strsplit(z, "[ ,]")))
+ return(c(mean(w), sd(w)))
+ }
> test()
1 2 3
[1] 2 1
1, 2, 3
[1] 2 1
> test()
1, 2, 3
[1] NA NA
最後の入力例では、４つのデータ区切り記号が入力されているので「"1" "" "2" "" "3"」という５つの文字列が入力されたとみなされるので、意図した計算が行われない。
scan 関数は、テキストファイルの保存したスペース区切りの数値データを読み込む。sep="," オプションを指定することでカンマ区切りのデータにも対応している。 what="" オプションを指定することで文字型データを読み込むこともできる。
> height
[1] 177.1 171.7 169.2 173.4 181.9 171.9 165.0 175.0 175.4
> write(height, "textht.txt")
> (data = scan("textht.txt"))
Read 9 items
[1] 177.1 171.7 169.2 173.4 181.9 171.9 165.0 175.0 175.4
# カンマ区切りデータの場合
> weight
[1] 66.8 72.6 77.8 75.7 81.2 70.3 75.4 81.5 72.2
> write(height, "textht.txt", sep=",")
> (commadata = scan("textht.txt", sep=","))
Read 9 items
[1] 177.1 171.7 169.2 173.4 181.9 171.9 165.0 175.0 175.4
データフレームは、統計分析の為のデータ構造で、各列ごとにデータ型を決めることができる行列（表）である。R で直接作成することもできるが、大量のデータを扱う場合は、Excelシートで整形してからRに取り込むとよい。read.table 関数は、横に項目を並べ、縦に標本を並べた統計データ表を入力する。
Excelシート上でデータ表が完成したら、項目名と一緒に選択してコピーし、クリップボードに一時保存したものを read.table() 関数で読み込む。
保存された Excelファイルを読み込む為には、その為のライブラリを必要とするので、その代わり、Excel ファイルを csv 形式で保存し、それを read.csv() 関数で読み込む。
カレントディレクトリを変更するには、R Consoleウィンドウの「ファイル」メニューから「ディレクトリの変更」サブメニューを選択、「フォルダの参照」ウィンドウで指定する。
> (df = as.data.frame(list(name=name, height=height, weight=weight)))
name height weight
1 相田 177.1 66.8
2 池田 171.7 72.6
3 太田 169.2 77.8
4 蒲田 173.4 75.7
5 島田 181.9 81.2
6 高田 171.9 70.3
7 富田 165.0 75.4
8 古田 175.0 81.5
9 松田 175.4 72.2
> write.table(df, "datafr.txt")
> (newdf = read.table("datafr.txt"))
name height weight
1 相田 177.1 66.8
2 池田 171.7 72.6
3 太田 169.2 77.8
4 蒲田 173.4 75.7
5 島田 181.9 81.2
6 高田 171.9 70.3
7 富田 165.0 75.4
8 古田 175.0 81.5
9 松田 175.4 72.2
# クリップボード経由
> df <- read.table("clipboard", header=FALSE) # Excelデータをコピーした後実行する
# csv ファイル入力
> dfs <- read.csv("testdata.csv", header=TRUE, skip=2)
# タブ区切りテキストファイル入力
> scan("data.txt", list(name="", height=1, weight=1), sep="\t")
注釈
header=TRUE/FALSE は、読み込むデータファイルのデータの前に項目名を入力した行があるかないかを指定するオプションである。skip= は、データテーブルの前に注釈行などが入っている場合などに、データ行でない行を読み飛ばすためのオプションである。sep="\t" はデータがタブ区切りであることを指定するオプションである。1変数関数の方程式 f(x) = 0 の解を見つける関数として uniroot(), polyroot() がある。
方程式の解：uniroot() は、連続関数の関数記号と、両端で関数値の符号が異なる区間を与えて、その区間に含まれる解（中間値の定理で保証される）を一つを計算する。結果は、誤差や反復回数などの情報を合わせてリストで与えられる。
> uniroot(sin, c(1,11)) $root [1] 3.141593 $f.root [1] 9.76346e-08 $iter [1] 7 $estim.prec [1] 6.103516e-05
> c(sin(1), sin(11)) [1] 0.8414710 -0.9999902
区間の端点で同符号となる場合は、エラーメッセージが表示される。
関数は、すでに定義されている関数の関数記号を書くか、短いものならば、その場で関数定義を書いてもよい。
例 関数記号の代わりに関数定義を使った例：標準正規分布の上側 0.1% 点の計算
> uniroot(function(x) pnorm(x)-0.999, c(1,5)) $root [1] 3.090241 $f.root [1] 3.026668e-08 $iter [1] 8 $estim.prec [1] 6.103516e-05 > qnorm(0.999) ## 計算の精度チェック [1] 3.090232
$iter は収束するまでの反復回数である。print() 関数を使って反復計算の過程を調べることも出来る。
> uniroot(function(x) pnorm(print(x))-0.999, c(1,5)) [1] 1 [1] 5 [1] 4.974795 [1] 2.987398 [1] 3.562238 [1] 3.178618 [1] 3.10404 [1] 3.089814 [1] 3.090241 [1] 3.09018 [1] 3.090241
多項式方程式の解：関数が多項式の場合は、polyroot() を使って、すべての解を同時に求めることが出来る。引数は、係数ベクトルとし、結果は複素数値ベクトルである。実数解でも、計算誤差に由来する「+0i, -0i」という複素数表示になる。
例 x^3 - 2x = 0 （実数解） > a = c(0, -2, 0, 1) > (r = polyroot(a)) [1] 0.000000+0i 1.414214-0i -1.414214+0i 検算（誤差） > sapply(1:3, function(i) sum(a * r[i]^(0:n))) [1] 0.000000e+00+0e+00i -4.440892e-16-0e+00i 4.440892e-16+0e+00i 例 x^3 - 1 = 0 （１の3乗根） > a = c(-1, 0, 0, 1) > (r = polyroot(a)) [1] 1.0-0.0000000i -0.5+0.8660254i -0.5-0.8660254i
検算（誤差）
> sapply(1:3, function(i) sum(a * r[i]^(0:n))) [1] 0.000000e+00-2.298900e-20i -2.220446e-16+1.110223e-16i [3] -2.220446e-16-1.110223e-16i
検算（ゼロ判定） > abs(sapply(1:3, function(i) sum(a * r[i]^(0:n)))) < 10^-10 [1] TRUE TRUE TRUE
連立一次方程式は、solve 関数を使って解くことが出来る。引数は、係数行列と定数ベクトルで、解ベクトルを返す。
%*% は（ベクトルも含めて）行列同士の積を計算する演算子である。> (A = rbind(c(-1,4,2), c(3,-2,1), c(1,1,1)))
[,1] [,2] [,3]
[1,] -1 4 2
[2,] 3 -2 1
[3,] 1 1 1
> b = c(1,2,-1)
> (r = solve(A, b))
[1] -3 -3 5
検算
> A %*% r
[,1]
[1,] 1
[2,] 2
[3,] -1
変数の数を n としたとき、定数ベクトルの代わりに n x m 行列を与えると、その行列の m 本の列ベクトルを定数ベクトルとする m 個の連立一次方程式を同時に解くことができる。
> (r = solve(A, cbind(c(1,2,-1), c(1,0,0))))
[,1] [,2]
[1,] -3 -0.6
[2,] -3 -0.4
[3,] 5 1.0
検算
> A %*% r
[,1] [,2]
[1,] 1 1.000000e+00
[2,] 2 0.000000e+00
[3,] -1 -2.220446e-16
定数ベクトル（行列）を単位行列にすると、係数行列の逆行列が計算出来る。この場合、単位行列を省略して、solve(A) と書けば、A の逆行列が計算される。
> solve(A) # solve(A, diag(3)) と同じ
[,1] [,2] [,3]
[1,] -0.6 -0.4 1.6
[2,] -0.4 -0.6 1.4
[3,] 1.0 1.0 -2.0
# 検算
> solve(A) %*% A
[,1] [,2] [,3]
[1,] 1.000000e+00 -2.220446e-16 -4.440892e-16
[2,] 2.220446e-16 1.000000e+00 -2.220446e-16
[3,] 0.000000e+00 0.000000e+00 1.000000e+00
係数行列の i 列目を定数ベクトルに置き換えた行列の行列式と係数行列の行列式の比が解になる、というクラメルの公式を利用した連立方程式を解くことも出来る。係数行列の行列式がゼロで ないことが必要である。
> A = rbind(c(-1,4,2), c(3,-2,1), c(1,1,1)) > b = c(0,0,1) > det(cbind(b,A[,2:3])) / det(A) [1] 1.6 > det(cbind(A[,1],b,A[,3])) / det(A) [1] 1.4 > det(cbind(A[,1:2],b)) / det(A) [1] -2
solve 関数を使わずに、連立一次方程式をガウスの消去法（掃き出し法ともいう）を使って解く手順を示す。これは、方程式の１つの式を定数倍しても、ある式の定数倍を別の式に加えても、式の順番を入れ換えても解は変わらない、という性質をだけを使って式を変形して解を導く、という方法である。
この方法は、行列の（行に関する）初等変形（ある行を定数倍する、ある行の定数倍を別の行に足す、ある行と別の行を入れ換える）だけを使って、係数行列を単位行列に変換することが出来る、ということを利用している。
> (A = rbind(c(1,5,3), c(2,3,-1), c(-2,2,5)))
[,1] [,2] [,3]
[1,] 1 5 3
[2,] 2 3 -1
[3,] -2 2 5
> x = c(3, -3, 4) # 問題作成（解から定数ベクトルを計算）
> b = A %*% x
> (B = cbind(A, b)) # 拡大係数行列
[,1] [,2] [,3] [,4]
[1,] 1 5 3 0
[2,] 2 3 -1 -7
[3,] -2 2 5 8
> k = 1; D = diag(3); D[,k] = -B[,k] / B[k,k]; D[k,k] = 1/B[k,k]
> (B = D %*% B) # １列目掃き出し計算
[,1] [,2] [,3] [,4]
[1,] 1 5 3 0
[2,] 0 -7 -7 -7
[3,] 0 12 11 8
> k = 2; D = diag(3); D[,k] = -B[,k] / B[k,k]; D[k,k] = 1/B[k,k]
> (B = D %*% B) # ２列目掃き出し計算
[,1] [,2] [,3] [,4]
[1,] 1 0 -2 -5
[2,] 0 1 1 1
[3,] 0 0 -1 -4
> k = 3; D = diag(3); D[,k] = -B[,k] / B[k,k]; D[k,k] = 1/B[k,k]
> (B = D %*% B) # ３列目掃き出し計算
[,1] [,2] [,3] [,4]
[1,] 1 0 0 3
[2,] 0 1 0 -3
[3,] 0 0 1 4
> # 検算
> A %*% B[,4] - b
[,1]
[1,] 0
[2,] 0
[3,] 0
掃き出し計算の際、対角要素が０になる場合は、その列で、０の下にある要素の中の最大値を見つけ、その行と入れ替えるという操作が必要になる。n 変数の場合、k 回目の掃き出しで対角要素が０でないかどうか、０ならば０でない行と入れ替える、という作業は次のようにすれば良い。
order 関数は配列を昇順あるいは降順に並べ替えた時の添え字を返す。decreasing=T オプションは降順を指定する。if(B[k,k] == 0) (B[k:n,] = B[order(abs(B[k:n,k]), decreasing=T)+(k-1),])
一般の n 元連立方程式を解くために関数としてまとめると、次のように書ける。
#### A x = b を解く。一次従属の場合は Inf を返す。
mySolve = function(A, b, n) { B = cbind(A, b) for(k in 1:n) {
if(B[k,k] == 0) (B[k:n,] = B[order(abs(B[k:n,k]), decreasing=T)+(k-1),])
if(B[k,k] == 0) return(rep(Inf, n)) D = diag(n); D[,k] = -B[,k] / B[k,k]; D[k,k] = 1/B[k,k] B = D %*% B } return(B[,n+1]) }
「x2の微分は 2x」というように、数式のまま微分するにはD(expression(), 変数名) 関数を使う。expression() は、簡易数式をタイプセットする関数である。微分する変数は "" で囲んだ文字列によって指定する。
例
> D(expression(x^2 * e^(-x)),"x") 2 * x * e^(-x) - x^2 * (e^(-x) * log(e)) > D(expression(1/sqrt(1 + r^2)),"r") # 分数は「/」を使う -(0.5 * (2 * r * (1 + r^2)^-0.5)/sqrt(1 + r^2)^2)
expression() は、簡易数式をタイプセットする関数である。たとえば、べき乗の代わりに「^」、平方根の代わりに「sqrt(..)」と書け ば、それを正しい数式と認識する。詳しくはここを参照。 integrate()関数を使って定積分を計算することができる。最初の引数として、被積分関数を定義し、2，3番目の引 数で積分範囲を指定する。無限はInfで表す。システムで定義されている数学関数、統計関数、あるいはそれらの組合せを指定しても 良い。結果は誤差評価とともにlistで返される。数式の中で使う 場合は「$value」を付けて関数値だけ取り出す必要があ る。 精度はあ まり 良くない。
> integrate(function(x) {x^2}, 0, 1)
0.3333333 with absolute error < 3.7e-15
> integrate(function(x) {exp(-x^2/2)/sqrt(2*pi)}, -Inf, Inf)
1 with absolute error < 9.4e-05
> upper25 = integrate(dnorm, 1.96, Inf)
> upper25
0.0249979 with absolute error < 1.9e-05
> upper25$value
[1] 0.0249979
> upper25$abs.error
[1] 1.894545e-05
あるいは
> integrate(dnorm, 1.96, Inf)$value
[1] 0.0249979
関数定義が何行にも渡るような場合に備えて、関数を定義してから、関数名を引数とすることもできる。
> g <- function(t, mu=2) mu * exp(-mu*t) > integrate(g, 0, Inf) ## 広義積分も計算できる 1 with absolute error < 5e-07 > integrate(g, mu=5, 0, Inf) ## パラメータの指定の仕方 1 with absolute error < 2.3e-06
注意：関数定義の中に条件分岐がある場合はうまくいかない。
> g <- function(t,mu=2) {
+ if(t > 0) return(mu * exp(-mu*t) / 2)
+ if(t <= 0) return(mu * exp(mu*t) / 2)
+ }
>
> g(2) // 計算だけならば受け付けられる
[1] 0.01831564
> g(-1)
[1] 0.1353353
> integrate(g,-2,2)
27.28992 with absolute error < 3e-13 ## 最初の関数だけを使って計算される
警告メッセージ：
1: In if (t > 0) mu * exp(-mu * t)/2 :
条件が長さが2以上なので，最初の一つだけが使われます
2: In if (t <= 0) mu * exp(mu * t)/2 :
条件が長さが2以上なので，最初の一つだけが使われます
関数定義が条件分岐によって定義されている場合、定義区間ごとに別の関数を定義して、区間毎に積分した値を足し合わせる必要があるが、そうする代わりに、ifelse()関 数を使うことにより、同じことを実現することが出来る。
# 別々に定義 > gp <- function(t,mu=2) mu * exp(-mu*t) / 2 ## [0,Inf) で定義 > gm <- function(t,mu=2) mu * exp(mu*t) / 2 ## (-Inf,0] で定義 > integrate(gp,0,Inf) 0.5 with absolute error < 7.7e-05 > integrate(gm,-Inf,0) 0.5 with absolute error < 7.7e-05 # この場合は、ifelseを使って次のように定義できる。 g <- function(t,mu=2) exp(ifelse(t>0,-1,1) * mu*t) * mu / 2 > integrate(g,-2,2) 0.9816844 with absolute error < 1.1e-14
数論の問題は、配列計算では対処できないことが多く、Cのようなプログラミング言語と同じ考え方でプログラムしなければいけない。したがって、逐次翻訳型のRでは計算は遅いので、本計算の前のアルゴリズムチェック用と考えたほうが良い。
自然数 a が自然数 n の原始根であるとは、am = 1 mod n となる最小の m が n-1 である、あるいは ak mod n (k=1,2,...,n-1) がすべて異なることをいう。例えば、3は19の原始根であることは次の計算でわかる。
> n = 19 > 3^(1:(n-1)) %% n [1] 3 9 8 5 15 7 2 6 18 16 10 11 14 4 12 17 13 1
a が n の原始根であるかどうか知るには、m = 1 から順番に am = 1 mod n となるまで順番に計算し、それが n-1 になるかどうかを調べるしかない。次は a, n に対して、am = 1 mod n となる最小の m を計算する関数である。
primRoot = function(a=3, n=19) {
b = a*a
m = 2
while(b %% n != 1) {
b = (b*a) %% n
m = m+1
}
return(m)
}
n のすべての原始根を求めるには sapply 関数を使えば良い。
> sapply(2:18, primRoot, n=19) [1] 18 18 9 9 9 3 6 9 18 3 6 18 18 18 9 9 2 > n = 19 > which(sapply(2:(n-1), primRoot, n=n) == n-1) + 1 [1] 2 3 10 13 14 15
ベーブルースとハンクアーロンのホームラン数がそれぞれ 714, 715本である。714, 715 を素因数分解すると 714 = 2 x 3 x 7 x 17, 715 = 5 x 11 x 13 であり、それらの素因数の総和はいずれも 29 になるという特異な性質があることが分かった。このことから、隣り合った整数で、その素因数の総和が等しいペアのことをルースアーロンペアと名付けられた。
10000 以下のペアは 20、次の10000個の中には６個しかない。しかし、無限にあるかどうかはわかっていない（が、その逆数の和は収束することがわかっている）。
ルースアーロンペアを知るには、個別に素因数を求めてその和が等しいかどうかをチェックするしかない。
# 素因数を（重複も含めて）リストアップ
> primefactordup = function(n) { factor = numeric(0) for(k in 2:sqrt(n)) { if(n %% k != 0) next repeat { factor = c(factor, k) n = n / k; if(n %% k != 0) break } } if(n != 1) factor = c(factor, n) return(factor) }
# ペアの計算
RuthAaronPair = function(n) { pf = primefactor(2) m = 0 for(i in 3:n) { pfg = primefactordup(i) if(sum(pf) == sum(pfg)) { m = m + 1 print(c(m, i-1, sum(pfg))) } pf = pfg } }
# 実行例
> RuthAaronPair(100) [1] 1 5 5 [1] 2 8 6 [1] 3 15 8 ## pfs(15) = 3 + 5, pfs(16) = 2 + 2 + 2 + 2 [1] 4 77 18
最小値を計算する関数として、optimize() がある。指定された解の探索区間に対して、黄金分割に基づくフィボナッチ探索で最適解を探索する。
> optimize(function(x) x^2*(x-1), c(0,2)) $minimum [1] 0.6666728 $objective [1] -0.1481481
フィボナッチ探索の収束の様子を見るために、print() 関数を利用することが出来る。
> optimize(function(x) x^2*(print(x)-1), c(0,2)) [1] 0.763932 [1] 1.236068 [1] 0.472136 [1] 0.6414298 [1] 0.6572841 [1] 0.6651937 [1] 0.6668159 [1] 0.6666728 [1] 0.6666321 [1] 0.6667135 [1] 0.6666728 $minimum [1] 0.6666728 $objective [1] -0.1481481
最初の2点は、黄金比を使って次の式で決まる。
> (phi = (sqrt(5) - 1) / 2) [1] 0.618034
> a = 0; b = 2 > c(a+(1-phi)*(b-a), a+phi*(b-a)) [1] 0.763932 1.236068
最大値を求めたい場合は、-f(x) の最小値を計算すれば良い。ただし、関数記号にマイナスを付けるだけではエラーになる。
> optimize(-dnorm, c(-1,1)) 以下にエラー -dnorm : 単項演算子に対する引数が不正です > optimize(function(x) -dnorm(x), c(-1,1)) $minimum [1] -8.326673e-17 $objective [1] -0.3989423
optim 関数も使える。
> f = function(x) x+1/x
> optim(2, f)
$par
[1] 1
$value
[1] 2
$counts
function gradient
30 NA
$convergence
[1] 0
$message
NULL
警告メッセージ：
In optim(2, f) : Nelder-Mead 法による一次元最適化は信頼できません:
"Brent" 法を用いるか optimize() を直接使ってください
最大値を計算する場合は、controlパラメータを設定することで対処するか、関数定義の符号を逆転させる。ただし、引数の関数名の前にマイナス符号を付 ける書き方は認められない。
> optim(3, dnorm, mean=5, sd=2, control=list(fnscale=-1)) $par [1] 5.000098 $value [1] 0.1994711 ...
線形計画法は boot ライブラリの simplex 関数を使う。「≦制約式」、「≧制約式」「＝制約式」の順番に、（係数行列、定数ベクトル）のペアで指定する。
> library(boot)
> A = rbind(c(1,2),c(3,1)) // 制約式の係数行列
> b = c(10,8) // 制約式の定数ベクトル
> c = c(1,1) // 目的関数の係数ベクトル
> simplex(c, A, b, maxi=T)
Linear Programming Results
Call : simplex(a = c, A1 = A, b1 = b, maxi = T)
Maximization Problem with Objective Function Coefficients
x1 x2
1 1
Optimal solution has the following values
x1 x2
1.2 4.4
The optimal value of the objective function is 5.6.
> sm = simplex(c, A, b, maxi=T) // 目的関数最大化
> cat("F(",sm$soln[[1]],",",sm$soln[[2]], ") =",sm$value,"\n")
F( 1.2 , 4.4 ) = 5.6
整数計画法も解けるパッケージ lpSolve を使うにはインストールする必要がある。インストールがうまくいかない場合は、次を試すと良い（R のアイコンを右クリックして、プロパティを表示させ、「リンク先」のパラメータに「--internet2」を追加する）。
> install.packages("lpSolve");
パッケージを ‘C:/Users/sakasegawa/Documents/R/win-library/2.15’ 中にインストールします (‘lib’ が指定されていないので)
--- このセッションで使うために、CRANのミラーサイトを選んでください ---
URL 'http://cran.md.tsukuba.ac.jp/bin/windows/contrib/2.15/lpSolve_5.6.6.zip' を試しています
Content type 'application/zip' length 678244 bytes (662 Kb)
開かれた URL
downloaded 662 Kb
パッケージ ‘lpSolve’ は無事に展開され、MD5 サムもチェックされました
ダウンロードされたパッケージは、以下にあります
C:\Users\sakasegawa\AppData\Local\Temp\Rtmp2lOMC5\downloaded_packages
> library(lpSolve)
使い方は以下の通り
> f.in = c(1,1) # max (c,x), s.t. Ax ≦ b を解く
> f.mat = rbind(c(1,2),c(3,1)) # = A
> f.dir = c("<=", "<=") # 「≦」
> f.rhs = c(10,8) # = b
> soln = lp ("max", f.in, f.mat, f.dir, f.rhs) ## "min" がデフォルト
> soln$solution
[1] 1.2 4.4
整数計画法ならば「all.int=TRUE」オプションを指定する。一部の変数だけの整数制約ならば、「int.vec=c(1,2,3)」 のように、整数条件の変数の添え字を指定する。
同じように、0-1計画問題ならば「all.bin=TRUE」オプションを指定する一部の変数だけが2値変数なら ば「binary.vec」によって、添え字集合を指定する。
0
> f.in = c(10,6,3,8,6)
> f.mat = matrix(c(5,3,2,3,4), nrow=1) // 制約式一つでも必ずmatrix指定、ベクトル指定は不可
> f.dir = "<="
> f.rhs = 8
> lp ("max", f.in, f.mat, f.dir, f.rhs, all.int=T)$solution // 整数計画法
[1] 0 0 1 2 0
> lp ("max", f.in, f.mat, f.dir, f.rhs, all.bin=T)$solution // 0-1 計画法
[1] 1 0 0 1 0
主要な確率分布については、分布固有の名称の前に「d」「p」「q」「r」を付与することで、密度関数（離散の場合は確率質量関数）、累積分布関数、その逆関数、乱数を生成する関数、を得ることが出来る。
正規分布の固有名称は norm なので、次のような4通りの関数が利用できる。
| 関数名 | 説明 |
|---|---|
pnorm(x, m, s) |
累積分布関数の値（m=0,s=1は省略可） |
dnorm(x, m, s) |
密度関数の値 |
qnorm(x, m, s) |
逆関数の値（分位点） |
rnorm(n, m, s) |
正規乱数（nは生成する乱数の個数） |
qxxx()関数を使うと、分布のパーセント点が計算できる。例えば、正規分布の場合、qnorm 関数を使って正規分布表が作れる。
p <- c(0.005, 0.025, 0.05, 0.975, 0.995) qnorm(p) # 標準正規分布の分位点 [1] -2.575829 -1.959964 -1.644854 1.959964 2.575829
Ｒで用意されている確率分布の主要なものをリストする。詳しくは「?xxx」と入力して、オンラインマニュアルを参照のこと。
|
|
密度関数のグラフは curve() 関数を使って描くことができる。確率質量関数は type="h" オプションを指定した plot 関数で描くことができるが、棒を太くみせたければ barplot() 関数を使うこともできる。
次は、正規分布と t分布の密度関数を重ねて描いた例である。
curve(dnorm, -4, 4, main="正規分布と t分布")
curve(dt(x,10), add=T, col=2)
curve(dt(x,5), add=T, col=3)
curve(dt(x,2), add=T, col=4)
legend(1.8, 0.4, c("normal","t(10)","t(5)","t(2)"), lwd=rep(1,4), col=1:4)
次は、２項分布とポアソン分布の確率質量関数の例である。２項分布のパラメータ n, p の n が大きく、np がほどほどの時、ポアソン分布で近似できることを視覚的に確認するものである。グラフを重ねるには低水準作図関数の lines 関数を使う。
dbinom 関数は2項分布の質量関数、qbinom 関数は2項分布の累積分布関数の逆関数で、この場合は下側 0.1%点、上側 99.9%点を求めている。plot 関数の type="h" オプションは棒グラフを描き、lwd= オプションは線の太さを指定する。qpois 関数はポアソン分布の質量関数。n = 1000
p = 0.02
na = qbinom(0.001, n, p)
nb = qbinom(0.999, n, p)
plot(na:nb, dbinom(na:nb,n,p), type="h", lwd=3, main=paste("n =", n, ", p =",p))
lines(na:nb, dpois(na:nb, n*p), col = 2)
abline(h=0)
text(30, 0.085, "折れ線はポアソン分布")
確率分布に従う乱数は、分布名の前に「r」を付ければ良い。
次は、標準正規分布に従う乱数を生成する例である。
> rnorm(10) [1] -0.14221818 -1.40954395 -1.09566066 0.82089141 0.78028615 -0.65467685 [7] -0.55899167 0.06699198 -1.83961413 -0.83795562
Ｒを立ち上げて最初に rnorm(10) を実行すると、結果は常に上のような数が表示される。並び方は一見ランダムのようであるが、何回でも同じ数列が再現される。
rnorm() は、あるプログラムによって生成される規則的数列を返す関数で、その規則性が分かってしまえば、次の数が計算でき、乱数でもなんでもないが、生成された数列だけをみてその規則性が分からなければそれを乱数と言っても良いだろう。そのような性質を持った数列を擬似乱数という。
擬似乱数を生成するプログラムを立ち上げたとき、初期状態が同じならば常に同じ数列を返す。その初期状態を変えることによって生成される数列を変えることが出来る。初期状態を変えるために、set.seed() 関数を使う。引数は適当に大きな数を指定する。
> set.seed(11212221) > rnorm(10) [1] -0.7741669 0.1801611 0.3575552 1.8654200 -1.8474021 -0.1752026 [7] -1.0366205 -0.7561760 -0.2727841 1.2816695 > hist(rnorm(1000), freq=F)
a1, ..., am をそれぞれ確率 p1, ..., pm で取るような、一般の離散分布に従う乱数は sample() 関数によって生成できる。
a1, ..., am をベクトルにしたものを最初の引数とし、prob= オプションで、p1, ..., pm をベクトルにしたものを与える。ただし、一様分布の場合は省略することが出来る。
replace=TRUE オプションを指定することにより、独立同分布のサンプリングが実行される。デフォルトでは replace=FALSE になっており、これは非復元サンプリングを実行する。
> sample(1:20, 20, replace=TRUE)
[1] 7 14 14 11 20 4 17 11 19 5 2 5 14 12 12 12 12 15 4 3
# 非復元抽出（この場合はランダム置換）
> sample(1:20, 20) # sample(20) としても同じ
[1] 1 17 2 4 13 3 19 5 9 15 18 14 20 12 16 10 11 7 6 8
# コイン投げ実験
> sample(c("T","H"), 10, replace=TRUE)
[1] "H" "T" "T" "T" "T" "H" "H" "H" "T" "H"
# ゆがんだサイコロ振り
> sample(1:6, 20, prob=c(2,1,1,1,1,2), replace=TRUE)
[1] 1 1 4 1 6 2 4 5 3 5 4 2 1 5 6 4 3 4 3 1
> barplot(table(sample(1:6, 1000, prob=c(2,1,1,1,1,2), replace=TRUE)))
prob= オプションは、確率ではなく、比率で与えてもよい。 
行列の行和、列和が与えられているとき、その行和列和に等しくなるような行列の各要素をランダムに生成する関数 r2dtable() がある。
最初の引数は、生成する行列の個数、2番目の引数で行和ベクトル、3番目の引数で列和ベクトルを指定する。
> r2dtable(2, c(10,12,18), c(20,10,10))
[[1]]
[,1] [,2] [,3]
[1,] 5 2 3
[2,] 8 2 2
[3,] 7 6 5
[[2]]
[,1] [,2] [,3]
[1,] 6 1 3
[2,] 5 5 2
[3,] 9 4 5
f(x,y) = f(x | y)f(y) と置き換え、条件付き分布が正規分布に従うことを利用する。
rnorm2 <- function(n, mx, sx, my, sy, r) {
x <- rnorm(n)
y <- rnorm(n)
x <- x * sqrt(1-r^2) + r * y
return(list(x = sx*x+mx, y = sy*y+my))
}
n <- 1000; mx <- 10; sx <- 3; my <- 20; sy <- 5; r <- -0.8
z <- rnorm2(n, mx, sx, my, sy, r)
plot(z$x, z$y, main=paste("2変量正規分布 ρ=", round(cor(z$x,z$y),3)))
return(list(...)) のようにすると良い。z <- rnorm2(... として計算した後、list= の左辺の変数を使って z$x, z$y のようにして二つのベクトルを参照することが出来る。cor() は、相関係数を計算する関数である。round() は、小数の桁数を指定して、四捨五入する関数である。paste() は文字列を連結する関数である。数値は文字列に置き換わる。 
n 個のうち m 個が１の 0, 1 がランダムに並んだベクトルを生成する。
sample() 関数を使って、m個の１と n-m 個の０の配列のランダム置換を生成すればよい。
> n = 10 > m = 4 > sample(c(rep(1,m), rep(0,n-m))) [1] 1 1 1 0 0 0 0 0 1 0 > sample(c(rep(1,m), rep(0,n-m))) [1] 1 0 1 0 1 0 0 0 0 1
パッケージ mvtnorm の中に、rmvnorm() 関数があり、平均ベクトルと共分散行列を与えて、多変量正規分布に従う乱数を生成することが出来る。
パッケージを利用するには、install.package() 関数と、library() 関数を使う。
## 多変量正規分布は、パッケージ mvtnorm に納められている
> install.packages("mvtnorm")
> library(mvtnorm)
## norm と書くところに mvnorm と書けばよい（d-, p-, r-)
## 2変量正規乱数の生成
> M <- c(1,4)
> S <- matrix(c(1,1.2,1.2,4),2)
> rmvnorm(10, M, S))
[,1] [,2]
[1,] -0.3069163 2.5247827
[2,] 2.3129524 5.5592589
[3,] 1.8486602 2.9057861
[4,] 0.4908641 7.9475700
（以下略）
## 3次元の散布図を描く：パッケージ scatterplot3d を使う
> install.packages("scatterplot3d")
> library("scatterplot3d")
> z <- rmvnorm(1000, M, S)
> scatterplot3d(z[,1], z[,2], dmvnorm(z,M,S))
## 結合密度関数を描く：persp 関数を利用
m = 60
M = c(1,4)
S = c(1,2)
R = matrix(c(1,0.6,0.6,1),2)
SGM = matrix(S,2) %*% matrix(S,1) * R
x = seq(M[1]-3*S[1], M[1]+3*S[1], length=m)
y = seq(M[2]-3*S[2], M[2]+3*S[2], length=m)
z = matrix(rep(0,m*m),m)
for(i in 1:m) for(j in 1:m) z[i,j] = dmvnorm(c(x[i],y[j]),M,SGM)
persp(x, y, z, theta=30, ticktype="detailed")
## 次のようにしたいのだが、エラーになる。
f = function(xx,yy) dmvnorm(c(xx,yy),M,SGM)
zz = outer(x, y, f)
ライブラリ MASS のなかにある mvrnorm() 関数を使ってもよい。平均値ベクトル、共分散行列を与えると、m次元乱数を表にして返してくる。
> library(MASS)
> x <- mvrnorm(1000, c(0,0), matrix(c(1,0.5,0.5,1),2))
> plot(x[,1], x[,2] main=paste("2変量正規分布 ρ=", round(cor(x[,1],x[,2]),3)))
文字列を扱う関数をまとめる。
| 関数名 | 説明 |
|---|---|
substr, substring |
文字列の部分列を取り出す(substring はベクトル化関数) |
strsplit |
文字列の要素を、与えられたパターンに従って取り出す |
paste |
文字列の結合 |
nchar |
文字列の長さをカウント、ベクトルでも良い |
charmatch |
部分文字列の照合 |
gsub, sub |
文字列の置き換え(sub は最初の文字列のみ） |
as.numeric |
文字列として与えられた数字を数に変換する |
noquote |
「"」を付けずに表示 |
幾つかの関数の使用例を示す。
> (text = date()) [1] "Sat Jan 28 22:40:45 2017"
> substr(text, 12, 19) [1] "22:40:45" > substring(date(), c(1,5,9,21), c(3,7,10,24)) [1] "Sat" "Jan" "28" "2017"
> strsplit(date(), " ") [[1]] [1] "Sat" "Jan" "28" "22:43:38" "2017" > paste("Today is", date()) [1] "Today is Sat Jan 28 22:46:04 2017" > charmatch(substr(date(),5,7), month.abb) # month.abb は月名の短縮形ベクトル定数
[1] 1
> a = "1,234,567" > gsub(",", "", a)
[1] "1234567"
> as.numeric(gsub(",", "", a)) * 0.9 [1] 1111110
アルファベットの大文字LETTERS、アルファベットの小文字letters、月の名前month.name、月の名前の短縮形month.abb、円周率pi、の５つが文字型定数として定義されている。
> LETTERS [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U" "V" "W" [24] "X" "Y" "Z" > letters [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s" "t" "u" "v" "w" [24] "x" "y" "z" > month.name [1] "January" "February" "March" "April" "May" "June" "July" [8] "August" "September" "October" "November" "December" > month.abb [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec" > pi [1] 3.141593
o
data 関数は日付時間情報を文字列（曜日、月、日、時分秒、年）で返す。
> date() [1] "Thu Feb 16 10:51:39 2017"