Beta 分发 R 代码的接受-拒绝

Question

我对 g(x) = 1, 0 ≤ x ≤ 1 的 beta 分布使用接受-拒绝法。功能是： 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) 如果您不打算在 g 中硬连线，那么您不应该在 runif() 中硬连线作为生成随机变量的函数，这些随机变量根据 g。函数 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？

进行这些更改会产生：

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 的密度直方图并将其与理论 beta 分布进行比较：

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

比赛相当不错：

Beta 分发 R 代码的接受-拒绝

Acceptance-rejection for beta distribution R code

r

statistical-sampling