乱数Tips大全の編集

SIZE(25){COLOR(red){乱数 Tips 大全}}

//間瀬(2003.08.17)

#contents
~

*Rにおける(疑似)乱数発生システムの概要 [#y9d2907f]

R は確率分布に関する豊富な関数群を持ち，従来統計学の利用で不可欠であった各種確率分布に関する数表は、よほど特殊なものを除き，もはや必要が無い．R では例えば正規分布は norm という名前を持ち、 dnorm はその密度関数、 pnorm はその分布関数， qnorm はそのクォンタイル関数，更に rnorm は正規(疑似)乱数の生成関数を表す，等の一貫した記法を使う．

*R の確率分布関数の使用例 [#w174b58f]

詳しい説明は help(rnorm) 又は ?rnorm  で得られる．

> dnorm(1, mean=2, sd=3) # N(2,9) の密度関数 f(1) の値
 [1] 0.1257944
 > pnorm(1, mean=2, sd=3) # 分布関数 F(1)
 [1] 0.3694413
 > qnorm(0.05, mean=2, sd=3) # クオンタイル関数 F(-2.934561)=0.05
 [1] -2.934561
 > rnorm(3, mean=2, sd=3)  # 疑似乱数を 3 つ生成
 [1]  5.287193 -0.211181  4.648124
 > rnorm(3, m=2, s=3)  # 省略形(引数名は他と一意に区別できる先頭部分だけでOK)
 > rnorm(3, 2, 2)           # 省略形

*R の基本確率関数一覧 [#t313e777]

*Rで使える疑似乱数発生法(R 1.7.0)以降 [#z7f501b5]

-.Random.seed は乱数のシードを含む整数ベクトル。保管したり、元に戻したり出来るが、ユーザーが変更すべきではない。
-RNGkind は RNG の種類を問い合わせたり、変更するためのより扱いやすいインタフェイスである。
-RNGversion 以前の R のバージョンの乱数発生法を設定するために使うことが出来る(一貫性のために)。
-set.seed は乱数種を指定するためのお勧めの方法である。

用法:
 .Random.seed <- c(rng.kind, n1, n2, ...)
 save.seed <- .Random.seed
 
 RNGkind(kind = NULL, normal.kind = NULL)
 RNGversion(vstr)
 set.seed(seed, kind = NULL)

引数:
 kind: 文字、もしくは NULL。もし kind が文字列ならば、R の RNG を指定した方法に 
         設定する。もし NULL ならば現在の手法を返す。"default" にすると現在の既定
         手法に戻す。  
 normal.kind:  文字、もしくは NULL。もし kind が文字列ならば、R の正規乱数発生
        法を指定した方法に設定する。もし NULL ならば現在の手法を返す。"default" 
        にすると現在の既定手法に戻す。  
 seed: 単一の値。一つの整数と解釈される。
 vstr: バージョン番号を含む文字列、例えば "1.6.2"
 rng.kind:  上の kind に対する 0:k 中のコード
 n1, n2, ...:  整数。rng.kind に依存する。詳細は各手法を参照

-"Mersenne-Twister" R 1.7.0 以降の既定の疑似乱数発生法。
Matsumoto and Nishimura (1998) の  twisted GFSR 法。周期 2^19937 - 1(=4.315 x 10^(6001))。すべての周期に関し 623 次元空間に均一分布。乱数種は32ビット整数の624次元ベクトル＋その中の現在位置を示す整数
-"Wichmann-Hill"'  乱数種, `.Random.seed[-1] == r[1:3]' は長さ３の整数ベクトルであり、
各 `r[i]' は `1:(p[i] - 1)' 中にある。ここで `p' は長さ３の素数のベクトル 
`p =(30269, 30307, 30323)' である。Wichmann-Hill 生成規則は長さ
6.9536e12 を持つ(= `prod(p-1)/4', 原論文を修正した Applied Statistics (1984) 33, 123 を見よ)。
-"Marsaglia-Multicarry"' キャリーつきの乗算 RNG を使い、George Marsaglia が彼の
メイリングリスト `sci.stat.math' への投稿記事で紹介した。2^60 以上の周期を持ち、
すべての検査をパスする (Marsaglia による)。乱数種は二つの整数である (あらゆる値が許される)。
-"Super-Duper"' Marsaglia の有名な７０年代の Super-Duper 法。これは MTUPLE test of the Diehard battery をパスしないオリジナル版である。多くの初期値に対し周期 4.6*10^18 を持つ。乱数種は二つの整数である(最初の種は任意、二番目は奇数でなければならない)。
二つの種はそれぞれ Tausworthe と congruence 倍長整数である。S の  `.Random.seed[1:12]' への一対一対応が可能であるが、好評はしないし、この生成法は S-PLUS の最近のバージョンとは全く同じではない。
-"Knuth-TAOCP" Knuth (1997) による。減算を伴う遅延つきフィボナッチ数列を用いる 
GFSR 法。つまり、使用される再帰式は X[j] = (X[j-100] - X[j-37]) mod 2^30 であり、種は最近の 100 個の数の集合である (実際は 101 個の数であり、最後のものはバッファのサイクリックなシフトである)。周期は約  2^129。
-"Knuth-TAOCP-2002" 以前のバージョンとは上位互換ではない 2002 年のバージョンであり、種からの GFSR の初期化が変更された。R は連続する種の選択を認めないという報告済みの弱点を持っており、種をすでに撹拌していた。
-"user-supplied" ユーザー提供の発生法を使用。詳細は `Random.user' を見よ。

*正規疑似乱数発生法 `normal.kind' [#w86e8ffc]
-"Inversion"' (現在の既定手法)
- "Kinderman-Ramage" (R 1.7.0 以前の既定手法、近似誤差を含み使うべきではない)
- "Buggy Kinderman-Ramage"
- "Ahrens-Dieter"
- "Box-Muller"
-"user-supplied"

`set.seed' はその単独の整数引数を用い、要求されただけの種を設定する。これは最小の整数引数を指定することにより全く異なった種を得、同時により複雑な手法（特に`"Mersenne-Twister"' と `"Knuth-TAOCP"'）に適正な乱数種集合を得る、簡単な方法を目指している。

値:

-`.Random.seed' は整数ベクトルであり、その最初の要素は RNG と正規乱数発生法の種類をコーディングしている。最後の二つの十進桁数は、k を利用可能な RNG の数として 0:(k-1)  中にある。この 0-99 の値は正規乱数の型を表す。
-背後にある C コードでは `.Random.seed[-1]' は `unsigned' であり、従って R では `.Random.seed[-1]' は負であっても良い。
-`RNGkind' は、呼び出しの前に使用中の RNG と正規乱数発生法の二文字の文字ベクトルを返す。もしどちらの引数も 'NULL' でなければコンソールには表示されない。
-`RNGversion' は同じ情報を返す。
-`set.seed' は `NULL' を返すが、コンソールには現れない。

注意:
-最初はシードは無い。必要な時現在の時間からつくり出される。したがって、既定では異なったセッションは異なったシミュレーション結果を与える。
-.Random.seed' は、少なくともシステムの発生法に対する、一様疑似乱数に対するシードのセットを保存する。他の発生法の状態を保存せず、特に、Box-Muller 正規疑似乱数の状態を保存しない。もし、後で状態を再生したければ .Rnadom.seed を設定するのではなく .set.seed を呼び出すべきである。

例
     runif(1); .Random.seed; runif(1); .Random.seed
     ## もしシードがなければ新しいものが"ランダム"につくり出される。
     rm(.Random.seed); runif(1); .Random.seed
     
     RNGkind("Wich") # (kind に対する部分文字列マッチング)
          
     ## 以下は runif(.) が  Wichmann-Hill 法に対しどのように動作するかを R 関数のみで示す
     p.WH <- c(30269, 30307, 30323)
     a.WH <- c(  171,   172,   170)
     next.WHseed <- function(i.seed = .Random.seed[-1])
       { (a.WH * i.seed) %% p.WH }
     my.runif1 <- function(i.seed = .Random.seed)
       { ns <- next.WHseed(i.seed[-1]); sum(ns / p.WH) %% 1 }
     rs <- .Random.seed
     (WHs <- next.WHseed(rs[-1]))
     u <- runif(1)
     stopifnot(
      next.WHseed(rs[-1]) == .Random.seed[-1],
      all.equal(u, my.runif1(rs))
     )
     
     ## ----
     .Random.seed
     ok <- RNGkind()
     RNGkind("Super")  # "Super-Duper" にマッチ
     RNGkind()
     .Random.seed # Super-Duper に対する新しいシード
       
     ## リセット:
     RNGkind(ok[1])

*疑似乱数発生法による違いを示す例 [#rcd56d4a]

2種類の疑似乱数発生法の比較。引き続く3組みの一様疑似乱数を座標とする3次元単位立方体中の1千万個の点の断面(実際は薄片)を座標面に射影する。MM法等の乗算合同法の癖を示す例(学芸大学の森さんの "発見" に基づく例)。

Mersenne-Twister法使用(R 1.7.0 以降の既定手法)
 random1 <- function () {
   old.par <- par(no.readonly = TRUE); on.exit(par(old.par))
   png("MTrand.png")
   plot(c(0,0,1,1), c(0,1,0,1), main="", xlab="", ylab="", type="n") # まず枠だけ描く
   RNGkind(kind = "Mersenne-Twister")
   for (i in 1:10^7) {
     x <- runif(3)     # 一千万個の点の座標を一様疑似乱数で発生
     # x 座標が 1/1000 以下のものだけ y,z 座標に射影
     if (x[1] < 0.001) points(x[2], x[3], pch = ".")
   }
   dev.off()
 }
#ref(乱数Tips大全/MTrand.png, ,left)

Marsaglia-Multicarry 法使用(R 1.7.0 以前の既定手法)。すべての点が平行な面の上に乗っていることがわかる！
 random2 <- function () {
   old.par <- par(no.readonly = TRUE); on.exit(par(old.par))
   png("MMrand.png")
   plot(c(0,0,1,1), c(0,1,0,1), main="", xlab="", ylab="", type="n")
   RNGkind(kind = "Marsaglia-Multicarry")
   for (i in 1:10^7) {    # 一千万個の点の座標を一様疑似乱数で発生
     x <- runif(3)
     # x 座標が 1/1000 以下のものだけ y,z 座標に射影
     if (x[1] < 0.001) points(x[2], x[3], pch = ".")
   }
   dev.off()
 }
#ref(乱数Tips大全/MMrand.png, ,left)

* ランダム抽出と並べ変え [#r225ebc9]

- 書式 COLOR(red){sample(x, size, replace = FALSE, prob = NULL)}
- 引数
-- COLOR(magenta){x} ベクトル(数値、複素数、文字列、論理値)
-- COLOR(magenta){size}  選び出される個数を表す正整数
-- COLOR(magenta){replace}  論理値。復元抽出を行なうか？
-- COLOR(magenta){prob}  ベクトルから抽出する際の重み確率
- 注意
-- x が長さ 1 なら 1:x から抽出
-- size の既定値は length(x)。sample(x) はランダムな並べ換えを意味する
-- prob 比率を意味し、総和が 1 でなくても良いが、非負で、すべてがゼロであってはいけない。もし replace = FALSE ならば、これらの確率は sequential に適用される(つまり、次の項目の選択確率は残った項目に対する比率に比例して選ばれる)

> sample(1:10)  # ランダムな並べ換え
  [1]  4  1  8  3 10  9  5  2  6  7
 > sample(1:10)
  [1]  9  3  6  4  7  8 10  1  2  5
 > sample(1:10, 5)  # ５個を非復元抽出
 [1] 10  9  7  2  5
 > sample(1:10, replace=TRUE)  # 復元抽出
  [1] 10  2  6  7  4  3  1  9  1  7
 > sample(1:10, 5, replace=TRUE)
 [1]  7  4 10 10  7

# 成功確率 0.8 の長さ百のベルヌイ試行データを生成
 > sample(c(0,1), 100,  replace = TRUE, prob = c(0.2, 0.8))
   [1] 1 1 1 0 1 1 0 1 1 1 0 1 1 1 1 0 1 0 0 1 1 0 1 1 1 0 1 1 1 0 1 0 0 1 0 1 0
 [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 0 1 1 1 1 1 1 0 1 0 1 0 0 0 0
 [75] 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 0 1 1 1 0 1 1 0 1 0 1

# 単純無作為抽出法のパラドックス
 # 10万個の母集団の平均値が100個の単純無作為標本の平均で十分推定可能なら
 # 100万個の母集団の平均値も100個の単純無作為標本の平均で十分推定可能!
 > x <- rnorm(1000000)
 > y <- numeric(10000)
 > for (i in 1:10000) y[i] <- mean(sample(x,100))
 > sd(y)   # 百万個から百個を単純無作為抽出した時の推定平均値の標準偏差
 [1] 0.1003183
 > x <- rnorm(100000)
 > for (i in 1:10000) y[i] <- mean(sample(x,100))
 > sd(y)  # 10万個から百個を単純無作為抽出した時の推定平均値の標準偏差
 [1] 0.1000383   # 殆んど一緒

//* 周辺和を与えたランダムな ２x２分割表 
// 2×2に限らないので
* 与えた周辺和を持つランダムな二次元分割表 [#hac9d045]
- 書式 COLOR(red){r2dtable(n, r, c)}
- 引数
-- COLOR(magenta){n}  生成する分割表の数
-- COLOR(magenta){r}  行和を与えるベクトル
-- COLOR(magenta){c}  列和を与えるベクトル (c、r の総和は一致する必要)
- 返り値  行・列和がそれぞれ r,c であるランダムな二次元分割表 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

> r2dtable(1,c(3,6,4), c(4,5,4))
 [[1]]
      [,1] [,2] [,3]
 [1,]    0    2    1
 [2,]    3    1    2
 [3,]    1    2    1

n 行 m 列の二次元分割表を対象にする検定で，計算量が大きい場合に[[モンテカルロ法により近似的な解を求める:http://aoki2.si.gunma-u.ac.jp/R/monte-carlo-u-test.html]]ときに使える。

* 乱数の再現 set.seed() 関数 (2003.12.27) [#e2a1ab34]
(疑似)乱数は R 起動時に適当に初期化(システム時間を利用？)され、それ以降毎回異なった乱数が発生される。もし数値計算の再現性が問題になる時は、乱数種を set.seed(整数値) 関数で指定する。与えられた整数値は内部的に各乱数生成法に適切な真の乱数種(それを知る必要はほとんど無いはず)に変換される。

> runif(5)
 [1] 0.8797957 0.7068747 0.7319726 0.9316344 0.4551206
 > runif(5) # 当然毎回違った乱数が得られる
 [1] 0.59031973 0.82043609 0.22411848 0.41166683 0.03861056
 > set.seed(101); runif(5)  # 乱数種を指定
 [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

* 逆関数法による乱数の作り方(2004.01.26) [#j3f49d96]

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

-P(X ≦ x) = P(F^{-1} (U)≦x) = P(U≦F(x)) = F(x)

となることより 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) を採択すればラプラス分布に従う乱数が得られる．

**両側指数分布の乱数を発生させる(2005.7.23) [#l663e7eb]

myrand <- function(n) {
   u <- ifelse(runif(n)>1/2, 1, -1)
   v <- log(runif(n))
   return(u*v)
 }
 mydata <- myrand(100000)                                                # 乱数を生成
 plot(density(mydata), xlim=c(-7, 7), ylim=c(0, 0.5), col="red", ann=F)  # 乱数のグラフ
 laplace <- function(x) 1/2*exp(-abs(x))                                 # 密度関数
 par(new=T)
 curve(laplace, xlim=c(-7, 7), ylim=c(0, 0.5))                           # 本当のグラフ

#ref(rand01.gif, center)

「平均と分散をいろいろ変えてみたい！」とおっしゃる場合は，逆関数法に従って上の関数定義を修正して下さい．

* 棄却法による乱数の作り方　(2004.01.15) [#j72f9ad2]

逆関数法は理論的には正しい方法なのだが，コンピュータは有限な値しか生成できないので，無限のサポートを持っている分布の裾の方の乱数の精度が悪くなる場合がある（例:棄却法で生成した指数乱数），そこで以下に棄却法による乱数生成の方法を紹介する．用いる関数を先に説明する．用いる関数・定数を先に説明する．

-- f(x)   ： 生成したい乱数が従う分布の密度関数
-- g(x)   ： 乱数生成が容易な分布の密度関数
-- c(>=1) ： f(x) <= c*g(x) が成り立つ定数(小さい方が良い)
-- h(x)   =  f(x)/{c*g(x)}

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

- u <- runif(1)
- v <- g(x)に従う乱数
- u <= h(v) ならば v を乱数として採択する

** f(x) = 1/2*sin(x)  (0 < x < pi) という分布に従う乱数を生成する(2005.7.23) [#a808384d]

--g(x) = 1/pi   (0 < x < pi)
--c    = pi/2

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

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)
 }
 > yy <- myrand(1000)  # 乱数を1000個生成
 > plot(density(yy))   # し、このデータについて密度推定

#ref(rand01.png, center)

* アドオンパッケージ gld (generalized (Tukey) lambda distribution) 2004.09.15 [#mc5e9dd8]
アドオンパッケージ gld (Tukey の一般化ラムダ分布) に関する関数がある。この分布は少数のパラメータで様々な形状の確率分布を表現できる。以下はそのインデックスである。
 rgl                 random numbers from the generalised lambda distribution
 qgl                 quantiles of the generalised lambda distribution
 pgl                 probabilities of the generalised lambda distribution
 dgl                 densities of the generalised lambda distribution
 qdgld               quantile densities of the generalised lambda distribution
 gl.check.lambda     Function to check the validity of parameters
                       of the generalized lambda distribution
 plotgl              Plots of density and distribution function for
                       the generalised lambda distribution
 plotgld             Plot probability density functions of the generalised 
 	              lambda distribution
 plotglc             Plot distribution functions of the generalised lambda  distribution
 qqgl	            Quantile-Quantile plot against the generalised lambda  distribution
 starship            Carry out the "starship" estimation method for
                       the generalised lambda distribution
 starship.adaptivegrid Carry out the "starship" estimation method for
                         the generalised lambda distribution using a   grid-based search
 starship.obj        Objective function that is minimised in starship  estimation method

*3変量正規分布の乱数 [#eac41cf8]
は以下のように生成することが出来る．以下で用いている行列の平方根についてはこちら．

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 分布の乱数 [#w1f3047a]
は以下のように生成することが出来る．まず，アルゴリズムを書く代わりに命題みたいなものを示す．

R               ： 自由度 m のカイ二乗分布に従う確率変数 
 Z               ： p 変量正規分布 N( 0 , I p ) に従う確率ベクトル（独立な標準正規乱数がズラッと縦に並んでいる） 
 V ： 正定値対称行列（ちらばりを表す行列で分散共分散行列とは少し異なるもの）
 V = C\%*\%t(C)  ： コレスキー分解 
 X = sqrt(m/R)*Z ： p 変量楕円 t 分布 met( p , m , 0 , I p ) に従う 
 Y = μ + C %*% X ： p 変量楕円 t 分布 met( p , m , μ , V )に従う

以下に3変量楕円t乱数に従う乱数を生成する関数を定義する．

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)
 }

タイムスタンプを変更しない