[R] Problem with uniroot.all

[Duncan Murdoch]

On 13-12-26 6:26 AM, Aurélien Philippot wrote:
> Dear R experts,
> I am trying to find all the solutions of an equation. Here is an example:
>
>
> integrand1<- function(x){1/x*dnorm(x)}
> integrand2<- function(x){1/(2*x-50)*dnorm(x)}
>
> integrand3<- function(x,C){
>    cd<- 1/(x+C)
>    return(cd)
> }
>
> res<- function(C){
>    ce<-integrate(integrand1, lower=1, upper=50)$value+ integrate(integrand2,
> lower=50, upper=100)$value-integrate(integrand3, C=C,lower=1,
> upper=100)$value
>    return(ce)
> }
>
> I want to find all the roots of res.
>
> First, I used uniroot to ensure that a solution existed, and it works.
> uniroot(res, c(1,1000))
>
> But,  there might be other roots, and I wanted to use the function uniroot.
> library(rootSolve)
> uniroot.all(res, c(1,1000))
>
> I obtain an error message:
> "Error in integrate(integrand3, C = C, lower = 1, upper = 100) :
>    evaluation of function gave a result of wrong length
> In addition: Warning message:
> In x + C : longer object length is not a multiple of shorter object length"
>
> I don't understand what is wrong? Thank you very much for any suggestion!

You can often use traceback() to find where the error came from:

3: integrate(integrand3, C = C, lower = 1, upper = 100) at #2
2: f(xseq, ...)
1: uniroot.all(res, c(1, 1000))

Not very informative, but notice that on line 2 the argument to f is 
called xseq.  If you call debug(res) and try again, it will break on the 
first call to res, and you can print C:

Browse[2]> C
   [1]    1.00   10.99   20.98   30.97   40.96   50.95   60.94   70.93 
  80.92   90.91  100.90  110.89  120.88  130.87  140.86  150.85
  [17]  160.84  170.83  180.82  190.81  200.80  210.79  220.78  230.77 
240.76  250.75  260.74  270.73  280.72  290.71  300.70  310.69
  [33]  320.68  330.67  340.66  350.65  360.64  370.63  380.62  390.61 
400.60  410.59  420.58  430.57  440.56  450.55  460.54  470.53
  [49]  480.52  490.51  500.50  510.49  520.48  530.47  540.46  550.45 
560.44  570.43  580.42  590.41  600.40  610.39  620.38  630.37
  [65]  640.36  650.35  660.34  670.33  680.32  690.31  700.30  710.29 
720.28  730.27  740.26  750.25  760.24  770.23  780.22  790.21
  [81]  800.20  810.19  820.18  830.17  840.16  850.15  860.14  870.13 
880.12  890.11  900.10  910.09  920.08  930.07  940.06  950.05
  [97]  960.04  970.03  980.02  990.01 1000.00

Aha!  uniroot.all() passes a vector of values to the function, whereas 
uniroot passes values one at a time.  You need to make sure res can 
handle a vector of C values.  The easiest method is to use Vectorize to 
vectorize it:

res <- Vectorize(res)

This isn't very fast, but in your function, it's fine.

Duncan Murdoch

P.S.  You might want to write to the maintainer of uniroot.all to point 
out that the documentation doesn't mention this difference from uniroot.