[R] numerical integration of a bivariate function

[Berend Hasselman]

On 22-04-2013, at 15:04, Hicham Mezouara <hicham_dess at yahoo.fr> wrote:

> 
> hello
> I work on
> the probabilities of bivariate normal distribution. I need
> integrate  the
> following function.
> f (x, y) = exp [- (x ^ 2 + y ^ 2 + x * y)] with - ∞ ≤ x ≤
> 7.44 and - ∞ ≤ y ≤ 1.44   , either software R or  matlab Version R 2009a

This discussion in a far and distant past may help you:

https://stat.ethz.ch/pipermail/r-help/2004-October/059494.html

I tried this

f <- function(x,y) exp(- (x ^ 2 + y ^ 2 + x * y))

integrate(function(y) {
  sapply(y, function(y) {
    integrate(function(x) f(x,y), -Inf, 7.44)$value
  })
}, -Inf, 1.44)
# 3.486503 with absolute error < 0.00013

library(cubature)              
h <- function(z) f(z[1],z[2])
adaptIntegrate(h,c(-10,-10),c(1.44,7.44))
# $integral
# [1] 3.486496
# 
# $error
# [1] 3.415381e-05
# 
# $functionEvaluations
# [1] 4267
# 
# $returnCode
# [1] 0

adaptIntegrate doesn't like -Inf so I used -10.


Berend