Rでの分布作成と自作分布に従う乱数生成

ピークが2つ以上ある分布をRで作成してみた。関数の合成です。そしてその自作分布に従う乱数を発生させました。

まず下のような正規分布の関数fxとgxを設定します。平均と分散はそれぞれ異なる値とします。

fx <- function(x) dnorm(x, mean=1, sd = 1)
gx <- function(x) dnorm(x, mean=5, sd= 2)

次にfxとgxを重ね合わせた関数hxを定義し、プロットしてみます。

hx <- function(x) 1/2*(fx(x)+gx(x))

plot(hx, xlim=c(-10,10))

hxをfx, gxと重ねて図に描くと下のようになります。このあたりは舟尾先生のHP「R-Tips」を参考にしました。

plot(hx, xlim=c(-10,10), ylim=c(0,0.5),ylab="hx")
par(new=T) 　　　　　　　　#重ね書きの指定
plot(gx, xlim=c(-10,10), ylim=c(0,0.5),ylab="", col="red")　　　　#y軸の値を上と揃える　プロットは赤色　軸のラベルを記載しない
par(new=T)
plot(fx, xlim=c(-10,10), ylim=c(0,0.5),ylab="", col="blue")

ここからは、作成した分布hxに従う乱数を生成します。やはり舟尾先生のHP「R-Tips」を参考にしました。ほぼ丸パクリ。特に「棄却法による乱数の作り方」という項です。

目標分布、つまり乱数を発生させたい分布は上記のhxです。提案分布は0.15の一定値（-10<x<10）としましょう。つまりhx(V)<0.3ならば乱数Vを採択させます。

hx.rand <- function(x) { 　　　#xは最終的に発生する乱数の数

y <- c(); i <- 1

while (i <= x) {

U <- runif(1); V <-runif(1, -10, 10); W <- hx(V)/0.3

if (U < W){

y <- append(y, V); i <- i+1

}

}

return(y)

}

hx.data <- hx.rand(10000) 　　　　#hxに従う乱数を10000個発生

乱数生成関数が正しいかどうか、分布hxと重ね合わせてプロットして確認します。

hist(hx.data, xlim=c(-5,10),ylim=c(0,0.25),prob=T,ann=F)
par(new=T)
plot(hx, xlim=c(-5,10), ylim=c(0,0.25),ylab="",xlab="",col="red")

重なり合っているので、乱数生成の関数が正しかったことが分かりました。

上記のように分布hxの最大値をグラフから求めるのも面倒なので、自動化してみました。太字部分が変更点。

hx.rand <- function(x) { 　　　#xは最終的に発生する乱数の数

y <- c(); i <- 1

while (i <= x) {

hx.max <- max(hx(seq(-10,10,0.1)))

U <- runif(1); V <-runif(1, -10, 10); W <- hx(V)/(ceiling(hx.max))

if (U < W){

y <- append(y, V); i <- i+1

}

}

return(y)

}

hx.data <- hx.rand(10000)

ちなみに分布の積分は下のようにintegrate関数を用いる。

integrate(hx,-10,20)
1 with absolute error < 4.5e-05

（追記: 2015.02.24）

二項分布で同様のことをおこなってみました。

fx <- function(x) dbinom(x, size=1, prob = 0.5)　　#sizeは試行回数

gx <- function(x) dbinom(x, size=20, prob = 0.5)

hx <- function(x) 1/12*fx(x)+1/4*gx(x)

plot(hx,xlim=c(-5,50),type="l")
50 件以上の警告がありました (最初の 50 個の警告を見るには warnings()

warnings()
警告メッセージ:
1: In dbinom(x, size = 1, prob = 0.5) : non-integer x = -4.450000
2: In dbinom(x, size = 1, prob = 0.5) : non-integer x = -3.900000
3: In dbinom(x, size = 1, prob = 0.5) : non-integer x = -3.350000

（以下略）