在不调用 polym() 的情况下，我如何计算它将 return 的交互次数？

Question

如何计算互动次数 poly 将 return？

如果我有两个变量，那么相互作用的数量 poly 将 return 在度数的函数中给出：

degree <- 2

dim(poly(rnorm(10), rnorm(10), degree = degree))[2]

这等同于：

(degree^2+3*degree)/2

有没有根据度数和变量的数量来计算交互的数量（如果我使用两个以上）？

Answer 1

组合的数学结果

假设您有 p 个变量，与 d 度相关的交互次数的计算公式为：

fd <- function (p, d) {
  k <- choose(p, d)
  if (d > 1) k <- k + p * sum(choose(p-1, 0:(d-2)))
  return(k)
  }

函数 poly（在本例中实际上是 polym），具有 p 个输入变量和一个 degree = D，将从 degree = 1 开始构建交互到 degree = D。所以下面的函数计算它：

fD <- function (p, D) {
   if (D < 1) return(0)
   component <- sapply(1:D, fd, p = p)
   list(component = component, ncol = sum(component))
   }

条目component给出了从1到D每个度数的交互次数，ncol分量给出了总交互次数。

快速测试：

a <- runif(50)
b <- runif(50)
c <- runif(50)
d <- runif(50)
X <- poly(a, b, c, d, degree = 3)
ncol(X)
# 34

fD(4, 3)
# component
# [1] 4  10  20
#
# ncol
# [1] 34

R 是怎么做到的？

polym 源代码的前几行解释了 R 如何解决这个问题。首先调用 expand.grid 来获取所有可能的交互，然后调用 rowSums 来计算所有可用交互的程度。最后，应用过滤器以仅保留度数在 1 和 D.

之间的交互项

Answer 2

三年多后，我不得不处理次数 >=3 的多项式。不幸的是，@李哲源解决方案对于大于 3 的度数失败。但是，我可以构建两个解决方案：

展开网格解决方案

此方法模拟 polym 原始行为，这对我们的目的来说不是很优雅，但它是一个自然的基准。

expand_grid_solution <- function(nd, degree){
  z <- do.call(expand.grid, c(rep.int(list(0:degree), nd), 
                              KEEP.OUT.ATTRS = FALSE))
  s <- rowSums(z)
  ind <- 0 < s & s <= degree
  z <- z[ind, , drop = FALSE]
  s <- s[ind]
  return(length(s))
}

结合重复解法

combination_with_repetition <- function(n, r){
  factorial(r+n-1)/(factorial(n-1)*factorial(r))
} 

poly_elements <- function(n, d) {
  x <- sapply(1:d, combination_with_repetition, n = n) 
  return(sum(x))
}

快速测试：

mapply(expand_grid_solution, c(2,2,2,3,3,3,4), c(2,3,4,2,3,4,4))
#[1]  5  9 14  9 19 34 69

mapply(poly_elements, c(2,2,2,3,3,3,4), c(2,3,4,2,3,4,4))
#[1]  5  9 14  9 19 34 69

在不调用 polym() 的情况下，我如何计算它将 return 的交互次数？

Without calling polym(), how can I count the number of interactions it will return?

r

polynomials