1

我正在使用接受-拒绝方法进行 beta 分布，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))
4

1 回答 1

1

你很接近，但有几个问题：

naccepts1) 与vs的错字。n.accepts

2）如果你不打算硬连线，g那么你不应该硬连线runif()作为生成随机变量的函数，这些随机变量根据g. 该函数rejection（为什么要硬连线beta？）也应该传递一个能够生成适当随机变量的函数。

3)u应取自[0,1]not [0,0.8]。0.8不应该参与价值观的产生，只是他们的解释。

4)c需要是 的上限f(y)/g(y)。2太小了。为什么不取 的导数f来找到它的最大值？3.5 作品。另外 --c不是 R 中变量的好名字（因为 function 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")

比赛还不错：

在此处输入图像描述

于 2017-11-14T00:45:47.953 回答