如何对R中的“系列”向量进行向量化操作

Question

我在 R 中有一个函数，它接受一个标量和一个向量作为参数，对它们执行一些操作，返回一个值。

给定一个“系列”标量（这里是向量mya）和一个“系列”向量（这里是矩阵myv），我如何向量化对myf的调用，以便每个元素在mya 中与myv 中的对应向量一起使用？

mya = 1:3
myv = matrix(1:30, 10, 3)

myf = function(a, v) {
  return(sum(a / (a/v + 1)))
}

sapply(1:3, function(x) {myf(mya[x], myv[,x])})
# [1]  7.980123 17.649590 26.809440

所以上面我想避免循环 sapply 操作直接做类似的事情：

myf(mya, myv)
# [1] 49.37443   <- Here I would like 3 values

这里最大的问题是性能：在我的实际情况下，mya 和 myv 将分别有超过 10e6 个值或向量，而myf 要复杂得多。

Answer 1

首先，您的myv 可能被组织为一系列向量，每个向量一列；许多工具最好将其转换为list 的向量。

asplit(myv, 2)
# [[1]]
#  [1]  1  2  3  4  5  6  7  8  9 10
# [[2]]
#  [1] 11 12 13 14 15 16 17 18 19 20
# [[3]]
#  [1] 21 22 23 24 25 26 27 28 29 30

基础 R

sapply/lapply 是单个向量/列表，因为 mapply/Map 是 n。

Map(myf, mya, asplit(myv , 2))
# [[1]]
# [1] 7.980123
# [[2]]
# [1] 17.64959
# [[3]]
# [1] 26.80944
mapply(myf, mya, asplit(myv , 2))
# [1]  7.980123 17.649590 26.809440

tidyverse

参数的顺序不同，而不是单独的参数，它需要将所有参数都放在 list 本身中。

purrr::pmap(list(mya, asplit(myv , 2)), myf)
# [[1]]
# [1] 7.980123
# [[2]]
# [1] 17.64959
# [[3]]
# [1] 26.80944
purrr::pmap_dbl(list(mya, asplit(myv , 2)), myf)
# [1]  7.980123 17.649590 26.809440

---

给定 cmets 的替代方法。

这种方法确实是矢量化的，但对函数进行了一点解构。

colSums(t(mya / (mya / t(myv) + 1)))
# [1]  7.980123 17.649590 26.809440

要达到这一点，需要识别transpose 的位置，这是必要的。我将从一些已知点开始：

mya[1] / myv[,1] + 1
#  [1] 2.000000 1.500000 1.333333 1.250000 1.200000 1.166667 1.142857 1.125000 1.111111 1.100000

为了用矩阵（而不仅仅是向量）模拟这一点，我们可以尝试

(mya / myv + 1)
#           [,1]     [,2]     [,3]
#  [1,] 2.000000 1.181818 1.142857
#  [2,] 2.000000 1.250000 1.045455
#  [3,] 2.000000 1.076923 1.086957
#  [4,] 1.250000 1.142857 1.125000
#  [5,] 1.400000 1.200000 1.040000
#  [6,] 1.500000 1.062500 1.076923
#  [7,] 1.142857 1.117647 1.111111
#  [8,] 1.250000 1.166667 1.035714
#  [9,] 1.333333 1.052632 1.068966
# [10,] 1.100000 1.100000 1.100000

但如果您注意到，mya 与 myv 的划分是按列划分的，因此它正在扩展到

c(mya[1] / myv[1,1], mya[2] / myv[2,1], mya[3] / myv[3,1], mya[1] / myv[4,1], ...)

我们希望它被转置的地方。好的，所以我们转置它，使myv 的行 垂直划分。

(mya / t(myv) + 1)[1,]
#  [1] 2.000000 1.500000 1.333333 1.250000 1.200000 1.166667 1.142857 1.125000 1.111111 1.100000

这样更好。现在我们需要为下一步做同样的事情。这给我们带来了

t(mya / (mya / t(myv) + 1))
#            [,1]     [,2]     [,3]
#  [1,] 0.5000000 1.692308 2.625000
#  [2,] 0.6666667 1.714286 2.640000
#  [3,] 0.7500000 1.733333 2.653846
#  [4,] 0.8000000 1.750000 2.666667
#  [5,] 0.8333333 1.764706 2.678571
#  [6,] 0.8571429 1.777778 2.689655
#  [7,] 0.8750000 1.789474 2.700000
#  [8,] 0.8888889 1.800000 2.709677
#  [9,] 0.9000000 1.809524 2.718750
# [10,] 0.9090909 1.818182 2.727273

因为您想对每个 mya 值求和。知道我们在mya 中有三列并且我们看到三列，我们可能会推断我们需要对每一列求和。我们可以通过经验证明：

sum(mya[1] / (mya[1] / myv[,1] + 1))
# [1] 7.980123
colSums(t(mya / (mya / t(myv) + 1)))
# [1]  7.980123 17.649590 26.809440

但实际上，当我们无法转置并对行求和时，我们不需要 tranpose 然后对列求和 :-)

rowSums(mya / (mya / t(myv) + 1))
# [1]  7.980123 17.649590 26.809440