efficiently approximate real solution for polynomial function

Question

I want to efficiently solve a degree-7 polynomial in k.

For example, with the following set of 7 unconditional probabilities,

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

the overall event probability is approximately 25% :

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

And if I want to approximate a constant k to proportionally change all elements of p so that the overall event probability is now approximately 30%, I can do this using an equation solver (such as Wolfram Alpha), which may use Newton's method or bisection to approximate k in:

here, k is approximately 1.23:

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

But what if I want to solve this for many different overall event probabilities, how can I efficiently do this in R?

I've looked at the function SMfzero in the package NLRoot, but it's still not clear to me how I can achieve it.

EDIT I've benchmarked the solutions so far. On the toy data p above:

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

NB, I'm not sure if its fair to compare the approximation_fun with the others but I did ask for an approximate solution so it does meet the brief.

The real problem is a degree-52 polynomial in k. Benchmarking on the real data:

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

This can be solved with the polynom library.

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

^{Created on 2021-04-28 by the reprex package (v1.0.0)}