乱数 Tips 大全


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

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

R の確率分布関数の使用例

詳しい説明は 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 の基本確率関数一覧

分布名R での関数名パラメータ引数名
ベータbetashape1, shape2, ncp
2項binomsize, prob
コーシーcauchylocation, scale
カイ自乗chisqdf, ncp
指数exprate
Ffdf1, df1, ncp
ガンマgammashape, scale
幾何geomprob
超幾何hyperm, n, k
対数正規lnormmeanlog, sdlog
ロジスティックlogislocation, scale
多項multinomn, size, prob 注1
負の2項nbinomsize, prob
正規normmean, sd
ポアソンpoislambda
符号付順位和signrankn
ttdf, ncp
ステューデント化した範囲tukeynmeans, df, nranges 注2
一様unifmin, max
ワイブルweibullshape, scale
ウィルコクソンwilcoxm, n

注1:dmultinom, rmultinom のみ
注2:ptukey, qtukey のみ

Rで使える疑似乱数発生法(R 1.7.0)以降

用法:
.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 に依存する。詳細は各手法を参照

正規疑似乱数発生法 `normal.kind'

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

値:

注意:

    例
    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])

疑似乱数発生法による違いを示す例

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()
}
 ,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()
}
 ,left

ランダム抽出と並べ変え

> 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)  # 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 = 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 列の二次元分割表を対象にする検定で,計算量が大きい場合にモンテカルロ法により近似的な解を求めるときに使える。

乱数の再現 set.seed() 関数 (2003.12.27)

(疑似)乱数は 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)

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

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

両側指数分布の乱数を発生させる(2005.7.23)

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(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)

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

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

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

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

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

とおくと 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(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)

アドオンパッケージ gld (generalized (Tukey) lambda distribution) 2004.09.15

アドオンパッケージ 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変量正規分布の乱数

は以下のように生成することが出来る.以下で用いている行列の平方根についてはこちら.

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 分布の乱数

は以下のように生成することが出来る.まず,アルゴリズムを書く代わりに命題みたいなものを示す.

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

添付ファイル: filerand01.png 2092件 [詳細] fileMTrand.png 4437件 [詳細] fileMMrand.png 4668件 [詳細] filerand01.gif 1937件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-03-25 (土) 11:19:16