[R] Numerical Integration

[Julio Rojas]

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?

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