[R] Inconsistent computation of an integral

[William Dunlap]

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


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of Aurélien Philippot
> Sent: Thursday, December 19, 2013 10:38 AM
> To: R-help at 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 at 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.