时间戳

首页 > 解决方案 > 使用接受-拒绝方法模拟随机变量

r - 使用接受-拒绝方法模拟随机变量

问题描述

我有以下算法

步骤 1. 用 qj=P(Y=j) 模拟 Y 的值

步骤 2. 生成统一变量

步骤 3. 如果 U<=Pj/(c*qj) 则 X=j 并停止。否则返回步骤 1。

而具体的例子：

X=1,2,3,4 与 P1=.2,P2=.15,P3=.25,P4=.4

为 X 生成值

  1. 让 Y~UD(1,4)

  2. c=.4/.25

这是我在 R 中实现此算法的方法：

f<-function(){
n<-4

probY<-c(.25,.25,.25,.25)
probX<-c(2,.15,.25,.4)

X<-rep(0,4)

U<-runif(n)

c<-.4/.25

for(j in 1:4)
{
if(U[j]<= probX[j]/(c*probY[j]))
   X<-j

}


return(X)

}

输出是4我认为不正确的。我不确定我是否应该写Y<-runif(n,1,4)这行probY<-c(.25,.25,.25,.25)。还有这一行“否则返回第 1 步”。循环中缺少，但总是相同的 .25

有人可以帮忙吗？

标签： ralgorithmsimulation

解决方案

我认为这里的问题与算法的工作原理有点混淆。

为了从您的分布中生成单个值 (X = 1,2,3,4 其中 P(X = 1) = .2, P(X = 2) = .15, P(X = 3) = .25 , P(X = 4) = .4)，我们需要遵循算法的步骤。假设我们选择c = .4/.25：

1.从 Y ~ UD (1,4)生成y 。

2.从 U ~ U(0,1)生成u 。

3.检查是否u≤f(y)/cg(y)。如果是，定义x = y就完成了！如果不是，请返回步骤 1。

在您提供的代码中，您实际上从未生成y变量。这是一个应该为您工作的功能！希望我的评论能很好地解释它！

accRej <- function(){

  #The probabilities for generating a r.v. from X
  probX <- c(.2,.15,.25,.4)

  #The Value for c
  c <- .4/.25

  #x is a placeholder for our final value of x 
  #and i is a counter variable which will allow us to stop the loop when we complete step 3
  x <- numeric(1)
  i <- 1

  #Now, start the loop!
  while(i <= 1){
    #Step 1, get y
    y <- sample(1:4,1)
    #Step 2, get u
    u <- runif(1)
    #Step 3, check to see if the inequality holds
    #If it does, assign y to x and add 1 to i (making it 2) to stop the while loop
    #If not, do nothing and try again!
    if(u <= probX[y]/c*.25){
      x[i] <- y
      i <- i+1
    }
  }
  #Return our value of x
  return(x)
}

请注意，在这段代码中，在我们的算法中probX[i]它等于f(y)，并且由于 Y~UD(1,4)，始终.25= g(y)。希望这可以帮助！

此外，这里是使用此方法生成n随机变量的代码。它与上面基本相同，只是可以将 1 更改为n.

accRej <- function(n){

  #The probabilities for generating a r.v. from X
  probX <- c(.2,.15,.25,.4)

  #The Value for c
  c <- .4/.25

  #x is a placeholder for our final value of x 
  #and i is a counter variable which will allow us to stop the loop when we complete step 3
  x <- numeric(n)
  i <- 1

  #Now, start the loop!
  while(i <= n){
    #Step 1, get y
    y <- sample(1:4,1)
    #Step 2, get u
    u <- runif(1)
    #Step 3, check to see if the inequality holds
    #If it does, assign y to x and add 1 to i
    #If not, do nothing and try again!
    if(u <= probX[y]/c*.25){
      x[i] <- y
      i <- i+1
    }
  }
  #Return our value of x
  return(x)
}

推荐阅读