[R] Numerical Integration

[Ravi Varadhan]

Here is a simple example using the "bkde" (binned kernel density estimator)
from the KernSmooth package that shows how to use the trapezoidal
integration:

> library(KernSmooth)
> data(geyser, package="MASS")
> x <- geyser$duration
> est <- bkde(x, grid = 101, bandwidth=0.25)
> trap.rule(est$x,est$y)
[1] 0.999961

The answer from the simple trapezoidal rule is pretty good with an error of
less than 1.e-04 (number of grid points used is 101, with equal spacing).
You could use higher order rules (e.g. Simpson's), but you can get
reasonable accuracy with Trapezoidal rule by just increasing the number of
grid points.

Ravi.
----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 

----------------------------------------------------------------------------
--------

-----Original Message-----
From: Xiaofan Cao [mailto:cao at stat.colostate.edu] 
Sent: Wednesday, November 08, 2006 4:21 PM
To: Ravi Varadhan
Cc: 'Doran, Harold'; r-help at stat.math.ethz.ch
Subject: RE: [R] Numerical Integration

Hi Ravi and Harold,

Thanks for the input. I'm using trapezoidal rule and like to know if 
there's other alternatives. This f(x) is the kernel density estimator and 
thus we can get an estimate of f(x) at any given x in theory.

Thanks again,
Martha

  On Wed, 8 
Nov 2006, Ravi Varadhan wrote:

> Can you get an estimate of f(x) at any given x?  If so, the Gaussian
> quadrature methods will work, but not otherwise since f(x) must be known
at
> all the nodes.  A rough approximation to the integral can be obtained
using
> the trapezoidal rule.  Here is a simple function to do that:
> trap.rule <- function(x,y) sum(diff(x)*(y[-1]+y[-length(y)]))/2
>
> However, the use of the word "knots" seems to indicate that some sort of
> spline is being fit to the data.  Martha - can you provide more
information
> about your function f(x)?
>
> Ravi.
>
----------------------------------------------------------------------------
> -------
>
> Ravi Varadhan, Ph.D.
>
> Assistant Professor, The Center on Aging and Health
>
> Division of Geriatric Medicine and Gerontology
>
> Johns Hopkins University
>
> Ph: (410) 502-2619
>
> Fax: (410) 614-9625
>
> Email: rvaradhan at jhmi.edu
>
> Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
>
----------------------------------------------------------------------------
> --------
>
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Doran, Harold
> Sent: Wednesday, November 08, 2006 1:04 PM
> To: Xiaofan Cao; r-help at stat.math.ethz.ch
> Subject: Re: [R] Numerical Integration
>
> You might try the statmod package which provides nodes and weights for
> gaussian quadrature.
>
>> -----Original Message-----
>> From: r-help-bounces at stat.math.ethz.ch
>> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Xiaofan Cao
>> Sent: Wednesday, November 08, 2006 12:43 PM
>> To: r-help at stat.math.ethz.ch
>> Subject: [R] Numerical Integration
>>
>> Hi everyone,
>>
>> I'm trying to integrate f(x) over x where f(x) does not have
>> a close form but only numerical values at certurn knots of x.
>> Is there a way that I can use any generical R function (such
>> as integrate) or any package to do so?
>>
>> Thanks! I appreciate your time.
>>
>> Best Regards,
>> Martha Cao
>>
>> ______________________________________________
>> R-help at stat.math.ethz.ch 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.
>>
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>