问题描述
我想有效地求解 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
解决方法
这可以通过 polynom 库解决。
library(polynom)
library(purrr)
p <- runif(3,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 年 4 月 28 日创建
,遵循@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 微秒的数量级，因此您将不得不运行很多次才能提高速度变得有问题...
另一种选择是只搜索段上的根而不专门计算多项式系数。这可以完成，例如使用 uniroot 函数。
在这里我们需要做的只有一件非常重要的事情就是指定段。 k 显然 >=0 - 所以这将是左边的点。然后我们知道所有 k*p 值都应该是概率，因此 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,max_k))
res$root
}
p <- runif(1000,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 年 4 月 29 日创建
即使对于 1000 个概率，这也非常快。