前回、乱数の意味と使い方を解説したのに続いて、今回はRで実際に乱数の発生させてみます。
・正規乱数の生成
正規分布に従う乱数を生成するときはrnorm()という関数を使います。デフォルトでは、平均値0、標準偏差１の標準正規分布に従った乱数を作ります。
n <- rnorm(100)
標本の大きさ、平均、標準偏差を弄りたいときは下のようにします。サイズ、meanが平均、sdで標準偏差を設定できます。
標本数1000、平均10、標準偏差５の正規乱数ならば下のように入力します。
n<-rnorm(1000,mean=10,sd=5)
hist(n)
なんか釣り鐘型とはちょっと違う気もしますが、これは計算ではなく、ちゃんと現実に近い乱数を発生させているということの証明ということにしておきましょう･･･。実際1000回くらいでは綺麗な正規分布にならないことは結構あります。
・一様乱数の生成
一様乱数とはどの値も、均等な確率で出てくる乱数を指します。代表的なのは、1から6のどの目も1/6の確率で出るサイコロの出目です。一様乱数の発生には、sample()関数かrunif()関数を使います。もしsample関数でサイコロを10000回転がすシュミレーションをするならば、下のように入力します。
dice<-1:6
R<-sample(dice,10000,replace=TRUE)
hist(R)
・べき分布に従う乱数生成
こいつはちょっとマニアックなのですが、べき分布は、株価の変動率などの経済現象によく出てくる分布なので、株式投資への統計学の応用を掲げている当ブログとしては見逃せない分布ということもあり、作りかたを紹介しておきます。
べき分布の乱数生成の関数は一様乱数や正規乱数のようにRでデフォルトに設定されていないのでfunction()を使って設定する必要があるらしいです。
コマンドは自分ではよく分からなかったので
http://alcuin.hatenablog.com/entry/2013/03/25/044624から引っ張ってきてます。
>rpowerlaw <- function (n, dist.power = 0, x0 = 0, x1 = 1) {
sapply(runif(n),
function (x.unif) (((x1 ^(dist.power + 1) - x0 ^(dist.power + 1))
* x.unif + x0 ^(dist.power + 1))
^(1 / (dist.power + 1)))
)
}
10000個の、0..1の範囲でべき分布に従った乱数を発生させる場合は下のように入力します。0,1の前の7の意味がよく分からないのですが、7以外だとヒストグラムがメチャクチャになります。べき分布については私も勉強不足なので、これから勉強していくつもりです。
>hist(rpowerlaw(10000,dist.power = 7, 0, 1))
関連記事
