# シミュレーションに使える関数を説明します。
# まずは関数sampel()です。 ?sample で使い方をチェックしてください。
?sample
# いつものように英語ですね。ちょっと読む気がしない？ しかし、少しは見てみましょう。
# サイコロの目についてシミュレーションしてみます。
# サイコロを４回振った場合が..
sample(1:6, 4, replace = TRUE)
## [1] 4 1 4 2
# 引数replace = TRUEは、同じ数の繰り返しを許す場合はTRUEになります。
# サイコロは当然、同じ数字が出る可能性があるので、TRUEを設定しました。
# 引数replace を指定しないときは、デフォルト値としてreplace = FALSE、つまり同じ値を繰り返さない、を指定したことになります。
sample(1:20, 10)
## [1] 2 4 11 17 19 7 20 15 3 9
# １から２０までの整数のうち、１０個の整数をサンプリングした結果が表示されました。
LETTERS #大文字のアルファベット
## [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q"
## [18] "R" "S" "T" "U" "V" "W" "X" "Y" "Z"
letters #小文字のアルファベット
## [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q"
## [18] "r" "s" "t" "u" "v" "w" "x" "y" "z"
# LETTERS を使ってサンプリングします。
sample(LETTERS)
## [1] "R" "A" "X" "E" "L" "P" "W" "V" "G" "D" "F" "Y" "N" "U" "B" "O" "I"
## [18] "S" "H" "J" "Q" "M" "Z" "T" "C" "K"
# 次のシミュレーションです。
# コインを投げたときに裏表の出る確率が半々ではなく、0.3と0.7のいかさまコインを考えてみます。
# このコインを１００回投げたときのシミュレーションです。
# その結果を、変数flipsに代入します。
flips <- sample(c(0,1), 100, replace = TRUE, prob = c(0.3, 0.7))
# このシミュレーションでは、上記の同じコードを実行しても、結果は変わってきます。
# それがシミュレーションです。しかし、同じ結果を再現したいときもあります。
# そのときに使う関数が、set.seed()です。
# 具体的には、()内に数値を設定します。具体例は、set.seed(10)、set.seed(2)のように。
# この関数をsample()を実行する前に実行します。
# 裏表の確率は、４つ目の引数 prob を c(0.3, 0.7)と指定しています(prob = c(0.3, 0.7))。
# flipsを表示します。
flips
## [1] 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 0 0 1 1 0 1
## [36] 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 1 1 1 1 1
## [71] 1 1 1 0 1 0 1 1 1 1 0 0 0 0 1 1 1 0 1 1 0 1 0 1 1 1 0 0 1 1
# 実際、このシミュレーションでは表（１）の数がいくつになりますか。
# さあ、どの関数を使えば確認できますか。考えてください。
# 次に、関数rbinom()を説明します。
# Rにおける確率分布は、r***関数（ランダム ramdom）、d***関数（密度 density）、p***関数（確率 probability）、q***関数（４分位 quantile）などがありますが、ここではr***関数のみを説明します。
# 関数rbinom()を使ってみます。
rbinom(1, size = 100, prob = 0.7)
## [1] 58
# これは表が出る確率が0.7 のインチキコインを１００回振った場合の、表の数を求めています。
# 100の大きさを１つの単位として、それを１回だけシミュレーションしています。
# 次に、それと同じコインを１００回振った場合の各々について、表は１、裏は０とし、その結果を変数flipsに代入しています。
# これは１回と１つの単位として、それを１００回シミュレーションしています。
flips2 <- rbinom(100, size = 1, prob = 0.7)
# flips2の表の数を数えます。
sum(flips2)
## [1] 70
# 次の関数の説明に移ります。
?rnorm
# usage: rnorm(n, mean = 0, sd = 1)
# 関数rnorm()は、引数を指定しない場合、標準正規分布（平均０、標準偏差１）の数値を生成します。
rnorm(10)
## [1] 0.4846600 -0.1937766 0.7698791 -0.4465303 0.4117201 1.4139706
## [7] -0.1691271 1.0722778 -0.8466025 0.5594808
#rnorm(10)は１０個の標準正規分布（平均０、標準偏差１）に従う乱数を生成しました。
# 平均を１００、標準偏差を２５とする乱数を１０個生成します。
rnorm(10, 100, 25)
## [1] 81.98693 123.54006 83.66111 84.17968 91.62236 82.97613 93.60896
## [8] 84.28467 68.24840 118.41983
# 次の説明に移ります。
# 平均値を１０とするポアソン分布に従う５つの値からなるグループを、１００グループ生成したい。
# まず、平均値を１０とするポアソン分布に従う５つの値を生成します。
rpois(5, 10)
## [1] 16 6 15 12 7
# 次に、これを１００個生成し、変数my_pois に代入します。
my_pois <- replicate(100, rpois(5, 10))
# my_poisの内容を見てみましょう。
my_pois
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 14 9 7 9 13 9 9 3 9 8 10 6 5
## [2,] 6 13 8 7 10 15 12 22 8 12 8 17 9
## [3,] 5 6 9 7 7 8 17 7 12 8 9 10 8
## [4,] 12 14 11 13 13 8 11 15 10 13 8 6 9
## [5,] 18 10 5 13 14 17 10 13 6 20 8 11 5
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 10 12 14 16 12 14 14 9 6 13 10
## [2,] 16 13 10 8 6 17 13 11 7 7 11
## [3,] 9 5 8 9 17 10 10 11 9 7 12
## [4,] 8 13 9 12 9 5 12 13 10 14 7
## [5,] 13 5 5 15 16 4 8 4 10 10 12
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,] 12 6 10 8 11 13 10 8 5 17 6
## [2,] 8 10 10 14 10 13 12 13 11 9 13
## [3,] 8 11 9 15 5 13 3 8 10 12 14
## [4,] 7 15 9 8 10 11 10 12 10 14 17
## [5,] 9 10 10 11 11 6 7 8 10 9 7
## [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
## [1,] 5 10 13 9 9 13 4 4 8 8 16
## [2,] 6 10 14 13 11 6 4 9 8 8 4
## [3,] 8 7 10 8 8 9 12 11 10 14 13
## [4,] 9 9 12 7 16 5 5 9 7 14 10
## [5,] 9 13 9 13 10 11 11 12 10 8 12
## [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57]
## [1,] 9 12 12 7 14 10 7 10 11 8 9
## [2,] 11 12 14 7 10 12 13 9 15 19 8
## [3,] 10 6 7 7 13 6 7 6 9 9 6
## [4,] 13 8 12 11 7 8 5 7 16 11 6
## [5,] 13 12 9 7 15 4 9 14 7 9 10
## [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [1,] 10 9 16 12 8 8 12 9 6 11 11
## [2,] 8 7 14 12 5 14 6 8 5 9 8
## [3,] 14 4 6 5 12 12 7 10 19 11 6
## [4,] 12 16 11 18 12 10 10 7 10 12 14
## [5,] 11 9 4 5 7 6 14 7 8 6 12
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79]
## [1,] 4 8 7 11 10 9 11 7 10 8 11
## [2,] 11 11 7 14 13 8 11 10 8 10 11
## [3,] 11 10 14 9 12 7 12 6 10 10 12
## [4,] 12 7 11 11 11 8 7 16 15 15 8
## [5,] 10 7 13 11 4 16 7 9 9 11 16
## [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
## [1,] 12 7 4 13 10 10 11 13 7 10 6
## [2,] 7 14 9 11 12 15 11 15 7 9 9
## [3,] 12 10 10 7 7 17 9 9 14 6 7
## [4,] 15 9 7 9 15 4 11 5 9 11 11
## [5,] 7 11 5 5 11 7 9 10 12 11 10
## [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## [1,] 9 13 8 13 6 6 10 12 8 11
## [2,] 11 16 12 8 16 12 4 9 14 15
## [3,] 9 11 12 9 8 8 13 8 15 11
## [4,] 15 13 7 7 10 10 10 14 9 8
## [5,] 9 10 8 7 13 9 7 6 11 15
# 関数replicate()はマトリックスを生成します。
# この場合は、平均値１０のポアソン分布に従う５つの乱数から構成される列を１００生成しました。
# 関数colMeans()を使って各々の列の平均値を求め、変数cmに代入します。
cm <- colMeans(my_pois)
# cmのヒストグラムを描き、列の平均値の分布を確認してみましょう。
hist(cm)
# ほぼ標準的な分布になっています。これが中央極限定理です。
# R にはこの他にも、指数関数rexp()、カイ二乗関数rchisq(), ガンマ関数 rgamma()などいろいろ揃っています。