Rで乱数を使う Rでシミュレーションする、Rで関数を作る

• 前の記事で、確率的に出力が変わる場合も扱った。そこではsample()関数を用いたが、もう少し一般的に乱数を用いてみる
• 一様乱数
  • 0から1までの一様乱数は0から1までの値のいずれかを取り、どの値も、等確率で得られるような乱数

# 発生させる乱数の数
N<-100000
R<-runif(N)
hist(R)

• a1,a2,a3をp1:p2:p3の確率で発生させる

as<-c("a1","a2","a3")
p1<-10
p2<-5
p3<-2
p<-c(p1,p2,p3)
p<-p/sum(p)
p
N<-10000
R<-sample(as,N,replace=TRUE,prob=p)
table(R)/N

> p
[1] 0.5882353 0.2941176 0.1176471
> N<-10000
> R<-sample(as,N,replace=TRUE,prob=p)
> table(R)/N
R
    a1     a2     a3 
0.5848 0.2957 0.1195 

• • 一様乱数を使って同じことをしてみよう
    • 一様乱数を発生させ、その値が0-p1のときはa1,p1-p2のときはa2,p2-1のときはa3を取ることにする

cumsum.p<-cumsum(p)
cumsum.p
R<-runif(N)
R2<-rep(length(p),N)
for(i in (length(cumsum.p)-1):1){
  tmp<-which(R<cumsum.p[i])
  R2[tmp]<-i
}
table(R2)/N

> table(R2)/N
R2
    1     2     3 
0.588 0.296 0.116 

• 正規分布に従う乱数

N<-10000
m<-10
sd<-5
R<-rnorm(N,mean=m,sd=sd)
hist(R)

• • これも、同様に一様乱数を使って発生させてみる

N<-10000
m<-10
sd<-5
R<-runif(N)
# 正規分布のクオンタイル値(0-1の値)に対応する値を取ればよい
R2<-qnorm(R,mean=m,sd=sd)
hist(R2)

• • Rで扱っている確率分布の一覧(こちら)を確認すれば、それに対応する乱数が発生できる
• 複数のカテゴリにいろいろな生起確率を割り当てたい
  • ２カテゴリの場合には、片方のカテゴリの生起確率を0から1の一様乱数を使うと、「あらゆる生起確率を均等に試せ」る
  • ３以上のカテゴリになると別の方法が必要になる

install.packages("MCMCpack")
library(MCMCpack)
k<-3
N<-1000
R<-rdirichlet(N,rep(1,k))
plot(as.data.frame(R))