接受-拒绝 beta 分发 R 代码答案

【问题标题】：Acceptance-rejection for beta distribution R code接受-拒绝 beta 分发 R 代码
【发布时间】：2018-04-26 09:48:03
【问题描述】：

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

【问题讨论】：

f(x) = 10x^2(1-x)^5 是一个怎样的算法？请更详细地解释您的问题，并展示您解决问题的尝试。
这是一个函数，而不是算法。目前尚不完全清楚您打算使用此功能做什么。
您可以从在[0,1] 上找到f(x) 的上限开始

标签： r statistical-sampling

【解决方案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")

比赛还不错：

【讨论】：