有效逼近多项式函数的实数解

Question

我想高效地求解 k 的 7 次多项式。

例如，对于以下一组 7 个无条件概率，

p <- c(0.0496772, 0.04584501, 0.04210299, 0.04026439, 0.03844668, 0.03487194, 0.03137491)

总体事件概率约为 25% :

> 1 - prod(1 - p)
[1] 0.2506676

如果我想近似常数 k 以按比例更改 p 的所有元素，以便现在总体事件概率约为 30%，我可以使用方程求解器 (例如 Wolfram Alpha)，它可能使用牛顿法或二分法来近似 k in:

这里，k大约是1.23:

> 1 - prod(1 - 1.23*p)
[1] 0.3000173

但是如果我想针对许多不同的整体事件概率解决这个问题，我如何才能在 R 中有效地做到这一点呢？

我查看了包NLRoot中的函数SMfzero，但我仍然不清楚如何实现它。

编辑到目前为止，我已经对解决方案进行了基准测试。上面的玩具数据p：

Unit: nanoseconds
              expr     min      lq      mean  median      uq      max neval
 approximation_fun     800    1700    3306.7    3100    4400    39500  1000
       polynom_fun 1583800 1748600 2067028.6 1846300 2036300 16332600  1000
      polyroot_fun  596800  658300  863454.2  716250  792100 44709000  1000
         bsoln_fun   48800   59800   87029.6   85100  102350   613300  1000
        find_k_fun   48500   60700   86657.4   85250  103050   262600  1000

注意，我不确定将 approximation_fun 与其他比较是否公平，但我确实要求了一个近似的解决方案，所以它确实符合要求。

真正的问题是 k 的 52 次多项式。对真实数据进行基准测试：

Unit: microseconds
              expr     min       lq       mean   median       uq     max neval
 approximation_fun     1.9     3.20     7.8745     5.50    14.50    55.5  1000
       polynom_fun 10177.2 10965.20 12542.4195 11268.45 12149.95 80230.9  1000
         bsoln_fun    52.3    60.95    91.4209    71.80   117.75   295.6  1000
        find_k_fun    55.0    62.80    90.1710    73.10   118.40   358.2  1000

Answer 1

这可以通过 polynom 库解决。

library(polynom)
library(purrr)

p <- runif(3, 0, 1)
p
#> [1] 0.1072518 0.5781922 0.3877427
# Overall probability
1 - prod(1 - p)
#> [1] 0.7694434

# Target overall probability
target_op <- 0.3

# calculate polynomial to solve for k
poly_list <- p %>% 
  map(~polynomial(c(1, -.))) %>% 
  as.polylist()

# List of linear polynomials to be multiplied:
poly_list
#> [[1]]
#> 1 - 0.1072518*x 
#> 
#> [[2]]
#> 1 - 0.5781922*x 
#> 
#> [[3]]
#> 1 - 0.3877427*x

# we want to solve this polynomial
poly <- 1 - prod(poly_list) - target_op
poly
#> -0.3 + 1.073187*x - 0.3277881*x^2 + 0.02404476*x^3
roots <- solve(poly)
good_roots <- 
  roots %>% 
  # keep only real values
  keep(~Im(.) == 0) %>% 
  Re() %>% 
  # only positive
  keep(~.>0)

good_roots
#> [1] 0.1448852

k <- good_roots[[1]]

1 - prod(1 - k*p)
#> [1] 0.3

^{由 reprex package (v1.0.0)}

于 2021-04-28 创建

Answer 2

遵循@IaroslavDomin 的解决方案，但手动构建此特定情况的系数，然后使用 polyroot():

这是一个包含三个函数的序列（计算各个系数，将它们放在一个向量中，找到正实根）：

## construct ith binomial coefficients: the sum of the products 
## of all i-element combinations
bcoef <- function(p,i) {
    sum(apply(combn(p,i),2,prod))
}
## compute all binomial coefficients and put them together
## into the vector of coeffs for 1-prod(1-k*p)
mypoly <- function(p,target=0.3) {
    c(-target,-1*sapply(seq_along(p), bcoef, p =-p))
}
## compute real positive solutions
soln <- function(p, target=0.3) {
    roots <- polyroot(mypoly(p))
    roots <- Re(roots[abs(Im(roots))<1e-16])
    roots <- roots[roots>0]
    if (length(roots)>1) warn(">1 solution")
    return(roots)
}

尝试几个案例：

p1 <- c(0.1072518,0.5781922, 0.3877427)
s1 <- soln(p1)
1-prod(1-s1*p1)

p2 <- c(0.0496772, 0.04584501, 0.04210299, 0.04026439, 0.03844668, 0.03487194, 0.03137491)
s2 <- soln(p2)
1-prod(1-s2*p2)

如果不想耍小聪明，那么蛮力就足够了（length(p) 为 52 时，我的机器上是 56 微秒）：

bsoln <- function(p, target=0.3) {
    f <- function(k) { (1-prod(1-k*p)) - target }
    return(uniroot(f, c(0,20))$root)
}
asoln <- function(p, target=0.3) {
    return(- log(1 - target) / sum(p))
}

我开始运行基准测试并放弃了；我不喜欢 microbenchmark 输出的格式，而且近似解太快 rbenchmark::benchmark() 无法准确计时。在任何情况下，bsoln() 和 length(p)==52 中的一个运行需要大约 50 微秒，所以你将不得不运行在速度出现问题之前重复了很多次...

Answer 3

另一种选择是只搜索线段上的根而不专门计算多项式系数。这可以做到，例如使用 uniroot 函数。

我们在这里只需要做一件非常重要的事情就是指定段。 k 显然是 >=0 - 所以那将是左边的点。然后我们知道所有 k*p 值应该是概率，因此 <=1。因此 k <= 1/max(p) - 这是正确的观点。

所以代码是：

find_k <- function(p, taget_op) {
  f <- function(x) 1 - prod(1 - x*p) - target_op
  max_k <- 1/max(p)
  res <- uniroot(f, c(0, max_k))
  res$root
}

p <- runif(1000, 0, 1)
target_op <- 0.3
k <- find_k(p, target_op)
k
#> [1] 0.000710281

1 - prod(1 - k*p)
#> [1] 0.2985806

^{由 reprex package (v1.0.0)}

于 2021-04-29 创建

即使是 1000 个概率，这也非常快。

有效逼近多项式函数的实数解

efficiently approximate real solution for polynomial function

r

probability

polynomials

newtons-method

bisection