a record of inner life

やったことや考えたこと・本・論文・音楽の備忘録。 特に環境科学・生態毒性に関して。

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

f:id:Kyoshiro1225:20150223192546j:plain

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

f:id:Kyoshiro1225:20150223195043j:plain

 

 

ここからは、作成した分布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")

f:id:Kyoshiro1225:20150223212342j:plain

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

 

上記のように分布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

(以下略)

 

合成した関数hxのプロットをすると上記のような警告が出て、しかも出てきた図も妙なものになりました。

離散分布なので、xの値をあらかじめ整数に指定しておかないといけないようです。修正版は下。

x<--5:50
plot(x,hx(x),xlim=c(-5,50),type="l")

f:id:Kyoshiro1225:20150224135508j:plain