【发布时间】:2020-06-15 07:39:48
【问题描述】:
我想应用拒绝采样方法来模拟来自单位圆盘D = { (X_1 , X_2) \in R^2: \sqrt{x^2_1 + x^2_2} ≤ 1} 的均匀分布的随机向量Y=(Y_1, Y_2),这样X = (X_1 , X_ 2) 是正方形S = [−1, 1]^2 中均匀分布的随机向量和联合密度f(y_1,y_2) = \frac{1}{\pi} 1_{D(y_1,y_2)}.
在拒绝方法中,如果f(x) \leq C * g(x),我们一般接受样品。我正在使用以下代码:
x=runif(100,-1,1)
y=runif(100,-1,1)
d=data.frame(x=x,y=y)
disc_sample=d[(d$x^2+d$y^2)<1,]
plot(disc_sample)
我有两个问题:
{使用上面的代码,从逻辑上讲,d 的大小应该大于disc_sample 的大小,但是当我调用它们时,我看到它们每个都有 100 个元素。这怎么可能。为什么尺寸相同。}这部分已解决,感谢下面的评论。
现在的问题
另外,我怎样才能重新编写我的代码,以提供让 100 个样本符合条件所需的样本总数。即在我获得所需的 100 个样本之前给我拒绝的样本数量?
感谢r2evans 的回答,但我希望写一些更简单的东西,一个 while 循环将所有可能的样本存储在矩阵或数据帧而不是列表中,然后从该数据帧调用只是样本跟随条件。我修改了答案中的代码,没有使用列表,也没有使用 sapply 函数,但它没有给出所需的结果,它只产生一行。
i=0
samps <- data.frame()
goods <- data.frame()
nr <- 0L
sampsize <- 100L
needs <- 100L
while (i < needs) {
samps <- data.frame(x = runif(1, -1, 1), y = runif(1, -1, 1))
goods <- samps[(samps$x^2+samps$y^2)<1, ]
i = i+1
}
我也想过这个:
i=0
j=0
samps <- matrix()
goods <- matrix()
needs <- 100
while (j < needs) {
samps[i,1] <- runif(1, -1, 1)
samps[i,2] <- runif(1, -1, 1)
if (( (samps[i,1])**2+(samps[i,2])**2)<1){
goods[j,1] <- samps[i,1]
goods[j,2] <- samps[i,2]
}
else{
i = i+1
}
}
但它不起作用。
如果能帮助我修改代码,我将不胜感激。
【问题讨论】:
-
您是否在其中看到了 100 个元素,或者您在其中看到了名为
100的行?当我运行这个(前面是set.seed(42))时,我看到nrow(disc_sample)是79。 -
@r2evans 是的,你是对的。您能否帮我获得 100 个样本而不是 79 个。即除了给出获得这个数字所需的样本总数之外,我还能做些什么来让代码继续运行,直到我得到 100 个样本遵循条件?
-
在询问包含随机性的问题时,请务必包含您的随机查看。我无法重现您的 79,但我可以提供一些您应该能够在您的系统上重现的内容(答案如下)。
-
“但它不起作用” 没有帮助。什么不起作用?错误?警告?数字错误?你是否更新过
j或needs,或者你真的想永远运行你的循环?
标签: r simulation probability sampling