[R] Numerical Integration

[William Dunlap]

> -----Original Message-----
> From: Julio Rojas [mailto:jcredberry at ymail.com] 
> Sent: Friday, December 18, 2009 9:06 AM
> To: William Dunlap; r-help at r-project.org
> Subject: RE: [R] Numerical Integration
> 
> Thanks a lot William. I'm sorry about the syntax problem. I 
> was working at the same time with R code and a document, so I 
> copied from the wrong document.
> 
> I saw that I had a problem with the conditions of the 
> integrand. A ">" instead of a ">=" on the second one. Again, thanks.
> 
> One last question: Is there a way to use "approx" as the integrand?

Try
   integrand2 <- approxfun(x=c(-1e20,fx,Inf),y=c(0,0,1,1,0,0))
instead of your
   Vector(integrand)

The -1e20 should logically be -Inf, but approxfun produces
bad code in that case so the resulting function produces NaN
instead of 0 for x<min(fx).

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com 

> 
> Best regards.
> 
> Julio
> 
> 
> --- El vie 18-dic-09, William Dunlap <wdunlap at tibco.com> escribió:
> 
> > De: William Dunlap <wdunlap at tibco.com>
> > Asunto: RE: [R] Numerical Integration
> > A: "Julio Rojas" <jcredberry at ymail.com>
> > Fecha: viernes, 18 diciembre, 2009, 4:03 pm
> > The immediate problem is that
> > integrand(0.5) is NULL
> > because none of the if clauses is TRUE if x==fx[2].
> > 
> > If you restructured the code such that it had no return
> > statements and instead assigned things to a variable
> > and returned that variable at the end it would be more
> > apparent what happened: you would get a 'variable not
> > found' error when returning.
> >    i <- function(x) {
> >       if (x<fx[1]) retval <- 0
> >       else if (x>fx[1] &&
> > x<=fx[2]) retval <- 1
> >       else if (x>fx[2]) retval <- 0
> >       retval
> >    }
> >    > i(.4)
> >    [1] 1
> >    > i(.3)
> >    Error in i(0.3) : object 'retval' not
> > found
> > 
> > When you convert your code for production use, look
> > at the ifelse and approx functions for doing this
> > sort of thing.  They are much faster than
> > Vector(yourIntegrand).
> > 
> > Also, when writing for help from R-help, please use
> > R syntax when defining variables, like
> >    fx<-c(0.3,0.5,0.5,0.6)
> > and not some other language, like
> >    fx<-[0.3,0.5,0.5,0.6]
> > If we can copy and paste your example into R instead
> > of translating the syntax by hand we will
> > be much more inclined to try to fix things up.
> > 
> > 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 Julio Rojas
> > > Sent: Friday, December 18, 2009 7:02 AM
> > > To: r-help at r-project.org
> > > Subject: [R] Numerical Integration
> > > 
> > > Dear @ll. I have to calculate numerical integrals for
> > 
> > > triangular and trapezoidal figures. I know you can
> > calculate 
> > > the exactly, but I want to do it this way to learn how
> > to 
> > > proceed with more complicated shapes. The code I'm
> > using is 
> > > the following:
> > > 
> > > integrand<-function(x) {
> > >    print(x)
> > >    if(x<fx[1]) return(0)
> > >    if(x>=fx[1] && x<fx[2])
> > return((x-fx[1])/(fx[2]-fx[1]))
> > >    if(x>fx[2] && x<=fx[3])
> > return(1)
> > >    if(x>fx[3] && x<=fx[4])
> > return((x-fx[4])/(fx[3]-fx[4]))
> > >    if(x>fx[4]) return(0)
> > > 
> > >  }
> > > 
> > > fx<-data[i,j,]
> > > reltol<-1e-07
> > >
> > integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions
> > > =200)$value
> > > 
> > > It works for most cases, but then, I tried for the
> > triangle 
> > > fx<-[0.3,0.5,0.5,0.6] and the following error was
> > presented:
> > > 
> > > >
> > integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions=200)
> > > [1] 0.5
> > > [1] 0.01304674
> > > [1] 0.9869533
> > > [1] 0.06746832
> > > [1] 0.9325317
> > > [1] 0.1602952
> > > [1] 0.8397048
> > > [1] 0.2833023
> > > [1] 0.7166977
> > > [1] 0.4255628
> > > [1] 0.5744372
> > > [1] 0.002171418
> > > [1] 0.9978286
> > > [1] 0.03492125
> > > [1] 0.9650787
> > > [1] 0.1095911
> > > [1] 0.8904089
> > > [1] 0.2186214
> > > [1] 0.7813786
> > > [1] 0.3528036
> > > [1] 0.6471964
> > > Error in integrate(Vectorize(integrand), 0, 1, rel.tol
> > = reltol,
> > > subdivisions = 200) :
> > >  evaluation of function gave a result of wrong
> > type
> > > 
> > > Does anybody know what happened? Thanks in advance!!!
> > > 
> > > 
> > >       
> > >
> > ______________________________________________________________
> > > ______________________
> > > ¡Obtén la mejor experiencia en la web!
> > > Descarga gratis el nuevo Internet Explorer 8. 
> > > http://downloads.yahoo.com/ieak8/?l=e1
> > > 
> > > ______________________________________________
> > > 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.
> > > 
> > 
> 
> 
>       
> ______________________________________________________________
> ______________________
> ¡Obtén la mejor experiencia en la web!
> Descarga gratis el nuevo Internet Explorer 8. 
> http://downloads.yahoo.com/ieak8/?l=e1
> 
>