R言語

概要

RはS言語のオープンソース実装である．バイオインフォをはじめとして統計学全般の分野で多用される．Matlab, Pythonと並び研究分野で多用されるスクリプト言語である．公開されているライブラリは10000を超え，大抵の計算であればライブラリの組み合わせによって計算することが出来る．

インストール

Rのインストールパッケージはここよりダウンロードする（baseと書かれているリンクをクリックする）．Windowsであれば，ダウンロードしたファイルをダブルクリックするとインストーラーが起動する．

試しに，y = x3 を0から10までプロットしてみる．まずしたのコマンドを打つ（最左の>はコンソールであることを表すので，入力はしない）．

> x <- seq(0, 10, by=0.1)

このコマンドは0から10まで0.1刻みで増えるベクトルを生成し，そのベクトルをxに格納する．

> y <- x^3

このコマンドでは，ベクトルxの各要素を3乗してyに格納する．yをxに関してプロットするには次のコマンドを打つ．

> plot(x, y)

上手くいくと以下の図が出力される．点ではなく直線で表示するにはplot(x,y,type="l")とする．

なお，plot関数の詳細を表示させるには，?を用いる．

?plot

ブラウザーが立ち上がり，plot関数の詳細な説明が表示される．一般的に，調べたい関数の前に?を付けると，その関数の詳細が表示される．

言語仕様

Rの簡単な使い方を記す．ここでは本当に基本的な事項のみを記す．詳細は教科書や入門サイトを参照せよ．

変数の代入は以下のコマンドで行う．

> x <- 10

この例では変数xに10を代入する．この変数を用いた計算は，例えば

> x^2

[1] 100

などで可能である（xの二乗，100が出力されている）．

Rではベクトルを扱うことが出来る．ベクトルは以下で定義できる．

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

このコマンドでは，長さ4で要素が1,2,3,4であるベクトルを生成し，それをxに代入している．ベクトルの要素にアクセスするには

> x[1]

[1] 1

のように四角かっこを用いる（この例ではxの第一要素にアクセスしている）．

等間隔のベクトルを生成するには，seqを用いる．例えば，1から10までのベクトルを生成するには

> seq(1,10)

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

長さを指定することも可能である．

> seq(1, 10, length=5) # 1から10までを5等分割したベクトル

[1] 1.00 3.25 5.50 7.75 10.00

差分を指定することも可能である．

> seq(1, 10, by=3) # 1から10までを3間隔のベクトル

[1] 1 4 7 10

Rは行列を扱うことも出来る．変数mに以下の2行2列の行列

[1,  4]

[10, 2]

を定義してみる．この場合，行列を列ごとに並べたベクトルを用い，以下のようにする．

> m <- matrix(c(1, 10, 4, 2),nrow = 2, ncol = 2)

nrow，ncolは行数，列数を表す．行で定義するには

> m <- matrix(c(1, 4, 10, 2),nrow = 2, ncol = 2, byrow = T)

見やすいように，改行しても良い．

m <- matrix(c(

  1, 4, 

  10, 2

),nrow = 2, ncol = 2, byrow = T)

行列mの各値にアクセスするには以下のように行う．

> m[1,2]

[1] 4

上の例では1行2列目の値にアクセスしている．行列に対しても様々な演算が可能である（足し算，引き算，掛け算等）．例えば，行列mの固有値，固有ベクトルを変数vに格納するには

> v <- eigen(m)

> print(v)

とする．結果は

$values

[1]  7.844289 -4.844289

$vectors

           [,1]       [,2]

[1,] -0.5045766 -0.5648066

[2,] -0.8633669  0.8252233

となる．vの型はlist型である．固有値と，その固有値に対応した固有ベクトルが計算されvに格納される（$vectorsのn列目（列ベクトル）がn番目の$valuesに対応している）．固有値にアクセスするには

> v$values

[1]  7.844289 -4.844289

とする．同様に固有ベクトルの場合はv$vectorsとする．

mの逆行列はsolve関数で出来る．

> solve(m)

            [,1]        [,2]

[1,] -0.05263158  0.10526316

[2,]  0.26315789 -0.02631579

制御構文

例えば，1からnまでの和を計算するプログラムを例にとる．このような繰り返しはfor文を用いる．

s <- 0

n <- 100

for (i in 1:n) {

  s <- s + i

}

print(s)

この例では和はsに足されている（上では1から100までの和である）．複数行のスクリプトの実行は下の「スクリプトの実行」を参照せよ．

if文は条件によって分岐するような処理を行うときに用いる．例えば1から100まで，偶数を表示するプログラムは以下のようになる．

n <- 100

for (i in 1:n) {

  if (i %% 2 == 0) {

    print(i)

  }

}

%%の演算子は割り算の余りを計算する（10%%3ならば1）．==は数学での等号の条件である．そのため，iを2で割ったあまりが0の場合，つまり偶数の場合に限り数字を出力する．実行すると，

2

4

...

98

100

のように出力される．

関数

ひとまとまり手続きを関数として定義してみる．for文の例で出てきたプログラムを関数化してみると以下のようになる．

mysum <- function(n) {

  s <- 0

  for (i in 1:n) {

    s <- s + i

  }

  return(s)

}

上のコードでは，mysumという関数を定義し，nまでの和を計算する．計算した結果はreturnによって返される．nは関数を実行するときに引数として与えることが出来る．コンソールで以下のように打ってみると，

> mysum(100)

[1] 5050

となり，正しく1から100までの和が計算される．

スクリプトの実行

一行一行コンソールで実行すると手間がかかるので，複数の命令をまとめて一度に実行してみる．

s <- 0

n <- 100

for (i in 1:n) {

  s <- s + i

}

print(s)

例えば上プログラムをスクリプトとして実行するには，Rのツールバーの「ファイル」から「新しいスクリプト」をクリックする．上のコードをスクリプトのウィンドウに張り付け，張り付けたコードを選択して右クリックし，「カーソル行または選 択中のRコード実行」を選ぶ．そうすることで，コードが一気に実行される．R-Fiddleであれば，スクリプトを張り付けて「Run Code」をクリックするとスクリプトが実行される．

C言語のように特定のファイルに保存して読み込む方法もある．上のコードを適当なファイル名で保存する（ここではtest.Rとする）．Rのツールバーの「ファイル」から「ディレクトリの変更」を選択し，先ほどのtest.Rが保存されたディレクトリを選択する．コンソールに以下のコマンドを打つ．

> source("test.R")

こうすることで，test.Rの中のコードが実行される．Rのスクリプトに対応したエディタはいくつかあるが，Visual Studio Codeなどが使いやすい．もちろんRに搭載されているエディタでも編集できるが，関数や制御構文の色付けがされないため見にくい．

deSolveパッケージ

Rでは大量のパッケージがオープンソースで公開されており，その数は9000を超える．そのため，一般的に使われる計算のほとんどはパッケージ化されており，公開されているパッケージを用いることで簡単に計算を行うことが可能である．課題では微分方程式の数値解法を用いるが，deSolveパッケージを用いる．deSolveパッケージのインストールは以下の要領で行う．

1. Rの上のメニューバーから「パッケージ」→「パッケージのインストール」をクリック
2. ミラーサイトは適当に選ぶ
3. PackagesからdeSolveを選択してＯＫを押す
4. deSolveを使うときは，library(deSolve)コマンドをコンソールで一度実行する

FitzHugh-Nagumo振動子をシミュレートしてみる．

dx/dt = x - x^3 - y;

dy/dt = a x;

初期値をx = 0, y = 1とし，パラメータをa = 1.0として実行するには以下のようなプログラムを書く．時間はt = 0からt = 100まで，刻み幅0.1で実行する．

library(deSolve)

model <- function (time, xs, parms) {

  with(as.list(c(xs, parms)), {

    dxdt <- x - x^3 - y

    dydt <- a * x

    list(c(dxdt, dydt))

  })

}

y <- c(x = 0.0, y = 1.0)

parms <- c(a = 1.0)

times <- seq(0, 100, 0.1)

out <- ode(y, times, model, parms)

plot(out)

実行すると以下のような振動する図が表示される．

tとxだけを表示するには

plot(out[ ,c("time","x")],type="l")

とする．xとyであれば

plot(out[,c("x","y")],type="l")

とする．

Links

• ドットインストール
  • 3分動画*13回で学べるサイト
• R-Tips
  • Rの基礎から発展が学べるサイト
• R-fiddle
  • オンラインでRを利用できるサイト
• biostatistics
  • Rの初歩から発展まで学習できるサイト
• Rの初歩
  • Rの初歩を学習できるサイト