-3

I have to find a double integral as a function of $ \alpha, \beta $ and $ \gamma $.

enter image description here

f(a,b,c) = \int_0^1 \int_0^{1-y1} \gamma(a+b+c)/\gamma(a)/\gamma(b)/\gamma(c) * y1^(a+1) * y2^(b-1) * (1-y1-y2)^(c-1) dy_2 dy_1

The output of this double integral should be an R function that takes three arguments a, b and c. How can I implement this in R?

MrFlick
  • 195,160
  • 17
  • 277
  • 295
user67275
  • 1
  • 9
  • 38
  • 64
  • 2
    http://stackoverflow.com/questions/8913603/calculating-double-integrals-in-r-quickly?rq=1 Didn't even have to search, this should have appeared when you typed in that title and warned you not to make this question – rawr May 25 '14 at 21:31
  • rawr/ I don't think your criticism is reasonable. My question is different from that. – user67275 May 26 '14 at 00:24

2 Answers2

1

The domain of integration is the triangle with vertices (0,0), (0,1) and (1,0). The SimplicialCubature package is the way to use for such a domain.

a=1
b=2
c=3
K <- gamma(a+b+c)/gamma(a)/gamma(b)/gamma(c)
f <- function(x){
  K * x[1]^(a+1)*x[2]^(b-1)*(1-x[1]-x[2])^(c-1)
}
S <- cbind(c(0,0),c(0,1),c(1,0)) # domain of integration (triangle)

> library(SimplicialCubature)
> adaptIntegrateSimplex(f, S)
$integral
[1] 0.04761905

$estAbsError
[1] 4.761905e-14

$functionEvaluations
[1] 32

$returnCode
[1] 0

$message
[1] "OK"
Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225
-1

In the past, I've used the cubature package. It has worked for me in cases similar to yours, and I'm sure you can use it!

user3616949
  • 53
  • 1
  • 2
  • 17