生成n个样本,R中的拒绝采样
Generate n samples, Rejection sampling in R
拒绝抽样
我正在处理具有截断正态分布的拒绝抽样,请参阅下面的 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 次之前使用它,我认为没有比使用 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
}
拒绝抽样
我正在处理具有截断正态分布的拒绝抽样,请参阅下面的 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 次之前使用它,我认为没有比使用 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
}