我正在寻找为 5 个变量创建一个矩阵,以便每个变量从seq(from = 0, to = 1, length.out = 500)
and中获取一个值rowSums(data) = 1
。
换句话说,我想知道如何创建一个矩阵来显示所有可能的数字组合以及每个row = 1
.
我正在寻找为 5 个变量创建一个矩阵,以便每个变量从seq(from = 0, to = 1, length.out = 500)
and中获取一个值rowSums(data) = 1
。
换句话说,我想知道如何创建一个矩阵来显示所有可能的数字组合以及每个row = 1
.
如果我理解正确,这至少可以带你走上正轨。
# Parameters
len_vec = 500 # vector length
num_col = 5 # number of columns
# Creating the values for the matrix using rational numbers between 0 and 1
values <- runif(len_vec*num_col)
# Creating matrix
mat <- matrix(values,ncol = num_col,byrow = T)
# ROunding the matrix to create only 0s and 1s
mat <- round(mat)
# Calculating the sum per row
apply(mat,1,sum)
这是一个使用循环的迭代解决方案。为您提供所有可能的数字排列,加起来为 1,它们之间的距离是 N 的倍数。这里的想法是取所有从 0 到 1 的数字(它们之间的距离是 N 的倍数),然后对于每个一个包括在新列中添加时不超过 1 的所有数字。冲洗并重复,除了在最后一次迭代中,您只需将完成行的数字添加到行的总和中。
就像人们在评论中指出的那样,如果你想要 N = 1/499*,它会给你一个非常大的矩阵。我注意到对于 N = 1/200,它已经花费了大约 2、3 分钟,所以对于 N = 1/499,它可能需要太长时间。
*seq(from = 0, to = 1, length.out = 500)
与seq(from = 0, to = 1, by = 1/499)
N = 1/2
M = 5
x1 = seq(0, 1, by = N)
df = data.frame(x1)
for(i in 1:(M-2)){
x_next = sapply(rowSums(df), function(x){seq(0, 1-x, by = N)})
df = data.frame(sapply(df, rep, sapply(x_next,length)))
df = cbind(df, unlist(x_next))
}
x_next = sapply(rowSums(df), function(x){1-x})
df = sapply(df, rep, sapply(x_next,length))
df = data.frame(df)
df = cbind(df, unlist(x_next))
> df
x1 unlist.x_next. unlist.x_next..1 unlist.x_next..2 unlist(x_next)
1 0.0 0.0 0.0 0.0 1.0
2 0.0 0.0 0.0 0.5 0.5
3 0.0 0.0 0.0 1.0 0.0
4 0.0 0.0 0.5 0.0 0.5
5 0.0 0.0 0.5 0.5 0.0
6 0.0 0.0 1.0 0.0 0.0
7 0.0 0.5 0.0 0.0 0.5
8 0.0 0.5 0.0 0.5 0.0
9 0.0 0.5 0.5 0.0 0.0
10 0.0 1.0 0.0 0.0 0.0
11 0.5 0.0 0.0 0.0 0.5
12 0.5 0.0 0.0 0.5 0.0
13 0.5 0.0 0.5 0.0 0.0
14 0.5 0.5 0.0 0.0 0.0
15 1.0 0.0 0.0 0.0 0.0
这正是包装partitions
的用途。基本上,OP 正在寻找总和为 499 的 5 个整数的所有可能组合。这可以通过以下方式轻松实现restrictedparts
:
system.time(combsOne <- t(as.matrix(restrictedparts(499, 5))) / 499)
user system elapsed
1.635 0.867 2.502
head(combsOne)
[,1] [,2] [,3] [,4] [,5]
[1,] 1.000000 0.000000000 0 0 0
[2,] 0.997996 0.002004008 0 0 0
[3,] 0.995992 0.004008016 0 0 0
[4,] 0.993988 0.006012024 0 0 0
[5,] 0.991984 0.008016032 0 0 0
[6,] 0.989980 0.010020040 0 0 0
tail(combsOne)
[,1] [,2] [,3] [,4] [,5]
[22849595,] 0.2024048 0.2004008 0.2004008 0.2004008 0.1963928
[22849596,] 0.2064128 0.1983968 0.1983968 0.1983968 0.1983968
[22849597,] 0.2044088 0.2004008 0.1983968 0.1983968 0.1983968
[22849598,] 0.2024048 0.2024048 0.1983968 0.1983968 0.1983968
[22849599,] 0.2024048 0.2004008 0.2004008 0.1983968 0.1983968
[22849600,] 0.2004008 0.2004008 0.2004008 0.2004008 0.1983968
而且由于我们处理的是数值,我们无法获得精确的精度,但是我们可以获得机器精度:
all(rowSums(combsOne) == 1)
[1] FALSE
all((rowSums(combsOne) - 1) < .Machine$double.eps)
[1] TRUE
有超过 2200 万条结果:
row(combsOne)
[1] 22849600