You should try to match up values of integrand() instead of the value of its integral on a single integral if you want to see if the function is being computed correctly. integrate(f) calls f(x) where x is a vector of values (typically 21 values). It expects that f() is vectorized: that f(x[i]) == f(x)[i] for i in seq_along(x). When you put max(x) into your function it messes that up: > integrand(0.5) [1] -0.07377961 > integrand(1.0) [1] -0.05514443 > integrand(c(0.5, 1.0)) # expect the above two numbers in output vector [1] -0.05106583 -0.05514443
Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com From: Aurélien Philippot [mailto:aurelien.philip...@gmail.com] Sent: Friday, December 20, 2013 8:59 AM To: William Dunlap; R-help@r-project.org Subject: Re: [R] Inconsistent computation of an integral Thanks William. I was convinced by pmax, until I played with the following example today. I tried to reproduce a computation from a paper. Here is the code: P<- function(x) { ab<-100*exp((0.0435-0.0269-0.5*0.3315^2)*4.3122+x*sqrt(4.3122)*0.3315) return(ab) } integrand<- function(x){ cd<- -1/2*((0.094+1.1465)*exp(0.0435*4.3122)+0.29/100*exp(0.0269*4.3122)*P(x)+0.89/100*max(P(x)-88.254,0))^(-2)*pnorm(x) return(cd) } Solution<- integrate(integrand, lower=-Inf, upper=Inf) Solution The above code gives: -0.1366377 The paper reports: -0.1349 If I use pmax instead of max, I get: -0.1965606, which is much worse. It may look like small discrepancies, but it makes a big difference to me. I am still very puzzled by these discrepancies. 2013/12/19 William Dunlap <wdun...@tibco.com<mailto:wdun...@tibco.com>> I think you want to use pmax(x-50, 0), which returns a vector the length of x, instead of max(x-50,0), which returns a scalar. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com<http://tibco.com> > -----Original Message----- > From: r-help-boun...@r-project.org<mailto:r-help-boun...@r-project.org> > [mailto:r-help-boun...@r-project.org<mailto:r-help-boun...@r-project.org>] On > Behalf > Of Aurélien Philippot > Sent: Thursday, December 19, 2013 10:38 AM > To: R-help@r-project.org<mailto:R-help@r-project.org> > Subject: [R] Inconsistent computation of an integral > > Dear R experts, > I computed the same integral in two different ways, and find different > values in R. > The difference is due to the max function that is part of the integrand. In > the first case, I keep it as such, in the second case, I split it in two > depending on the values of the variable of integration. > > 1) First computation > > # Function g > g<- > function(x){1/(x*0.20*sqrt(10)*sqrt(2*pi))*exp(-0.5*((log(x/50)- > 0.1*10)/(0.20*sqrt(10)))^2)} > > ####### Function f1 > f1<- function(x) {1/(5000000+100000*x+10000*max(x-50,0))} > > integrand1<- function(x) { > out<- f1(x)*g(x) > return(out) > } > > i2<- integrate(integrand1, lower=0, upper=Inf )$value > > It gives me: i2= 3.819418e-08 > > 2) Second computation > I break the max function in two, depending on the values of the variable of > integration x (and I use the same density g as before): > > f11<- function(x) {1/(5000000+100000*x)} > f12<- function(x) {1/(5000000+100000*x+10000*(x-50))} > > > integrand11<- function(x) { > out<- f11(x)*g(x) > return(out) > } > > integrand12<- function(x) { > out<- f12(x)*g(x) > return(out) > } > > > i21<- integrate(integrand11, lower=0, upper=50 )$value > +integrate(integrand12, lower=50, upper=Inf)$value > > I get i21=5.239735e-08 > > The difference makes a huge difference for the computations I do. Does > anyone know where it comes from? > Thanks in advance! > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@r-project.org<mailto: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. [[alternative HTML version deleted]]
______________________________________________ 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.
