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

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在:

在这里， k大约是 1.23 :

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

但是如果我想针对许多不同的整体事件概率解决这个问题，我如何在 R 中有效地做到这一点？
我看过函数 SMfzero包裹内 NLRoot ，但我仍然不清楚如何实现它。
编辑
到目前为止，我已经对解决方案进行了基准测试。上玩具资料 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

创建于 2021-04-28 由 reprex package (v1.0.0)