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>
> 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-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
> > 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 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.
