[R] integration over a simplex

[Duncan Murdoch]

On 7/10/2007 6:57 AM, Robin Hankin wrote:
> Hello
> 
> The excellent adapt package integrates over multi-dimensional
> hypercubes.
> 
> I want to integrate over a multidimensional simplex.  Has anyone
> implemented such a thing in R?
> 
> I can transform an n-simplex to a hyperrectangle
> but the Jacobian is a rapidly-varying (and very lopsided)
> function and this is making adapt() slow.
> 
> [
> A \dfn{simplex} is an n-dimensional analogue of a triangle or  
> tetrahedron.
> It is the convex hull of (n+1) points in an n-dimensional Euclidean  
> space.
> 
> My application is a variant of the Dirichlet distribution:
> With p~D(a), if length(p) = n+1 then the requirement that
> all(p>0) and sum(p)=1 mean that the support of the
> Dirichlet distribution is an n-simplex.

I don't know what shape of simplex you're working with, but I believe 
the subset of an n-cube with coordinates ordered x[1] < x[2] < ... < 
x[n] is a simplex, and the cube can be tiled with n! of those, by 
permuting the order of the coordinates.  So if your function is smooth 
enough at the edges you might be able to map n! copies of it onto a 
cube, and use adapt to integrate over that.

That is:  if f() is your function, defined on 0 < x[1] < x[2] < ... < 
x[n] < 1, define g <- function(x) f(sort(x)), and the integral you want 
is (1/n!) times the integral of g over the unit cube.

Duncan Murdoch