0

拒绝抽样

我正在使用截断正态分布的拒绝抽样，请参阅下面的 r 代码。如何使采样停止在特定的 n？例如 1000 个观察值。即当接受的样本数量达到n（1000）时，我想停止采样。

有什么建议么？任何帮助是极大的赞赏 ：）

#Truncated normal curve    
curve(dnorm(x, mean=2, sd=2)/(1-pnorm(1, mean=2, sd=2)),1,9)

#create a data.frame with 100000 random values between 1 and 9

sampled <- data.frame(proposal = runif(100000,1,9))
sampled$targetDensity <- dnorm(sampled$proposal, mean=2, sd=2)/(1-pnorm(1, mean=2, sd=2))

#accept proportional to the targetDensity

maxDens = max(sampled$targetDensity, na.rm = T)
sampled$accepted = ifelse(runif(100000,0,1) < sampled$targetDensity / maxDens, TRUE, FALSE)

hist(sampled$proposal[sampled$accepted], freq = F, col = "grey", breaks = 100, xlim = c(1,9), ylim = c(0,0.35),main="Random draws from skewed normal, truncated at 1")
curve(dnorm(x, mean=2, sd=2)/(1-pnorm(1, mean=2, sd=2)),1,9, add =TRUE, col = "red", xlim = c(1,9),  ylim = c(0,0.35))



X <- sampled$proposal[sampled$accepted]

采样时如何将 X 的长度设置为特定数字？

4

1 回答 1

0

睡在上面之后，如果您决定使用拒绝抽样并且只在 1,000 次过去之前这样做，我认为没有比仅使用 while 循环更好的选择了。这明显低于

sampled$accepted = ifelse(runif(100000,0,1) < sampled$targetDensity / maxDens, TRUE, FALSE)
X <- sampled$proposal[sampled$accepted][1:1000]

上述代码所用时间为0.0624001s. 下面的代码所花费的时间是0.780005s。我包括它是因为它是对您所问的特定问题的答案，但这种方法效率低下。如果有其他选择，我会使用它。

#Number of samples
N_Target <- 1000
N_Accepted <- 0

#Loop until condition is met
i = 1
sampled$accepted = FALSE
while( N_Accepted < N_Target ){

    sampled$accepted[i] = ifelse(runif(1,0,1) < sampled$targetDensity[i] / maxDens, TRUE, FALSE)
    N_Accepted = ifelse( sampled$accepted[i], N_Accepted + 1 , N_Accepted )
    i = i + 1
    if( i > nrow( sampled ) ) break

}
于 2018-01-05T10:48:38.237 回答