乱数Tips大全
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
SIZE(25){COLOR(red){乱数 Tips 大全}}
//間瀬(2003.08.17)
#contents
~
*Rにおける(疑似)乱数発生システムの概要 [#y9d2907f]
R は確率分布に関する豊富な関数群を持ち，従来統計学の利用...
*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.9345...
[1] -2.934561
> rnorm(3, mean=2, sd=3) # 疑似乱数を 3 つ生成
[1] 5.287193 -0.211181 4.648124
> rnorm(3, m=2, s=3) # 省略形(引数名は他と一意に区別で...
> rnorm(3, 2, 2) # 省略形
*R の基本確率関数一覧 [#t313e777]
|分布名 | R での関数名 | パラメータ引数名 |
|ベータ | beta |shape1, shape2, ncp|
|2項 | binom | size, prob|
|コーシー | cauchy | location, scale|
|カイ自乗 | chisq | df, ncp|
|指数 | exp | rate |
|F | f | df1, df1, ncp|
|ガンマ | gamma | shape, scale|
|幾何 | geom | prob|
|超幾何 | hyper | m, n, k|
|対数正規 | lnorm | meanlog, sdlog|
|ロジスティック | logis | location, scale|
|多項 | multinom | n, size, prob 注1|
|負の2項 | nbinom | size, prob|
|正規 | norm | mean, sd|
|ポアソン | pois | lambda|
|符号付順位和 | signrank | n |
|t | t | df, ncp|
|ステューデント化した範囲 | tukey | nmeans, df, nranges ...
|一様 | unif | min, max|
|ワイブル | weibull | shape, scale|
|ウィルコクソン | wilcox | m, n|
注1：dmultinom, rmultinom のみ~
注2：ptukey, qtukey のみ
*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 ...
設定する。もし NULL ならば現在の手法を返す。"def...
手法に戻す。
normal.kind: 文字、もしくは NULL。もし kind が文字列な...
法を指定した方法に設定する。もし NULL ならば現在...
にすると現在の既定手法に戻す。
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 法。周期 ...
-"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', 原論文を修正した Applie...
-"Marsaglia-Multicarry"' キャリーつきの乗算 RNG を使い、G...
メイリングリスト `sci.stat.math' への投稿記事で紹介した。...
すべての検査をパスする (Marsaglia による)。乱数種は二つの...
-"Super-Duper"' Marsaglia の有名な７０年代の Super-Duper ...
二つの種はそれぞれ Tausworthe と congruence 倍長整数であ...
-"Knuth-TAOCP" Knuth (1997) による。減算を伴う遅延つきフ...
GFSR 法。つまり、使用される再帰式は X[j] = (X[j-100] - X[...
-"Knuth-TAOCP-2002" 以前のバージョンとは上位互換ではない ...
-"user-supplied" ユーザー提供の発生法を使用。詳細は `Rand...
*正規疑似乱数発生法 `normal.kind' [#w86e8ffc]
-"Inversion"' (現在の既定手法)
- "Kinderman-Ramage" (R 1.7.0 以前の既定手法、近似誤差を...
- "Buggy Kinderman-Ramage"
- "Ahrens-Dieter"
- "Box-Muller"
-"user-supplied"
`set.seed' はその単独の整数引数を用い、要求されただけの種...
値:
-`.Random.seed' は整数ベクトルであり、その最初の要素は RN...
-背後にある C コードでは `.Random.seed[-1]' は `unsigned'...
-`RNGkind' は、呼び出しの前に使用中の RNG と正規乱数発生...
-`RNGversion' は同じ情報を返す。
-`set.seed' は `NULL' を返すが、コンソールには現れない。
注意:
-最初はシードは無い。必要な時現在の時間からつくり出される...
-.Random.seed' は、少なくともシステムの発生法に対する、一...
例
runif(1); .Random.seed; runif(1); .Random.seed
## もしシードがなければ新しいものが"ランダム"につく...
rm(.Random.seed); runif(1); .Random.seed
RNGkind("Wich") # (kind に対する部分文字列マッチング)
## 以下は runif(.) が Wichmann-Hill 法に対しどのよ...
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) %%...
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組みの一様疑似乱数...
Mersenne-Twister法使用(R 1.7.0 以降の既定手法)
random1 <- function () {
old.par <- par(no.readonly = TRUE); on.exit(par(old.pa...
png("MTrand.png")
plot(c(0,0,1,1), c(0,1,0,1), main="", xlab="", ylab=""...
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.pa...
png("MMrand.png")
plot(c(0,0,1,1), c(0,1,0,1), main="", xlab="", ylab=""...
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 =...
- 引数
-- COLOR(magenta){x} ベクトル(数値、複素数、文字列、論理値)
-- COLOR(magenta){size} 選び出される個数を表す正整数
-- COLOR(magenta){replace} 論理値。復元抽出を行なうか？
-- COLOR(magenta){prob} ベクトルから抽出する際の重み確率
- 注意
-- x が長さ 1 なら 1:x から抽出
-- size の既定値は length(x)。sample(x) はランダムな並べ...
-- prob 比率を意味し、総和が 1 でなくても良いが、非負で、...
> 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 ...
[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...
[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 であるランダムな二次元分...
> 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) [#e2a1ab34]
(疑似)乱数は R 起動時に適当に初期化(システム時間を利用？)...
> 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]
上に紹介している分布に無い分布に従う乱数を生成したくなる...
-P(X ≦ x) = P(F^{-1} (U)≦x) = P(U≦F(x)) = F(x)
となることより X は分布関数 F に従う．よって，分布関数の...
-指数分布の分布関数 F(x) = 1-exp(-x) の逆関数は -log(1-x)...
-ラプラス分布（両側指数分布）の密度関数は f(x) = 1/2 exp(...
**両側指数分布の乱数を発生させる(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...
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) という分布に従う乱数を...
--g(x) = 1/pi (0 < x < pi)
--c = pi/2
とおくと f(x) ≦ c*g(x) が成り立ち，この場合は h(x) = sin(...
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 dist...
アドオンパッケージ gld (Tukey の一般化ラムダ分布) に関す...
rgl random numbers from the generalised ...
qgl quantiles of the generalised lambda ...
pgl probabilities of the generalised lam...
dgl densities of the generalised lambda ...
qdgld quantile densities of the generalise...
gl.check.lambda Function to check the validity of pa...
of the generalized lambda distribu...
plotgl Plots of density and distribution fu...
the generalised lambda distribution
plotgld Plot probability density functions o...
lambda distribution
plotglc Plot distribution functions of the g...
qqgl Quantile-Quantile plot against the gene...
starship Carry out the "starship" estimation ...
the generalised lambda distribution
starship.adaptivegrid Carry out the "starship" estimatio...
the generalised lambda distribut...
starship.obj Objective function that is minimised...
*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 ...
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...
}
return(w)
}
終了行:
SIZE(25){COLOR(red){乱数 Tips 大全}}
//間瀬(2003.08.17)
#contents
~
*Rにおける(疑似)乱数発生システムの概要 [#y9d2907f]
R は確率分布に関する豊富な関数群を持ち，従来統計学の利用...
*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.9345...
[1] -2.934561
> rnorm(3, mean=2, sd=3) # 疑似乱数を 3 つ生成
[1] 5.287193 -0.211181 4.648124
> rnorm(3, m=2, s=3) # 省略形(引数名は他と一意に区別で...
> rnorm(3, 2, 2) # 省略形
*R の基本確率関数一覧 [#t313e777]
|分布名 | R での関数名 | パラメータ引数名 |
|ベータ | beta |shape1, shape2, ncp|
|2項 | binom | size, prob|
|コーシー | cauchy | location, scale|
|カイ自乗 | chisq | df, ncp|
|指数 | exp | rate |
|F | f | df1, df1, ncp|
|ガンマ | gamma | shape, scale|
|幾何 | geom | prob|
|超幾何 | hyper | m, n, k|
|対数正規 | lnorm | meanlog, sdlog|
|ロジスティック | logis | location, scale|
|多項 | multinom | n, size, prob 注1|
|負の2項 | nbinom | size, prob|
|正規 | norm | mean, sd|
|ポアソン | pois | lambda|
|符号付順位和 | signrank | n |
|t | t | df, ncp|
|ステューデント化した範囲 | tukey | nmeans, df, nranges ...
|一様 | unif | min, max|
|ワイブル | weibull | shape, scale|
|ウィルコクソン | wilcox | m, n|
注1：dmultinom, rmultinom のみ~
注2：ptukey, qtukey のみ
*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 ...
設定する。もし NULL ならば現在の手法を返す。"def...
手法に戻す。
normal.kind: 文字、もしくは NULL。もし kind が文字列な...
法を指定した方法に設定する。もし NULL ならば現在...
にすると現在の既定手法に戻す。
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 法。周期 ...
-"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', 原論文を修正した Applie...
-"Marsaglia-Multicarry"' キャリーつきの乗算 RNG を使い、G...
メイリングリスト `sci.stat.math' への投稿記事で紹介した。...
すべての検査をパスする (Marsaglia による)。乱数種は二つの...
-"Super-Duper"' Marsaglia の有名な７０年代の Super-Duper ...
二つの種はそれぞれ Tausworthe と congruence 倍長整数であ...
-"Knuth-TAOCP" Knuth (1997) による。減算を伴う遅延つきフ...
GFSR 法。つまり、使用される再帰式は X[j] = (X[j-100] - X[...
-"Knuth-TAOCP-2002" 以前のバージョンとは上位互換ではない ...
-"user-supplied" ユーザー提供の発生法を使用。詳細は `Rand...
*正規疑似乱数発生法 `normal.kind' [#w86e8ffc]
-"Inversion"' (現在の既定手法)
- "Kinderman-Ramage" (R 1.7.0 以前の既定手法、近似誤差を...
- "Buggy Kinderman-Ramage"
- "Ahrens-Dieter"
- "Box-Muller"
-"user-supplied"
`set.seed' はその単独の整数引数を用い、要求されただけの種...
値:
-`.Random.seed' は整数ベクトルであり、その最初の要素は RN...
-背後にある C コードでは `.Random.seed[-1]' は `unsigned'...
-`RNGkind' は、呼び出しの前に使用中の RNG と正規乱数発生...
-`RNGversion' は同じ情報を返す。
-`set.seed' は `NULL' を返すが、コンソールには現れない。
注意:
-最初はシードは無い。必要な時現在の時間からつくり出される...
-.Random.seed' は、少なくともシステムの発生法に対する、一...
例
runif(1); .Random.seed; runif(1); .Random.seed
## もしシードがなければ新しいものが"ランダム"につく...
rm(.Random.seed); runif(1); .Random.seed
RNGkind("Wich") # (kind に対する部分文字列マッチング)
## 以下は runif(.) が Wichmann-Hill 法に対しどのよ...
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) %%...
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組みの一様疑似乱数...
Mersenne-Twister法使用(R 1.7.0 以降の既定手法)
random1 <- function () {
old.par <- par(no.readonly = TRUE); on.exit(par(old.pa...
png("MTrand.png")
plot(c(0,0,1,1), c(0,1,0,1), main="", xlab="", ylab=""...
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.pa...
png("MMrand.png")
plot(c(0,0,1,1), c(0,1,0,1), main="", xlab="", ylab=""...
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 =...
- 引数
-- COLOR(magenta){x} ベクトル(数値、複素数、文字列、論理値)
-- COLOR(magenta){size} 選び出される個数を表す正整数
-- COLOR(magenta){replace} 論理値。復元抽出を行なうか？
-- COLOR(magenta){prob} ベクトルから抽出する際の重み確率
- 注意
-- x が長さ 1 なら 1:x から抽出
-- size の既定値は length(x)。sample(x) はランダムな並べ...
-- prob 比率を意味し、総和が 1 でなくても良いが、非負で、...
> 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 ...
[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...
[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 であるランダムな二次元分...
> 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) [#e2a1ab34]
(疑似)乱数は R 起動時に適当に初期化(システム時間を利用？)...
> 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]
上に紹介している分布に無い分布に従う乱数を生成したくなる...
-P(X ≦ x) = P(F^{-1} (U)≦x) = P(U≦F(x)) = F(x)
となることより X は分布関数 F に従う．よって，分布関数の...
-指数分布の分布関数 F(x) = 1-exp(-x) の逆関数は -log(1-x)...
-ラプラス分布（両側指数分布）の密度関数は f(x) = 1/2 exp(...
**両側指数分布の乱数を発生させる(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...
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) という分布に従う乱数を...
--g(x) = 1/pi (0 < x < pi)
--c = pi/2
とおくと f(x) ≦ c*g(x) が成り立ち，この場合は h(x) = sin(...
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 dist...
アドオンパッケージ gld (Tukey の一般化ラムダ分布) に関す...
rgl random numbers from the generalised ...
qgl quantiles of the generalised lambda ...
pgl probabilities of the generalised lam...
dgl densities of the generalised lambda ...
qdgld quantile densities of the generalise...
gl.check.lambda Function to check the validity of pa...
of the generalized lambda distribu...
plotgl Plots of density and distribution fu...
the generalised lambda distribution
plotgld Plot probability density functions o...
lambda distribution
plotglc Plot distribution functions of the g...
qqgl Quantile-Quantile plot against the gene...
starship Carry out the "starship" estimation ...
the generalised lambda distribution
starship.adaptivegrid Carry out the "starship" estimatio...
the generalised lambda distribut...
starship.obj Objective function that is minimised...
*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 ...
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...
}
return(w)
}
ページ名:
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
