[R] Numerical Integration

[Xiaofan Cao]

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@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] numerical integration problem

[przeszczepan]

Hi,

I have got problems integrating the following function using "integrate": 

lambdat<-function(t){
tempT<-T[k,][!is.na(T[k,])]#available values from k-th row of matrix T
tempJ<-J[k,][!is.na(J[k,])]

hg<-length(tempT[tempT<=t & tempJ==0])#counts observations satisfing the 
conditions
ag<-length(tempT[tempT<=t & tempJ==1])

lambdaXY[hg+1,ag+1]#takes  values from a 10x10 matrix
}

I keep receiving this message:

1: longer object length
is not a multiple of shorter object length in: tempT <= t 
2: longer object length
is not a multiple of shorter object length in: tempT <= t & tempJ == 0 

What I suspect is that the "integrate" function submits the whole vector of 
points at which the integral is to be evaluated at once. For my function to be 
integrated it would rather have to be evaluated at each point after another in 
a loop of some kind.

Can you think of a way to solve this problem without me having to write the 
integrating procedure from scratch (I have no idea about FORTRAN and this is 
what the "integrate" description refers to)?

Thank you.

Kind Regards,
Lukasz Szczepanski
Student

__
R-help@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

Re: [R] Numerical Integration

[Doran, Harold]

You might try the statmod package which provides nodes and weights for
gaussian quadrature. 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Xiaofan Cao
> Sent: Wednesday, November 08, 2006 12:43 PM
> To: r-help@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@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@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.

Re: [R] Numerical Integration

[Gabor Grothendieck]

You could approximate it with splines and then integrate that:

x <- 0:10/10
f <- function(x) 1 + x + x^2 + x^3 + x^4
y <- f(x)
fs <- splinefun(x, y)
integrate(fs, 0, 1)
integrate(f, 0, 1)

On 11/8/06, Xiaofan Cao <[EMAIL PROTECTED]> wrote:
> 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@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@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.

Re: [R] Numerical Integration

[Ravi Varadhan]

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: [EMAIL PROTECTED]

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

 




-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Doran, Harold
Sent: Wednesday, November 08, 2006 1:04 PM
To: Xiaofan Cao; r-help@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: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Xiaofan Cao
> Sent: Wednesday, November 08, 2006 12:43 PM
> To: r-help@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@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@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@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.

Re: [R] Numerical Integration

[Xiaofan Cao]

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: [EMAIL PROTECTED]
>
> Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
> 
> 
>
> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Doran, Harold
> Sent: Wednesday, November 08, 2006 1:04 PM
> To: Xiaofan Cao; r-help@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: [EMAIL PROTECTED]
>> [mailto:[EMAIL PROTECTED] On Behalf Of Xiaofan Cao
>> Sent: Wednesday, November 08, 2006 12:43 PM
>> To: r-help@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@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@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@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.

Re: [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.61

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: [EMAIL PROTECTED]

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

 




-Original Message-
From: Xiaofan Cao [mailto:[EMAIL PROTECTED] 
Sent: Wednesday, November 08, 2006 4:21 PM
To: Ravi Varadhan
Cc: 'Doran, Harold'; r-help@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: [EMAIL PROTECTED]
>
> Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
>

> 
>
> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Doran, Harold
> Sent: Wednesday, November 08, 2006 1:04 PM
> To: Xiaofan Cao; r-help@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: [EMAIL PROTECTED]
>> [mailto:[EMAIL PROTECTED] On Behalf Of Xiaofan Cao
>> Sent: Wednesday, November 08, 2006 12:43 PM
>> To: r-help@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@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@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@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.

Re: [R] numerical integration problem

[Sundar Dorai-Raj]

przeszczepan wrote:
> Hi,
> 
> I have got problems integrating the following function using "integrate": 
> 
> lambdat<-function(t){
> tempT<-T[k,][!is.na(T[k,])]#available values from k-th row of matrix T
> tempJ<-J[k,][!is.na(J[k,])]
> 
> hg<-length(tempT[tempT<=t & tempJ==0])#counts observations satisfing the 
> conditions
> ag<-length(tempT[tempT<=t & tempJ==1])
> 
> lambdaXY[hg+1,ag+1]#takes  values from a 10x10 matrix
> }
> 
> I keep receiving this message:
> 
> 1: longer object length
> is not a multiple of shorter object length in: tempT <= t 
> 2: longer object length
> is not a multiple of shorter object length in: tempT <= t & tempJ == 
> 0 
> 
> What I suspect is that the "integrate" function submits the whole vector of 
> points at which the integral is to be evaluated at once. For my function to 
> be integrated it would rather have to be evaluated at each point after 
> another in a loop of some kind.
> 

You suspect correctly. Best to read ?integrate too.

> Can you think of a way to solve this problem without me having to write the 
> integrating procedure from scratch (I have no idea about FORTRAN and this is 
> what the "integrate" description refers to)?
> 

Just put a "for"-loop in your function to iterate over t.

   n <- length(t)
   hg <- ag <- vector("numeric", n)
   for(i in seq(n)) {
 hg[i] <- length(tempT[tempT <= t[i] & tempJ == 0])
 ag[i] <- length(tempT[tempT <= t[i] & tempJ == 1])
   }

I doubt this will work because integrate is expecting a vector of 
n=length(t) from lambdat. The last line of the function returns a nxn 
matrix. Please submit data to run the function plus your call to 
integrate in any subsequent postings.

HTH,

--sundar

> Thank you.
> 
> Kind Regards,
> Lukasz Szczepanski
> Student
> 
> __
> R-help@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

__
R-help@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

Re: [R] numerical integration problem

[Duncan Murdoch]

On 6/29/2006 7:38 AM, Sundar Dorai-Raj wrote:
> 
> przeszczepan wrote:
>> Hi,
>> 
>> I have got problems integrating the following function using "integrate": 
>> 
>> lambdat<-function(t){
>> tempT<-T[k,][!is.na(T[k,])]#available values from k-th row of matrix T
>> tempJ<-J[k,][!is.na(J[k,])]
>> 
>> hg<-length(tempT[tempT<=t & tempJ==0])#counts observations satisfing the 
>> conditions
>> ag<-length(tempT[tempT<=t & tempJ==1])
>> 
>> lambdaXY[hg+1,ag+1]#takes  values from a 10x10 matrix
>> }
>> 
>> I keep receiving this message:
>> 
>> 1: longer object length
>> is not a multiple of shorter object length in: tempT <= t 
>> 2: longer object length
>> is not a multiple of shorter object length in: tempT <= t & tempJ == 
>> 0 
>> 
>> What I suspect is that the "integrate" function submits the whole vector of 
>> points at which the integral is to be evaluated at once. For my function to 
>> be integrated it would rather have to be evaluated at each point after 
>> another in a loop of some kind.
>> 
> 
> You suspect correctly. Best to read ?integrate too.
> 
>> Can you think of a way to solve this problem without me having to write the 
>> integrating procedure from scratch (I have no idea about FORTRAN and this is 
>> what the "integrate" description refers to)?
>> 
> 
> Just put a "for"-loop in your function to iterate over t.

Or use Vectorize().

vlambdat <- Vectorize(lambdat)

should give a function that can be passed to integrate(), assuming that 
lambdat works when given a length 1 vector as input.

Duncan Murdoch


> 
>n <- length(t)
>hg <- ag <- vector("numeric", n)
>for(i in seq(n)) {
>  hg[i] <- length(tempT[tempT <= t[i] & tempJ == 0])
>  ag[i] <- length(tempT[tempT <= t[i] & tempJ == 1])
>}
> 
> I doubt this will work because integrate is expecting a vector of 
> n=length(t) from lambdat. The last line of the function returns a nxn 
> matrix. Please submit data to run the function plus your call to 
> integrate in any subsequent postings.
> 
> HTH,
> 
> --sundar
> 
>> Thank you.
>> 
>> Kind Regards,
>> Lukasz Szczepanski
>> Student
>> 
>> __
>> R-help@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
> 
> __
> R-help@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

__
R-help@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

[R] numerical integration of x,y data

[Rajarshi Guha]

Hi, I was looking for a routine that would do numerical integration
using x,y data rather than integrating a function. 

Searching CRAN showed me the sfsmisc package, which contains
integrate.xy. However the documentation mentions that its not good for
noisy data and plans to implement the Romberg method.

Is this the package that is generally used to perform this type of
integration or are there other routines available?

Thanks,


---
Rajarshi Guha <[EMAIL PROTECTED]> 
GPG Fingerprint: 0CCA 8EE2 2EEB 25E2 AB04 06F7 1BB9 E634 9B87 56EE
---
Say it with flowers - give her a triffid

__
R-help@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