昨日，名古屋に出張しました．新幹線は，コロナ前と比べればまだ空いていますが，海外の旅行者も戻ってきて，席が埋まってきたように感じます．ということで，今回は，Rで乱数を生成するコマンドの整理です．疑似乱数を使って下の図のような，中心極限定理の数値実験もやってみます．
疑似乱数生成コマンドは頭文字 r
Rで確率分布を扱うコマンド表記には基本ルールがあります．乱数を生成するときは，頭文字はrです．コマンドの頭文字は以下のルールになってます．
| 頭文字 | 意味 |
r |
疑似乱数の生成 |
d |
確率密度関数 |
p |
確率分布関数 |
q |
クォンタイル関数 (確率分布関数の逆関数) |
この頭文字の後に，分布を表す文字列が続きます．例えば，正規乱数を生成したいときは，正規分布を意味する文字列normの頭にrを付けて， rnormとなります．その他の分布については，以下を参考にしてください．
| 分布 | Rでの表記 | パラメタ |
| 正規分布 | norm |
平均mean，標準偏差sd |
| 一様分布 | unif |
最小値min，最大値max |
| t分布 | t |
自由度df，非心度ncp |
| F分布 | f |
自由度df1，df2 |
| 指数分布 | exp |
イベント発生率rate |
| 対数正規分布 | lnorm |
対数平均meanlog，対数標準偏差sdlog |
| ワイブル分布 | weibull |
位置location，スケールscale |
| カイ2乗分布 | chisq |
自由度df，非心ncp |
| ガンマ分布 | gamma |
形状shape，スケールscale |
| ベータ分布 | beta |
形状shape1，shape2，非心ncp |
| コーシー分布 | cauchy |
位置location，スケールscale |
| ロジスティック分布 | logis |
位置location，スケールscale |
パッケージを使って利用可能になる分布もあります．パッケージを使いたいときは，Rで一度だけ，
install.packages("パッケージ名")を実行して，パッケージをインストールしてください．実際に使うときは，Rを起動するたびに，
require(パッケージ名)
を実行してください．
パッケージを使って利用可能な分布は以下のようなものがあります (私が知っているものだけです)．
| 分布 | Rでの表記 | パッケージのインストール |
| 安定分布 | stable |
install.packages("stabledist") |
| ラプラス分布 (両側指数) | dexp |
install.packages("nimble") |
| Gumbel分布 | gumbel |
install.packages("gumbel") |
| フレシェ分布 | frechet |
install.packages("VGAM") |
| 逆正規分布 (逆ガウス) | invgauss |
install.packages("statmod") |
| 逆ガンマ分布 | invgamma |
install.packages("invgamma") |
| 逆カイ2乗分布 | invchisq |
install.packages("invgamma") |
| 逆指数分布 | invexp |
install.packages("invgamma") |
乱数を生成するときは，上のRでの表記の前にrをつけてください．必要なパラメタは，
?コマンド
を実行して調べてください．
中心極限定理の数値
乱数を生成して，中心極限定理を体験してみます．中心極限定理は，分散が有限な確率変数をたくさん足すと，正規分布で近似できるようになるというものです．
ここでは，確率変数をの独立なコピーを
として，
で計算される，の分布を描いてみます．
サンプルのRスクリプトです．
# 和をとる変数の数 (Nの値を増やしてみてください)
N <- 1
###################
# 繰り返し回数
N.rep <- 10000
###
x.means <- c()
for(i in 1:N.rep){ # 繰り返し
###################
# 一様乱数(-1, 1)
# (この部分の分布を替えて実験してみてください)
x <- runif(N,-1, 1)
x.means[i] <- sum(x)/sqrt(N)
} # 繰り返し終わり
###################
# xの平均値x.meansの分布
hist(x.means,probability=TRUE,breaks="FD",main=paste("Sum of",N,"variables"),border=4,col=4)
# 正規分布の当てはめ
mu.est <- mean(x.means)
sig.est <- sd(x.means)
curve(dnorm(x,mu.est,sig.est),col=2,add=TRUE,lwd=2,lty=2)
