Rでシミュレーション練習：乱数生成
このnoteでは、統計ソフトウェアRを利用して簡単なシミュレーションをする方向けに、私が学習した基礎的なRの内容を解説します。
1. Rについて
Rで利用するオブジェクト変数は、数値型、文字型、論理型、（因子型）の4種類に分かれます。ここで、オブジェクト変数をRが扱う数値、計算結果、関数などの名前を付けられた対象のことを指します。
例えば、Rで以下のオブジェクトaに数値型、オブジェクトbに文字型、オブジェクトcに論理型を格納しました。ここで<-の記号は”代入”を表します。例えば、以下のプログラムでは、aというオブジェクトに数値型変数として5を代入しています。=記号を利用しても代入を表しますが、私は<-で統一して利用しています。また、PCで<-と入力するのは手間なので、<-のショートカットキーとしてAlt＋-を利用するのがが便利です。
a <- 5
b <- "stat"
c <- TRUEこのa,b,cがそれぞれ数値型、文字型、論理型のオブジェクトであるかを確認する方法は下記となります。class関数は、格納されているオブジェクトを明らかにしてくれます。このプログラムの結果はそれぞれ、[1]"numeric"、[1]”character”、[1]"logical"と表示されます。よってa,b,cはそれぞれ数値型、文字型、論理型であることが確認できます。
class(a)
class(b)
class(c)> class(a)
[1] "numeric"
> class(b)
[1] "character"
> class(c)
[1] "logical"因子型は、カテゴリカルなデータを扱うときなどに使われます。例えば、以下のようなベクトルデータに対して、1を男性、2を女性として入力されたデータがあったとします。class関数によってdは数値型であると分かります。
d <- c(1,2,1,1)> class(d)
[1] "numeric"これを因子型にするには、as.factor関数を利用します。as.factorのasは変換を意味しており、as.factor関数以外にも、後述するベクトル形式やdata.frame型に変換できるas.vector関数やas.data.frameも存在します。以下の結果より、数値型が因子型に変換されていることが分かります。
x <- as.factor(d)
print(x)
levels(x)
class(x)> print(x)
[1] 1 2 1 1
Levels: 1 2
> levels(x)
[1] "1" "2"
> class(x)
[1] "factor"2. ベクトル表示と乱数生成
Rでシミュレーション実験をする際は、生成したデータをベクトル型や行列型（＋data.frame型）で保持することが多いと思います。例えば、標準正規分布に従う確率変数100個のベクトルを生成するRプログラミングはrnorm関数を利用して以下のように書けます。
x <- rnorm(100)rnorm関数のrnorm(n,mean=0,sd=1)について、n,mean,sdが引数で、デフォルトでmean=0,sd=1となりますので、nだけ指定することで標準正規分布に従う乱数をn個作成できます。rnormのrはrandomだと思いますので、randomにnormal（正規分布）から乱数を生成する関数と覚えておくとよいと思います。そう考えると、正規分布以外にも確率分布の英語表記を知っていれば、1変数の確率密度関数の乱数生成は以下を参考に容易に行えます。
#代表的な離散確率分布の乱数生成
#二項分布
rbinm(n,x,p)
#多項分布
rmultinom(n,size,prob)
#ポアソン分布
rpois(n,lambda)
#幾何分布
rgeom(n,prob)
#負の二項分布
rbinom(n,size,prob)
#代表的な連続分布の乱数生成
#正規分布
rnorm(n,mean=0,sd=1)
#一様分布
runif(n,min=0,max=1)
#指数分布
rexp(n,rate=1)
#ベータ分布
rbeta(n,shape1,shape2,ncp = 0)
#ガンマ分布
rgamma(n,shape,rate=1,scale=1/rate)ここで、パラメータの引数で注意が必要な確率分布はガンマ分布です。教科書によってガンマ分布の確率密度関数の定義は異なり、Rではガンマ分布の確率密度関数を下記で定義しています。
$${f(x)=\frac{1}{\sigma\Gamma(\alpha)}x^{a-1}e^{-x/\sigma}}$$
よってガンマ分布の期待値は$${E[X]=\alpha\sigma}$$で分散は$${V[X]=\alpha\sigma^2}$$となるので、引数を指定する際には注意が必要です。
3. 行列表示とリスト表示
行列について、行列はmatrix関数を利用して、以下のように作成できます。matrix関数には、業悦にするデータと行数と列数を指定します。データを行から入力するか列から入力するかをbyrowで指定します。デフォルトでは行から入力されるようです。
matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE,
dimnames = NULL)リスト表示について、リスト表示は様々なオブジェクトを詰め込める箱と思うとよいかと思います。例えば、シミュレーション条件を記した変数（ベクトル）とシミュレーション結果（ベクトル）、利用したデータ（行列 ）をまとめて1つの変数のようにlist型で保持することができます。例えば、以下のようなコードになります。
a <- 5
b <- "stat"
c <- TRUE
x <- rnorm(100)
list<-list(a,b,c,x)list関数を用いて、a,b,c,dをまとめてlistとして格納しています。こうすることで、これまでの結果を1つのオブジェクトとして格納することができます。リストの要素にアクセスするためには、リストの何番目の要素にアクセスしたいかを考えて、以下のように記述します。こうすることによって、解析結果等の情報を全てリスト型として保存することができます。
list[[1]]
list[[2]]
list[[3]]
list[[4]]> list[[1]]
[1] 5
> list[[1]]
[1] 5
> list[[2]]
[1] "stat"
> list[[3]]
[1] TRUE
> list[[4]]
[1] 0.27329690 -0.93859939 -1.23969852 0.84926657 -0.80390340 0.38316615 -1.16824447 -0.15568960
[9] 1.92199925 -0.58325587 0.61551326 -0.03324432 -0.29845277 -0.98317322 0.80757980 0.50832968
[17] 0.53516186 -1.05780987 0.06063103 -1.46826970 0.86956256 -0.13838560 -0.43793469 -2.17091301
[25] 1.27977024 -0.26948510 0.47206434 0.72882642 0.03432766 0.24426625 0.48268525 -0.31862249
[33] -0.12547027 -0.23289232 -1.00496878 0.95015900 0.19653412 0.17546243 -0.80079231 -0.77258472
[41] -1.57811873 -1.32666304 0.48547020 -0.60371080 -0.26183886 -0.94348008 -0.41967547 -1.66710785
[49] 0.96228424 -0.12048306 0.75649305 -0.69782709 0.34020032 1.44483900 0.89300996 0.39468935
[57] 0.45782261 -0.19787514 1.65051958 -0.03660476 -0.43882793 2.78720674 -0.88047025 -0.47640251
[65] 0.30847152 -0.02017465 -0.04159986 -0.85958388 0.59576608 -0.17709183 -2.28235200 0.68484393
[73] 2.36119822 -0.58981171 1.29408104 0.53322570 1.50142459 -0.13618502 -0.42911151 0.16976167
[81] 0.66311451 0.61210340 -0.63380295 -0.02263158 -0.08955137 -0.64158978 0.57738415 -2.08141977
[89] -0.15137456 0.47633514 -0.08859354 1.18939749 0.92771688 0.58765277 0.76750430 -1.74553874
[97] 0.78970091 0.86203201 -0.02194766 1.07919386もしくは、リスト表示の各要素に変数名を割り当てて、以下のように格納することも可能です。ここで、$記号はlistというオブジェクトの要素の変数にアクセスすることを示しています。Rでは$記号は頻繁に利用します。
list<-list(num=a,char=b,logi=c,data=x)list$num
list$char
list$logi
list$data3. data.frame型
シミュレーション研究を目的とせず、一般的なデータ分析に対しては、data.frame型のデータセットが一般的かと思います。データフレーム型のメリットは私はイマイチわかりませんが、Excel表みたくデータを保持できることが魅力なのではないかと思っています。Matrixでは、要素に文字変数を格納できないので、そういう点ではdata.frame型は重宝されるのかもしれません。実データの分析においては、data.frame型で行に対象者ID、列に変数のn×p行列の形で保持することが多いと思います。data.frame型において、列ベクトルがp個あるとみなせるので、ベクトルが列でつながっているデータフレームと解釈することもでき、下記のようにベクトルを繋げて作成することができます。
id <- 1:5
x1 <- rpois(5,lambda = 2)
x2 <- runif(5,min=0,max=1)
x3 <- c("男性","女性","男性","男性","女性")
y <- rnorm(5)
data <- data.frame(id,x1,x2,x3,y)class(data)> class(data)
[1] "data.frame"データフレーム型の情報を表示するための関数として以下のようなものがあります。data.frame型が想定している通りに作成できているかを確認するのに有用だと思います。
#data .frameの次元を返す
dim(data)
#data .frameの行数を返す
nrow(data)
#data .frameの列数を返す
ncol(data)
#data .frameの列のラベルをベクトルで返す
colnames(data)
#data .frameの行のラベルをベクトルで返す
rownames(data)
#列名
names(data)
#data .frameの構造を表示する
str(data)
#データフレームの冒頭が表示される
head(data)
#データフレームの2行を表示させる
head(data,2)> #data .frameの次元を返す
> dim(data)
[1] 5 5
>
> #data .frameの行数を返す
> nrow(data)
[1] 5
>
> #data .frameの列数を返す
> ncol(data)
[1] 5
>
> #data .frameの列のラベルをベクトルで返す
> colnames(data)
[1] "id" "x1" "x2" "x3" "y"
>
> #data .frameの行のラベルをベクトルで返す
> rownames(data)
[1] "1" "2" "3" "4" "5"
>
> #列名
> names(data)
[1] "id" "x1" "x2" "x3" "y"
>
> #data .frameの構造を表示する
> str(data)
'data.frame': 5 obs. of 5 variables:
$ id: int 1 2 3 4 5
$ x1: int 3 0 2 1 1
$ x2: num 0.0983 0.5388 0.2338 0.0933 0.8549
$ x3: chr "男性" "女性" "男性" "男性" ...
$ y : num 0.163 0.152 1.156 -0.618 -0.864
>
> #データフレームの冒頭が表示される
> head(data)
id x1 x2 x3 y
1 1 3 0.09827769 男性 0.1625317
2 2 0 0.53883258 女性 0.1517656
3 3 2 0.23376628 男性 1.1560645
4 4 1 0.09327212 男性 -0.6175629
5 5 1 0.85492076 女性 -0.8643795
>
> #データフレームの2行を表示させる
> head(data,2)
id x1 x2 x3 y
1 1 3 0.09827769 男性 0.1625317
2 2 0 0.53883258 女性 0.1517656データフレームの要素を抽出する場合は、list型と類似した考え方で、行や列の情報を指定したり、変数名を指定したりすることで抽出します。この時に、抽出の方法によっては、data.frame型が数値型に変換されます。
#データの2 ,3列の値だけ取り出す
data[c(2,3)]
class(data[c(2,3)])
#データの2 ,3列（x1,x2）の値だけ取り出す
data[c("x1","x2")]
class(data[c("x1","x2")])
#dataのid列を取り出す
data$id
class(data$id)
#dataのid列の4番目の値を取り出す
data$id[4]
class(data$id[4])> data[c(2,3)]
x1 x2
1 3 0.09827769
2 0 0.53883258
3 2 0.23376628
4 1 0.09327212
5 1 0.85492076
> class(data[c(2,3)])
[1] "data.frame"
>
> #データの2 ,3列（x1,x2）の値だけ取り出す
> data[c("x1","x2")]
x1 x2
1 3 0.09827769
2 0 0.53883258
3 2 0.23376628
4 1 0.09327212
5 1 0.85492076
> class(data[c("x1","x2")])
[1] "data.frame"
>
> #dataのid列を取り出す
> data$id
[1] 1 2 3 4 5
> class(data$id)
[1] "integer"
>
> #dataのid列の4番目の値を取り出す
> data$id[4]
[1] 4
> class(data$id[4])
[1] "integer"4. まとめ
本noteでは、Rでシミュレーションを実行するためのRの基礎事項として、Rでのデータ表示の方法について解説しました。本記事が皆様の勉強に少しでも参考になれば嬉しく思います。