Loop While 条件为 TRUE

Question

我正在尝试生成 n 个总和小于 1 的随机数。

所以我不能只运行runif(3)。但是我可以将每次迭代的条件设置为到该点生成的所有值的总和。

这个想法是开始一个空向量v，并设置一个循环，以便对于每次迭代i，生成一个runif()，但在它被接受为v的元素之前，即v[i] <- runif()，测试sum(v) < 1，而FALSE最后一个条目v[i]最终被接受，但是如果TRUE，即总和大于1，v[i]被扔出向量，重复迭代i。

我远未实现这个想法，但我想按照类似于以下内容的方式解决它。这与其说是一个实际问题，不如说是一个练习，以了解一般循环的语法：

n <- 4
v <- 0

for (i in 1:n){
    rdom <- runif(1)
    if((sum(v) + rdom) < 1) v[i] <- rdom
    }
    # keep trying before moving on to iteration i + 1???? i <- stays i?????
} 

我研究了while（实际上我在标题中加入了while 函数）；但是，我需要向量具有n 元素，所以如果我尝试一些基本上告诉R 添加随机统一实现作为向量v while sum(v) < 1 的元素的方法，我会卡住，因为我可以结束v 中的元素少于 n。

Answer 1

这是一个可能的解决方案。 它不使用while，而是更通用的repeat。我将它编辑为使用while并保存了几行。

set.seed(0)
n <- 4
v <- numeric(n)
i <- 0
while (i < n) {
  ith <- runif(1)
  if (sum(c(v, ith)) < 1) {
    i <- i+1
    v[i] <- ith
  }
}
v
# [1] 0.89669720 0.06178627 0.01339033 0.02333120

使用repeat 块，无论如何您都必须检查条件，但是，消除日益严重的问题，它看起来非常相似：

set.seed(0)
n <- 4
v <- numeric(n)
i <- 0
repeat {
  ith <- runif(1)
  if (sum(c(v, ith)) < 1) {
    i <- i+1
    v[i] <- ith
  }
  if (i == 4) break
} 

Answer 2

如果您真的想保持与您发布的完全相同的过程（也就是从标准均匀分布中一次一个迭代地采样 n 值，拒绝任何导致总和超过 1 的样本），那么以下代码在数学上是等效的、更短且更高效：

samp <- function(n) {
  v <- rep(0, n)
  for (i in 1:n) {
    v[i] <- runif(1, 0, 1-sum(v))
  }
  v
}

基本上，这段代码使用的数学事实是，如果向量的和当前为sum(v)，那么从标准均匀分布中采样直到得到不大于1-sum(v) 的值与在均匀分布中采样完全相同分布从 0 到 1-sum(v)。使用后一种方法的优点是效率更高——我们不需要一直拒绝样本并再次尝试，而是只需为每个元素采样一次。

要了解运行时差异，请考虑使用 n=10 抽样 100 个观察值，与您帖子中代码的工作实现进行比较（复制自我对此问题的其他答案）：

OP <- function(n) {
  v <- rep(0, n)
  for (i in 1:n){
    rdom <- runif(1)
    while (sum(v) + rdom > 1) rdom <- runif(1)
    v[i] <- rdom
  }
  v
}
set.seed(144)
system.time(samples.OP <- replicate(100, OP(10)))
#    user  system elapsed 
# 261.937   1.641 265.805 
system.time(samples.josliber <- replicate(100, samp(10)))
#    user  system elapsed 
#   0.004   0.001   0.004

在这种情况下，新方法的速度快了近 100,000 倍。

Answer 3

听起来您正试图从具有以下约束的n 变量空间中统一采样：

x_1 + x_2 + ... + x_n <= 1
x_1 >= 0
x_2 >= 0
...
x_n >= 0

"hit and run" algorithm 是一种数学机制，可以让您做到这一点。在二维空间中，算法会从下面的三角形中均匀采样，阴影区域中的每个位置被选中的可能性相同：

该算法在 R 中通过 hitandrun 包提供，它要求您通过约束矩阵、方向向量和右侧向量指定定义空间的线性不等式：

library(hitandrun)
n <- 3
constr <- list(constr = rbind(rep(1, n), -diag(n)),
               dir = c(rep("<=", n+1)),
               rhs = c(1, rep(0, n)))
set.seed(144)
samples <- hitandrun(constr, n.samples=1000)
head(samples, 10)
#             [,1]       [,2]       [,3]
#  [1,] 0.28914690 0.01620488 0.42663224
#  [2,] 0.65489979 0.28455231 0.00199671
#  [3,] 0.23215115 0.00661661 0.63597912
#  [4,] 0.29644234 0.06398131 0.60707269
#  [5,] 0.58335047 0.13891392 0.06151205
#  [6,] 0.09442808 0.30287832 0.55118290
#  [7,] 0.51462261 0.44094683 0.02641638
#  [8,] 0.38847794 0.15501252 0.31572793
#  [9,] 0.52155055 0.09921046 0.13304728
# [10,] 0.70503030 0.03770875 0.14299089

稍微分解一下这段代码，我们生成了以下约束矩阵：

constr
# $constr
#      [,1] [,2] [,3]
# [1,]    1    1    1
# [2,]   -1    0    0
# [3,]    0   -1    0
# [4,]    0    0   -1
# 
# $dir
# [1] "<=" "<=" "<=" "<="
# 
# $rhs
# [1] 1 0 0 0

阅读constr$constr 的第一行，我们有 1, 1, 1 表示“1*x1 + 1*x2 + 1*x3”。 constr$dir的第一个元素是<=，constr$rhs的第一个元素是1；把它放在一起我们有x1 + x2 + x3 <= 1。从constr$constr 的第二行我们读取 -1, 0, 0 表示“-1*x1 + 0*x2 + 0*x3”。 constr$dir的第二个元素是<=，constr$rhs的第二个元素是0；把它放在一起我们有-x1 <= 0，这和x1 >= 0是一样的。其余行中遵循类似的非负约束。

请注意，hit and run 算法的优点是每个变量的分布完全相同：

hist(samples[,1])

hist(samples[,2])

hist(samples[,3])

同时，您的程序中的样本分布将非常不均匀，并且随着n 的增加，此问题将变得越来越严重。

OP <- function(n) {
  v <- rep(0, n)
  for (i in 1:n){
    rdom <- runif(1)
    while (sum(v) + rdom > 1) rdom <- runif(1)
    v[i] <- rdom
  }
  v
}
samples.OP <- t(replicate(1000, OP(3)))

hist(samples.OP[,1])

hist(samples.OP[,2])

hist(samples.OP[,3])

一个额外的优势是 hit-and-run 算法看起来更快 - 我在我的计算机上使用 hit-and-run 在 0.006 秒内生成了这 1000 个重复，而使用 OP 中的修改代码则需要 0.3 秒。

Answer 4

我会这样做，没有任何循环，if 或 while：

set.seed(123)
x <- runif(1) # start with the sum that you want to obtain
n <- 4 # number of generated random numbers, can be chosen arbitrarily 
y <- sort(runif(n-1,0,x)) # choose n-1 random points to cut the range [0:x]
z <- c(y[1],diff(y),x-y[n-1]) # result: determine the length of the segments
#> z
#[1] 0.11761257 0.10908627 0.02723712 0.03364156
#> sum(z)
#[1]  0.2875775
#> all.equal(sum(z),x)
#[1] TRUE

这里的好处是您可以准确地确定要获得的总和以及要为此生成多少个数字n。例如，如果在第二行设置x <- 1，则存储在向量z 中的n 随机数将加起来为1。