我有两个矩阵'time'和'scaleFactor'。"时间"矩阵中的每一行代表一个人,列代表某个事件发生的时间。该事件更改放大因子,该因子在"比例因子"矩阵中可用。例如,对于人 1,时间间隔 (1,3) 中的因子为 1.2。
我想创建一个矩阵 C,在每个整数时间点都有因子,其中第一列是 time=0。我已经编写了生成矩阵"C"的代码,但我想知道是否可以避免"for"循环,因为实际矩阵非常大。
scaleFactor <- matrix(c(1,1.1,1.2,1.4,
1,1.3,1.4,1.6,
1,1.2,1.6,2.1),nrow = 3,ncol = 4, byrow = T)
times <- matrix(c(0,1,3,99,
0,2,5,99,
0,1,4,99),nrow = 3,ncol = 4,byrow = T)
> scaleFactor
[,1] [,2] [,3] [,4]
[1,] 1 1.1 1.2 1.4
[2,] 1 1.3 1.4 1.6
[3,] 1 1.2 1.6 2.1
> times
[,1] [,2] [,3] [,4]
[1,] 0 1 3 99
[2,] 0 2 5 99
[3,] 0 1 4 99
C <- matrix(0,nrow = 3,ncol = 6)
for (i in 1:ncol(C)){
indices <- max.col(i-1<=times,'first')
C[,i] <- scaleFactor[cbind(1:3,indices)]
}
> C
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 1.1 1.2 1.2 1.4 1.4
[2,] 1 1.3 1.3 1.4 1.4 1.4
[3,] 1 1.2 1.6 1.6 1.6 2.1
非常感谢您的帮助。
编辑:我后来意识到"times"矩阵可以有非整数时间值。对于这种情况,akrun建议的方法有效。谢谢!
times2 <- matrix(c(0,1.2,3.6,99,
0,2.1,5.3,99,
0,1,4,99),nrow = 3,ncol = 4,byrow = T)
> times2
[,1] [,2] [,3] [,4]
[1,] 0 1.2 3.6 99
[2,] 0 2.1 5.3 99
[3,] 0 1.0 4.0 99
# Matrix C would be -
> C
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 1.1 1.2 1.2 1.4 1.4
[2,] 1 1.3 1.3 1.4 1.4 1.4
[3,] 1 1.2 1.6 1.6 1.6 2.1
这是一种rep
纯Vectorized
代码许可的方法
- 我们
rep
C 中的列序列(减去 1)与"时间"的length
- 对
t
的"时代"进行比较 - 转换为
matrix
以创建dim
属性 - 应用
max.col
以获取每行最大值的列索引 - 然后,使用行索引
cbind
- 根据 5 中的行/列索引获取相应的"比例因子"值
- 使用
[]
将输出分配回 C,以便保留矩阵属性
C[] <- scaleFactor[cbind(1:3, max.col(matrix(rep(seq_len(ncol(C)) - 1,
each = length(times)) <= c(t(times)), ncol = ncol(times),
byrow = TRUE), "first"))]
-输出
C
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 1.1 1.2 1.2 1.4 1.4
[2,] 1 1.3 1.3 1.4 1.4 1.4
[3,] 1 1.2 1.6 1.6 1.6 2.1
我们可以使用rep
和asplit
的技巧(将数组拆分为向量列表,在这里很方便)。
C <- t(
mapply(rep, asplit(scaleFactor,1),
times = asplit(cbind(1, t(apply(times, 1, diff))), 1))
)
dim(C)
# [1] 3 100
C[,1:8]
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# [1,] 1 1.1 1.2 1.2 1.4 1.4 1.4 1.4
# [2,] 1 1.3 1.3 1.4 1.4 1.4 1.6 1.6
# [3,] 1 1.2 1.6 1.6 1.6 2.1 2.1 2.1
它的稳定性取决于最后一列times
相等(以便所有人都迭代相同数量的时间段)。
这通过首先将times
(即某些内容发生变化的时间段)转换为"持续时间"(每列中的天数)来工作,其中
cbind(1, t(apply(times, 1, diff)))
# [,1] [,2] [,3] [,4]
# [1,] 1 1 2 96
# [2,] 1 2 3 94
# [3,] 1 1 3 95
从那里,我们使用这些数字作为rep
的times=
参数。为了保持尺寸正确,我们为每个人独立执行此操作,这就是为什么我们需要使用asplit(., 1)
按行拆分矩阵:
asplit(scaleFactor, 1)
# [[1]]
# [1] 1.0 1.1 1.2 1.4
# [[2]]
# [1] 1.0 1.3 1.4 1.6
# [[3]]
# [1] 1.0 1.2 1.6 2.1
然后,我们将比例因子与相应的时间配对成rep
:
head(
rep(c(1, 1.1, 1.2, 1.4), times = c(1, 1, 2, 96)),
8)
# [1] 1.0 1.1 1.2 1.2 1.4 1.4 1.4 1.4
为了同时在两个矩阵的所有行中有效地迭代它,我们asplit
矩阵并将它们输入mapply
,它通过第一行scaleFactor
的第一行"持续时间"(times
差),第二行与第二行,依此类推。