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.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
