[R] Double Integration

[Ravi Varadhan]

Here is another approach using `integrate':

fvec = function(x, y) sapply(x, function(z, y) exp(-0.549451*(z^2+y^2-0.6*z*y)), y=y)
 
gvec = function(x) sapply(x, function(y) integrate(fvec, lower=-1.51, upper=2.696,  subdivisions=10000, rel.tol=1.e-08, y=y)$val)

> integrate(gvec, lower=1.98, upper=3.54, subdivisions=10000, rel.tol=1.e-06) 
0.1376102 with absolute error < 1.5e-15
> 

Ravi. 

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Hans W Borchers <hwborchers at googlemail.com>
Date: Thursday, July 1, 2010 2:47 pm
Subject: Re: [R] Double Integration
To: r-help at stat.math.ethz.ch


> Sarah Sanchez <sarah_sanchez09 <at> yahoo.com> writes:
>  > 
>  > Dear R helpers
>  > 
>  > I am working on the Bi-variate Normal distribution probabilities.
>  > I need to double integrate the following function
>  > (actually simplified form of bivariate normal distribution)
>  > 
>  > f(x, y) = exp [ - 0.549451 * (x^2 + y^2 - 0.6 * x * y) ]
>  > 
>  > where 2.696 < x < 3.54 and -1.51 < y < 1.98
>  > 
>  > I need to solve something like
>  > 
>  > INTEGRATE (2.696 to 3.54) dx INTEGRATE [(-1.51 to 1.98)] f(x, y) dy
>  > 
>  > I have referred to stats::integrate but it deals with only one 
> variable. 
>  > 
>  > This example appears in Internal Credit Risk Model by Michael Ong
>  > (page no. 160).
>  
>  There has been the same request last month, the answer is still valid:
>  
>      library(cubature)
>      fun <- function(x) exp(-0.549451*(x[1]^2+x[2]^2-0.6*x[1]*x[2]))
>      adaptIntegrate(fun, c(-1.51,1.98), c(2.696,3.54), tol=1e-8)
>  
>      # $integral
>      # [1] 0.1376102
>     # $error
>      # [1] 1.372980e-09
>      # ...
>  
>  Hans Werner
>  
>  >
>  > Kindly guide.
>  > 
>  > Regards
>  > 
>  > Sarah
>  >
>  
>  ______________________________________________
>  R-help at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.