【问题标题】:Comparison of two vectors resulted after simulation模拟后两个向量的比较
【发布时间】: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,但我可以提供一些您应该能够在您的系统上重现的内容(答案如下)。
  • “但它不起作用” 没有帮助。什么不起作用?错误?警告?数字错误?你是否更新过jneeds,或者你真的想永远运行你的循环?

标签: r simulation probability sampling


【解决方案1】:

至于您的第二个问题...您无法重新编写代码以准确知道获得(至少)100 个结果组合需要多少。您可以使用while 循环并连接结果,直到您至少有 100 个这样的行,然后截断超过 100 个的行。因为分段(按比例)使用熵是“昂贵的”,您可能更喜欢总是高估行您需要并一次获取所有内容。

(根据作业限制进行了编辑以减少“复杂性”。)

set.seed(42)
samps <- vector(mode = "list")
goods <- vector(mode = "list")
nr <- 0L
iter <- 0L
sampsize <- 100L
needs <- 100L
while (nr < needs && iter < 50) {
  iter <- iter + 1L
  samps[[iter]] <- data.frame(x = runif(sampsize, -1, 1), y = runif(sampsize, -1, 1))
  rows <- (samps[[iter]]$x^2 + samps[[iter]]$y^2) < 1
  goods[[iter]] <- samps[[iter]][rows, ]
  nr <- nr + sum(rows)
}
iter                # number of times we looped
# [1] 2
out <- head(do.call(rbind, goods), n = 100)
NROW(out)
# [1] 100
head(out) ; tail(out)
#            x          y
# 1  0.8296121  0.2524907
# 3 -0.4277209 -0.5668654
# 4  0.6608953 -0.2221099
# 5  0.2834910  0.8849114
# 6  0.0381919  0.9252160
# 7  0.4731766  0.4797106
#               x          y
# 221 -0.65673577 -0.2124462
# 231  0.08606199 -0.7161822
# 251 -0.37263236  0.1296444
# 271 -0.38589120 -0.2831997
# 28  -0.62909284  0.6840144
# 301 -0.50865171  0.5014720

【讨论】:

  • @r2evans 非常感谢您的回答。我尝试编写更简单的东西但没有成功,只是一个while循环将所有可能的样本存储在一个矩阵而不是一个列表中,然后从该矩阵调用只是样本遵循条件。你有什么建议我可以在没有列表和 sapply 函数的情况下简化你的代码吗?
  • 类似于上述问题编辑中的代码。
  • 好的,所以我颠倒了逻辑(&gt; 而不是&lt;)。尽管如此,结果是一个单帧(不是列表),正好有 100 行满足条件。这不是你需要的吗?
  • 谢天谢地,您的回答给出了所需的结果,但问题是我仅限于几个简单的函数,即我需要在不使用函数 listwith、@ 的情况下编写代码987654327@。所以我厌倦了像上面的编辑那样简化你的代码,但现在它只给出一行
  • 是的,那些“限制”不在您原来的问题中。 sapply 的使用只是为了演示某些东西,它与功能无关(注意 doingshowing 某物之间的区别)。 with 是减少代码的便利,您可以轻松地将with(...) 替换为(samps[[iter]]$x^2 + samps[[iter]]y^2) &gt; 1。不能使用listvector(mode="list") 呢? ( 你确定你可以使用 R 吗?你可以使用 c 吗? :-) (我厌倦了教导不良做法的家庭作业。)
猜你喜欢
  • 1970-01-01
  • 2010-12-25
  • 1970-01-01
  • 1970-01-01
  • 2016-06-29
  • 2016-02-26
  • 1970-01-01
  • 2014-02-24
  • 2010-11-26
相关资源
最近更新 更多