在 R 中生成满足约束的随机数

Question

我需要有关根据约束生成随机数的代码方面的帮助。具体来说，我试图分别从正态分布和伽马分布中模拟随机数 ALFA 和 BETA，使得 ALFA - BETA < 1.

这是我写的，但根本不起作用。

set.seed(42) 
n <- 0 
repeat {
  n <- n + 1
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) 
  alfa[n] <- a
  beta[n] <- b
  l = length(alfa)
  if (l == 10000) break
}

Answer 1

由于矢量化，生成数字 "all at once" 比循环生成数字更快：

set.seed(42)
N = 1e5
a = rnorm(N, 10, 2)
b = rgamma(N, 8, 1)
d = a - b
alfa = a[d < 1]
beta = b[d < 1]
length(alfa)
# [1] 36436

这产生了 100,000 名候选人，其中 36,436 名符合您的标准。如果您想生成 n 个样本，请尝试设置 N = 4 * n，您可能会生成足够多的样本，保留第一个 n.

您的循环有两个问题：(a) 您需要花括号将 if 语句后的多行括起来。 (b) 您正在使用 n 作为尝试计数器，但它应该是成功计数器。正如所写，只有在第 10000 次尝试成功时，您的循环才会停止。将 n <- n + 1 移到 if 语句中以修复：

set.seed(42) 
n <- 0 
alfa = numeric(0)
beta = numeric(0)
repeat {
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) {
    n <- n + 1
    alfa[n] <- a
    beta[n] <- b
    l = length(alfa)
    if (l == 500) break
  }
}

但是第一种方式更好... 由于"growing" alfa 和 beta 在循环中，一次生成一个数字，这种方法生成时间较长500 个数字比上面的代码生成 30,000 个数字要多。

Answer 2

正如@Gregor Thomas 所评论的，您的尝试失败是由于缺少花括号来括起 if 语句。如果您想跳过 {} 以获得 if 控制，也许您可以尝试下面的代码

set.seed(42) 
r <- list()
repeat {
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) r[[length(r)+1]] <- cbind(alfa = a, beta = b)
  if (length(r) == 100000) break
}
r <- do.call(rbind,r)

这样

> head(r)
          alfa      beta
[1,]  9.787751 12.210648
[2,]  9.810682 14.046190
[3,]  9.874572 11.499204
[4,]  6.473674  8.812951
[5,]  8.720010  8.799160
[6,] 11.409675 10.602608

在 R 中生成满足约束的随机数

Generate random numbers in R satisfying constraints

random

r

constraints