logo标签列表

对乘以设计矩阵后仍然有序的向量进行采样

ralgorithmmatrixrandom
4

我想以一种方式对长度为C的向量a、b、c的元素进行抽样，使得每个向量a、b、c的元素以及D %*% A1的元素都被排序。这里，A1是通过行绑定向量a、b和c得到的矩阵。D是一个设计矩阵（见代码）。

我不能使用暴力方法来寻找解决方案，因为实际问题涉及超过三个向量（a、b、c、d等）和潜在的更高的C值。因此，暴力方法需要太长时间才能找到解决方案。

最小示例：

# Size of the vectors
C <- 3

# Sample sorted vectors
a <- sort(runif(C, 0, 1))
b <- sort(runif(C, 0, 1))
c <- sort(runif(C, 0, 1))

# Row-bind the vectors a, b, c to matrix A1
A1 <- rbind(a, b, c)

# Create a design matrix D
D <- matrix(c(1, 1, 0, -1, 0, 1, 0, -1, -1), nrow = 3)

# Multiply
A2 <- D %*% A1

# Is the result as desired?
all(t(apply(A2, 1, sort)) == A2)

通用示例

创建赏金中提到的“更大案例”的代码。

# Size of the vectors
C <- 4

# Number of vectors
nvec <- 100

# Create matrix A1 with sorted rows
A1 <- vector("list", nvec)
for (i in 1:nvec) {
  A1[[i]] <- sort(runif(C, 0, 1))
}
A1 <- do.call(rbind, A1)

# Create a design matrix D
DPosition <- combn(nvec, 2)
D <- matrix(0, ncol(DPosition), nvec)
for (i in 1:nrow(D)) {
  D[i, DPosition[, i]] <- c(1, -1)
}

# Multiply
A2 <- D %*% A1

# Is the result as desired?
all(t(apply(A2, 1, sort)) == A2)

你有想法如何在R中实现这个吗？
编辑：添加了一个更一般的例子以澄清问题的规模。
你希望行排在列的前面吗？ - Maël
我只希望将A1和A2的每一行进行排序。列并不重要。 - tc90kk
1
我发布了一些代码，但是如果我们了解更广泛的任务，我们可能会给你一个更好的答案。 - David Eisenstat
2个回答
4

更新

这里有一个更新，针对f函数，使其返回一个包含给定C和nr的A1和D列表，并且A2是不可忽略的

f <- function(C, nr) {
  res <- matrix(nrow = nr, ncol = C)
  res[1, ] <- sort(runif(C, 1 - 1 / nr, 1))
  for (k in 2:nr) {
    d <- diff(res[k - 1, ])
    v <- c(runif(1, 0, d[1]), d)
    res[k, ] <- 1 - k / nr + cumsum(v - runif(1, 0, v[1]))
  }

  DPosition <- combn(nr, 2)
  D <- matrix(0, ncol(DPosition), nr)
  for (i in 1:nrow(D)) {
    D[i, DPosition[, i]] <- c(1, -1)
  }

  list(A1 = res, D = D)
}

例如

> C <- 4

> nr <- 6

> (r <- f(C, nr))
$A1
             [,1]         [,2]       [,3]       [,4]
[1,] 0.9074367087 0.9136961138 0.92941219 0.94173666
[2,] 0.6696354939 0.6726631943 0.68514756 0.69424033
[3,] 0.5017587194 0.5035297048 0.51475736 0.52259342
[4,] 0.3333633168 0.3347584666 0.34561028 0.35307051
[5,] 0.1669981771 0.1675170401 0.17749257 0.18407651
[6,] 0.0002043442 0.0006086951 0.01046971 0.01693914

$D
      [,1] [,2] [,3] [,4] [,5] [,6]
 [1,]    1   -1    0    0    0    0
 [2,]    1    0   -1    0    0    0
 [3,]    1    0    0   -1    0    0
 [4,]    1    0    0    0   -1    0
 [5,]    1    0    0    0    0   -1
 [6,]    0    1   -1    0    0    0
 [7,]    0    1    0   -1    0    0
 [8,]    0    1    0    0   -1    0
 [9,]    0    1    0    0    0   -1
[10,]    0    0    1   -1    0    0
[11,]    0    0    1    0   -1    0
[12,]    0    0    1    0    0   -1
[13,]    0    0    0    1   -1    0
[14,]    0    0    0    1    0   -1
[15,]    0    0    0    0    1   -1


> all(apply(r$A1, 1, Negate(is.unsorted)))
[1] TRUE

> (A2 <- with(r, D %*% A1))
           [,1]      [,2]      [,3]      [,4]
 [1,] 0.2378012 0.2410329 0.2442646 0.2474963
 [2,] 0.4056780 0.4101664 0.4146548 0.4191432
 [3,] 0.5740734 0.5789376 0.5838019 0.5886662
 [4,] 0.7404385 0.7461791 0.7519196 0.7576602
 [5,] 0.9072324 0.9130874 0.9189425 0.9247975
 [6,] 0.1678768 0.1691335 0.1703902 0.1716469
 [7,] 0.3362722 0.3379047 0.3395373 0.3411698
 [8,] 0.5026373 0.5051462 0.5076550 0.5101638
 [9,] 0.6694311 0.6720545 0.6746778 0.6773012
[10,] 0.1683954 0.1687712 0.1691471 0.1695229
[11,] 0.3347605 0.3360127 0.3372648 0.3385169
[12,] 0.5015544 0.5029210 0.5042876 0.5056543
[13,] 0.1663651 0.1672414 0.1681177 0.1689940
[14,] 0.3331590 0.3341498 0.3351406 0.3361314
[15,] 0.1667938 0.1669083 0.1670229 0.1671374

> all(apply(A2, 1, Negate(is.unsorted)))
[1] TRUE

之前的回答

我猜在真正“随机”抽样之前，你需要一些人工干预。这里有一个你可以尝试的例子。

f <- function(C, nr) {
  res <- matrix(nrow = nr, ncol = C)
  res[1, ] <- sort(runif(C, 0, 1))
  for (k in 2:nr) {
    d <- diff(res[k - 1, ])
    v <- c(runif(1, 0, d[1]), d)
    res[k, ] <- cumsum(v - runif(1, 0, v[1]))
  }
  res
}

这样

> (A1 <- f(3, 3))
          [,1]      [,2]      [,3]
[1,] 0.3800352 0.7774452 0.9919061
[2,] 0.1579321 0.5128165 0.6847518
[3,] 0.0979778 0.4387943 0.5966617

> all(t(apply(A1, 1, sort)) == A1)
[1] TRUE

> (A2 <- D %*% A1)
          [,1]       [,2]       [,3]
[1,] 0.2221031 0.26462868 0.30715428
[2,] 0.2820574 0.33865089 0.39524441
[3,] 0.0599543 0.07402221 0.08809012

> all(t(apply(A2, 1, sort)) == A2)
[1] TRUE
感谢您的绝佳创意！如上所述，A2的元素不一定要分布均匀随机。然而，对于C = 4和增加的nr，A2中堆叠向量之间的偏差变得极小，这对我的应用程序来说是次优的。也许这是无法避免的。尽管如此，我会等待几天，因为仍有可能有人提出不同的方法。 - tc90kk
4

这不是一个完整的答案，但对于评论来说太长了。

据我理解，目标是采样 n 个向量 v1, v2, …, vn ∈ (0, 1)d，使得所有 1 ≤ i < j ≤ n 的向量 vi − vj 具有排序条目。

将差分映射 Δ : Rd → Rd−1 定义为

(x1, x2, …, xd) ↦ (x2 − x1, x3 − x2, …, xd − xd−1).

一个向量 x ∈ Rd 如果且仅如果 Δx ≥ 0，则被排序。特别地，当且仅当 Δ(vi − vj) ≥ 0 时，向量 vi − vj 被排序。由于 Δ 是线性的，这个条件等价于 Δvi ≥ Δvj，或者，扩展全称并包含命题 vn 被排序，

Δv1 ≥ Δv2 ≥ … ≥ Δvn ≥ 0。

Thomas 的（第一个）建议是

  • 从(0,1)d中随机抽取一个向量v1，并对其进行排序（即抽取一个随机向量并对其进行排序，因为理论上退化情况的概率为零）。
  • 对于i从2到n，从(0，(Δvi-1)1) × (0，(Δvi-1)2) × … × (0，(Δvi-1)d-1)中均匀随机采样yi，然后随机采样vi∈(0,1)d，使得Δvi=yi（即从适当的范围内采样第一个元素并计算累积和）。

我已经多次说了“均匀随机”，但总体结果不是均匀随机的。与暴力拒绝采样相比，Δvn的条目倾向于小。

我怀疑，对于每个j∈{1, 2, …, d}，独立地对样本进行均匀随机(yn)j ≤ (yn−1)j ≤ … ≤ (y2)j ≤ (Δv1)j，可能会有所改善。但这仍然不正确，因为（直观上）我们缩小差异的程度越大，我们就有越多的选项来选择第一个元素。

我似乎找不到一种好的方法来随机均匀地采样v1, v2, …, vn。我找到了一种采样Δv1, Δv2, …, Δvn的方法，通过获取n(d-1)+1个指数分布的随机值，规范化额外的自由度，并计算行和列的累积和（参见从狄利克雷分布中采样）。由于您不关心采样器是否均匀，因此我们可以从产生这些增量的v1, v2, …, vn中均匀采样。在（可能很差的）R中：

C <- 4
nvec <- 20
x <- rexp(nvec * (C - 1))
x <- x/(sum(x) + rexp(1))
A1 <- matrix(0, nvec, C)
A1[1:nvec, 2:C] <- matrix(x, nvec, C - 1)
A1 <- apply(A1, 2, cumsum)
A1 <- apply(A1, 2, rev)
A1 <- t(apply(A1, 1, cumsum))
A1 <- A1 + (1 - A1[1:nvec, C]) * runif(nvec)
很好的分析，+1！是的，你说得对，我没有限制自己在“均匀”分布内（因为我没有看到OP问题中的硬性要求），但只是提出了一种构建矩阵（或一堆向量，如果按行读取矩阵）的随机性想法。 - ThomasIsCoding
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接