【问题标题】:Generate n samples, Rejection sampling in R生成n个样本,R中的拒绝采样
【发布时间】:2018-06-14 06:29:01
【问题描述】:

拒绝抽样

我正在使用截断正态分布的拒绝抽样,请参阅下面的 r 代码。如何使采样停止在特定的 n?例如 1000 个观察值。 IE。我想在接受的样本数达到 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 的长度设置为特定的数字?

【问题讨论】:

  • 您是要执行所有计算,然后选择前 1,000 个通过的,还是只计算直到 1,000 个通过?前者很简单(sampled$proposal[sampled$accepted][1:1000]),后者需要额外的步骤,因为您必须在每次采样后检查您是否已达到 1,000 次通过。除非数据库非常大,否则我认为这不会更有效。
  • 我想做后者。你知道它可以如何执行吗? @AodhanOL
  • 我会写一些方法来做到这一点。我认为它们不是很好的 R 代码。
  • 是否有必要使用拒绝抽样?有一些相当简单的方法可以做到这一点而不会被拒绝。

标签: r statistics normal-distribution resampling statistical-sampling


【解决方案1】:

在使用它之后,如果您决定使用拒绝采样并且只在 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

}

【讨论】:

  • 我明白你所说的效率是什么意思。但是,在我的情况下,后一个选项似乎没有返回一个包含 1000 个观察值的向量。当我写 length(sampled$accepted) 我得到 100000。length(sampled$proposal[sampled$accepted]) 给我 755。
  • 哪个向量有 1000 个观测值?
  • 这表明只有 755 人通过了测试(这似乎很不寻常)。 sum(sampled$accepted) 的结果是什么?对我来说,我创建的数据是 43364。这也是我运行 lenght(sampled$proposal[sampled$accepted]) 时得到的结果。
  • 好吧,如果我将 N_Target 更改为 500,我只会得到一个包含 422 个观测值的向量。所以我不确定循环返回 1000 个观察值的向量,而是在测试 1000 个时停止,其中大约 750 个被接受。我希望它在获得 1000 个接受的观察时停止,而不是在测试 1000 个时停止。 @AodhanOL
  • 这是奇怪的行为。我也明白了——让我看看我能不能找出原因。
猜你喜欢
  • 2018-12-14
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2022-01-13
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2020-01-08
相关资源
最近更新 更多