60. 確率分布と乱数


R ではさまざまな理論分布の確率密度 f(x) ,累積分布 P(X ≦ x) ,確率点 min{ x : P(X ≦ x) > q } ,及びその分布に従う乱数を求めることが出来る.

確率分布と乱数に関する関数の使い方

後に紹介する確率分布に関する関数を使用する場合は,以下のような対応表を用いればよい. この表は,確率分布名が xxx である分布の確率点は “qxxx” で求めることが出来ることを表す.例えば t 分布ならば,この累積分布は pt ,確率点は qt で求めることが出来る.

用途

関数名
(コード名:xxx)

説明

確率密度
(pdf)

dxxx(q)

q は確率点を表す.例えばコード名が norm ならば dnorm(q) となる.

累積分布
(cdf)

pxxx(q)

q は確率点を表す.例えばコード名が norm ならば pnorm(q) となる.

確率点
(quantile)

qxxx(p)

p は確率を表す.例えばコード名が norm ならば qnorm(p) となる.

乱数

rxxx(n)

n は生成する乱数の個数を表す.例えばコード名が norm ならば rnorm(n) となる.

Rに用意されている確率分布

R では以下の理論分布が用意されている.ただし,スチューデント化された分布は qtukey (確率点) と ptukey (累積分布) を求める関数のみ,多項分布は dmultinom(確率密度)と rmultinom(乱数)のみ,誕生日問題の分布(近似解)は pbirthday(一致確率)と qbirthday(一致確率に必要な観測数)のみしか用意されていない.

分布名

コード名

パラメータ

ベータ分布

beta

shape1, shape2, ncp

二項分布

binom

size, prob

コーシー分布

cauchy

location, scale

カイ二乗分布

chisq

df, ncp

指数分布

exp

rate

F分布

f

df1, df2, ncp

ガンマ分布

gamma

shape, scale

幾何分布

geom

prob

超幾何分布

hyper

n, m, k

対数正規分布

lnorm

meanlog, sdlog

ロジスティック分布

logis

location, scale

多項分布

multinom

n, size, prob

負の二項分布

nbinom

size, prob

正規分布

norm

mean, sd

ポアソン分布

pois

lambda

ウィルコクソンの符号付順位和統計量の分布

signrank

m, n

t 分布

t

df, ncp

一様分布

unif

min, max

スチューデント化された分布

tukey

nmeans, df

ワイブル分布

weibull

shape, scale

ウィルコクソンの順位和統計量の分布

wilcox

m, n

MASS パッケージには多変量正規乱数生成関数 mvrnorm() が,mvtnorm パッケージには多変量正規乱数生成関数 rmvt() か用意されている.mvtnorm パッケージには他にも 多変量正規分布に関する関数 dmvnorm(),pmvnorm(),rmvnorm() が用意されている.

確率分布 binom についての書式は以下の通りである.引数 lower.tail には論理値を指定し,TRUE ならば確率は P[X ≦ x] で,FALSE ならば P[X > x] が返される.ちなみに,引数には log, log.p というものもあり,TRUE を指定すると結果は対数を取ったものが返される.

 
 dbinom(x, size, prob)                     # x : ベクトル, prob : 成功の確率
 pbinom(q, size, prob, lower.tail = TRUE)  # q : quantile のベクトル
 qbinom(p, size, prob, lower.tail = TRUE)  # p : 確率ベクトル
 rbinom(n, size, prob)                     # n : 観測数,  size : 試行数

自由度 4 の t 分布において有意水準 0.05 で両側検定を行なう場合は qt(0.025, 4) ,qt(0.975, 4) として確率点を求めれば良いことが分かる.

正規分布 norm について確率密度 f(x) ,その分布に従う乱数を求めるには以下の様にすればよい.

 
 dnorm(0)    # 確率密度 (離散分布の場合は確率関数)

[1] 0.3989423

 rnorm(5)    # 正規分布に従う乱数を 5 個生成

[1] -2.0024918 -0.5996763 -0.3108348 1.2590405 0.3661534 

確率密度のグラフを描くには以下のようにする.

 
 curve(dnorm, -4, 4, type="l")                        # 正規分布
 plot(0:10, dbinom(0:10, 10, 0.5), type="h", lwd=5)   # 二項分布

自由度 (p, n - p - 1) の F 分布の上側 100×α 点は qf(1-α, p, n-p-1 ) で求められるので,例えば自由度 (2, 3) の F 分布の上側 5% 点は以下の様にして求めることが出来る.

 
 qf(0.95, 2, 3)           # qf(0.05, 2, 3, lower=F) でも可

[1] 9.552094 

乱数の再現:set.seed()

疑似乱数は R 起動時に初期化され,毎回異なった乱数が発生される (ちなみに R は脅威のメルセンヌ・ツイスター法により一様乱数を生成している) .乱数を再現したい場合は,乱数の種を関数 set.seed() で指定すればよい.

 
 runif(5)                 # 一様乱数を 5 個生成すると

[1] 0.8039566 0.6593698 0.8120159 0.1279622 0.1115031

 runif(5)                 # 当然毎回違った乱数が得られる

[1] 0.5694816 0.4504077 0.1823374 0.4573758 0.9061499

 set.seed(101); runif(5)  # 乱数の種(seed)を指定

[1] 0.37219838 0.04382482 0.70968402 0.65769040 0.24985572

 set.seed(101); runif(5)  # 乱数の種を同じにすれば乱数を再現することが出来る

[1] 0.37219838 0.04382482 0.70968402 0.65769040 0.24985572

2 次元正規乱数

相関が r だけある 2 次元正規乱数を生成する関数 r2norm() を定義する.ちなみに,MASS パッケージの多変量正規乱数生成関数 mvrnorm() や mvtnorm パッケージの多変量正規乱数生成関数 rmvt() を使う方法もある.

 
 r2norm <- function(n, mu, sigma, rho) {
   tmp <- rnorm(n)
   x   <- mu+sigma*tmp
   y   <- rho*x + sqrt(1-rho^2)*rnorm(n)
   return(data.frame(x=x,y=y))
 }
 mydata <- r2norm(100, 0, 1, 0.5)
 cor(mydata$x,mydata$y)

[1] 0.5756098

3 変量正規乱数

3 変量正規分布に従う乱数を生成する関数 r3norm() を定義する.

 
 r3norm <- function(mu, A, n) {
   U <- svd(A)$u;  V <- svd(A)$v;  D <- diag(sqrt(svd(A)$d))
   B <- U %*% D %*% t(V)                                     # 行列 A の平方根
   w <- c()
   for (i in 1:n) w <- append(w, list(mu + B%*%cbind(rnorm(3))))
   return(w)
 }
 mu <- cbind(c(1,1,1))                                       # 平均ベクトル(縦ベクトル)
 A <- array(c(2,1,1,1,2,1,1,1,2), dim=c(3,3))
 w <- r3norm(mu, A, 2000)

3 変量 t 乱数

まず,以下を準備する.

次に,3 変量楕円 t 分布に従う乱数を生成する関数 met3() を定義する.

 
 met3 <- function(m, mu, V, n) {
  # m  : 自由度, mu : 平均ベクトル, V : 散らばり行列, n  : 乱数の個数
  U  <- svd(V)$u;  V1 <- svd(V)$v;  D  <- diag(sqrt(svd(V)$d))
  B  <- U %*% D %% t(V1)
  w <- c()
  for (i in 1:n) {
     R <- 0
     for (j in 1:m) R <- R + rnorm(1)^2
     w <- append(w, list(mu + B %*% (cbind(rnorm(3))*sqrt(m/R)))) 
   }
   return(w)
 }

この分布の平均ベクトルは μ,分散共分散行列は m/(m-2) * V  (m>2) なので,多変量楕円 t 乱数がうまく出来ているかどうかはこれで確かめることが出来る.

 
 mu <- cbind(c(1,1,1))
 V  <- array(c(2,1,1,1,2,1,1,1,2), dim=c(3,3))
 n  <- 1000
 w  <- met3(5, mu, V, n)
 sm <- 0
 for (i in 1:n) { sm <- sm + w[[i]] }
 sm <- sm/n
 sv <- 0
 for (i in 1:n) { sv <- sv + (w[[i]]-sm) %*% t(w[[i]]-sm) }
 sv <- sv/n

逆関数法による乱数の作り方

上に紹介している分布に無い分布に従う乱数を生成したくなる場合がある.まず,U を一様分布に従う確率変数とし,X = F-1 (U) は

となることより X は分布関数$ F に従う.よって,分布関数の逆関数の引数に一様乱数を入れてやれば欲しい分布の乱数が得られる.

指数分布の分布関数 F(x) = 1-exp(-x) の逆関数は -log(1-x) なので,X = -log(1-U) とすれば X は平均 1 の指数分布に従う.さらに,1-U と U は同分布なので,X = -log(U) の U に一様乱数を入れてやることで平均 1 の指数乱数が得られる  (ただし x > 0 ) .

ラプラス分布 (両側指数分布) の密度関数は f(x) = 1/2 exp(-|x|) だが,この場合は一様乱数 U1,U2 を生成し,U1 が 1/2 以上ならば -log(U2) を,U1 が 1/2 未満ならば log(U2) を採択すればラプラス分布に従う乱数が得られる.

例としてラプラス分布に従う乱数を生成する関数を定義する.

 
 rlaplace <- function(n) {
   u <- log(runif(n))
   v <- ifelse(runif(n)>1/2, 1, -1)
   return(u*v)
 }
 x <- rlaplace(1000)

棄却法による乱数の作り方

逆関数法は理論的には正しい方法なのだが,コンピュータは有限な値しか生成できないので,無限のサポートを持っている分布の裾の方の乱数の精度が悪くなる場合がある (例:逆関数法で生成した指数乱数) ,

そこで以下に棄却法による乱数生成の方法を紹介する.用いる関数を先に説明する.

生成するアルゴリズムはこちら.

  1. u <- runif(1)
  2. v <- 「 g(x) に従う乱数 」
  3. u ≦ h(v) ならば v を乱数として採択する

例として f(x) = 1/2 sin(x)  (0 < x < π) という分布に従う乱数を生成することを考える.このとき

とおくと f(x) ≦ c・g(x) が成り立ち,この場合は h(x) = sin(x) となる.よって以下の関数 myrand(生成する乱数の数) で生成することが出来る.

 
 myrand <- function(x) {
   y <- c();  i <- 1
   while (i <= x) {
     u <- runif(1);  v <- pi*runif(1);  w <- sin(v)
     if (u < w){
       y <- append(y, v);   i <- i+1
     }
   }
   return(y)
 }
 mydata <- myrand(1000)

周辺和を与えたランダムな 2 x 2 分割表

周辺和を与えて 2 x 2 分割表を関数 r2dtable() で作成することが出来る.引数は以下の通り.

返り値は 「行・列和がそれぞれ r,c であるランダムな 2×2 分割表 n 個のリスト」 となる.

 
 x <- r2dtable(2, c(10,10), c(15,5))
 x

[[1]]
     [,1] [,2]
[1,]    8    2
[2,]    7    3
[[2]]
     [,1] [,2]
[1,]    6    4
[2,]    9    1

ランダム抽出・ブートストラップサンプリング

長さ n の与えられたベクトルの要素から長さ m の部分ベクトルをランダムに取り出す場合は,関数 sample() を用いればよい.デフォルトでは replace=FALSE となっている.

 
 n <- 20;  m <- 5
 x <- 1:n                           # 単に n としても良い
 sample(x, m, replace=TRUE)         # 同じ要素が選ばれても良い

[1]  5 18 20 18  5

 sample(x)                          # m = n でランダムな置換

[1]  6 11  9  2  1 14  4  7 15 20 16 12  3  5 18 17 13 10  8 19

 sample(c(0,1), 100, replace=TRUE)  # 100個のベルヌーイ試行

 [1] 1 1 1 1 1 1 0 0 0 1 1 0 1 1 0 0 0 0 0 1 0 0 1 1 1 1 0 1 0 0 1 1 0 1 0 0 0
[38] 0 0 1 0 0 0 0 1 1 0 0 1 1 1 0 1 0 0 1 1 0 1 0 1 0 0 0 1 0 0 1 0 1 1 0 0 1
[75] 1 1 0 0 1 0 0 1 1 1 0 0 1 1 1 1 0 1 1 0 0 0 0 0 1 1

オプションで確率ベクトルを与えると,各要素が選ばれる確率を指定することが出来る (既定値は等確率 1/n ) .

 
 p <- runif(n);  p <- p/sum(p)
 sample(x, m, replace=TRUE, prob=p) 

[1] 17  4 16 12 20