R help - Dec 2013 - 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]]

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.