ベータ版の受け入れ拒否Rコード

Question

g（x）= 1,0,0≦x≦1のベータ版の受け入れ拒否方法を使用しています。 関数は次のとおりです。 f（x）= 100x^3（1-x）^ 2である。

この密度関数からデータを生成するアルゴリズムを作成します。

k = 1000の繰り返し（n = 1000）でP（0≦X≦0.8）を推定するにはどうすればよいですか？ Rでこれをどのように解決できますか？

私が既に持っている：

beta.rejection <- function(f, c, g, n) { 
naccepts <- 0 
result.sample <- rep(NA, n) 

    while (naccepts < n) { 
    y <- runif(1, 0, 0.8) 
    u <- runif(1, 0, 0.8) 

    if (u <= f(y)/(c*g(y))) { 
     naccepts <- naccepts + 1 
     result.sample[n.accepts] = y 
    } 
    } 

    result.sample 
} 

f <- function(x) 100*(x^3)*(1-x)^2 
g <- function(x) x/x 
c <- 2 

result <- beta.rejection(f, c, g, 10000) 

for (i in 1:1000){ # k = 1000 reps 
    result[i] <- result[i]/n 
} 

print(mean(result)) 

Answer 1

あなたが接近しているが、いくつかの問題：あなたはハードに行くされていない場合は

1）naccepts対n.accepts

2）とタイプミス-wireをgに設定した場合、gに従って配布されるランダム変数を生成する関数としてrunif()をハードワイヤしてはいけません。関数rejection（betaのハードワイヤはなぜですか？）にも、適切なランダム変数を生成できる関数を渡す必要があります。

3）uは、[0,1]からではなく、[0,0.8]から引き出されます。 0.8は値の生成には関与せず、それらの解釈のみが必要です。

4）cは、f(y)/g(y)の上限である必要があります。 2が小さすぎます。それが最大だと見つけるためにfの派生物を取ってみませんか？ 3.5作品。また、 - cは、Rの変数の名前にはよくない（関数c()のため）。なぜそれをMと呼びませんか？

これらの変更の利回りを作る：あなたはまた、resultの濃度ヒストグラムを作成し、理論的なベータ分布と比較することができ0.9016

：

> hist(result,freq = FALSE) 
> points(seq(0,1,0.01),dbeta(seq(0,1,0.01),4,3),type = "l") 

rejection <- function(f, M, g, rg,n) { 
    naccepts <- 0 
    result.sample <- rep(NA, n) 

    while (naccepts < n) { 
    y <- rg(1) 
    u <- runif(1) 

    if (u <= f(y)/(M*g(y))) { 
     naccepts <- naccepts + 1 
     result.sample[naccepts] = y 
    } 
    } 

    result.sample 
} 

f <- function(x) 100*(x^3)*(1-x)^2 
g <- function(x) 1 
rg <- runif 
M <- 3.5 #upper bound on f (via basic calculus) 

result <- rejection(f, M, g,rg, 10000) 

print(length(result[result <= 0.8])/10000) 

典型的な出力を一致はかなり良いです：