베타 배포에 대한 수락 거부 R 코드

Question

g (x) = 1, 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?에있는 hard-wire?)도 적절한 임의 변수를 생성 할 수있는 함수를 전달해야합니다.

3) u은 [0,1]이 아닌 [0,0.8]에서 가져와야합니다. 0.8은 값의 생성에 참여하지 않아야하며 해석 만하면됩니다.

4) c은 f(y)/g(y)의 상한이어야합니다. 2가 너무 작습니다. 왜 f의 파생물을 가지고 그것이 최대라고 생각하지 않습니까? 3.5 작품. 또한 - c은 R의 변수에 적합하지 않습니다 (함수가 c()이므로). 왜 전화하지 M?

만들기 이러한 변화 수율 :

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) 

일반적인 출력 : 0.9016

또한 result의 밀도 히스토그램을하고 이론적 인 베타 분포를 비교할 수 있습니다

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

경기가 아주 좋습니다 :