[R] integration over a simplex
Duncan Murdoch
murdoch at stats.uwo.ca
Tue Jul 10 19:25:06 CEST 2007
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
More information about the R-help
mailing list