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)}

Answer 2

Following @IaroslavDomin's solutions, but constructing the coefficients for this particular case by hand, then using polyroot():

Here's a sequence of three functions (compute individual coeffs, put them together into a vector, find positive real roots):

## 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)
}

Try it out for a couple of cases:

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)

If you don't want to be clever, then brute force is perfectly adequate (56 microseconds on my machine when length(p) is 52):

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))
}

I started to run benchmarks and gave up; I don't like the format of microbenchmark output and the approximate solution is too fast for rbenchmark::benchmark() to time accurately. In any case, one run of bsoln() with length(p)==52 takes on the order of 50 microseconds, so you're going to have to run this a whole bunch of times before speed becomes problematic ...

Answer 3

Another option would be to just search for a root on a segment without specifically calculating polynomial coefficients. This can be done e.g. with the uniroot function.

Only one not-so-trivial thing we need to do here is to specify the segment. k is obviously >=0 - so that would be the left point. Then we know that all the k*p values should be probabilities, hence <=1. Therefore k <= 1/max(p) - that's the right point.

And so the code is:

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

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

This works pretty fast even for 1000 probabilities.