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.
]
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
Hi Robin, A Monte-Carlo approach could be attempted, if one could generate samples that are either uniformly distributed over the simplex. There is a small section in Luc Devroye's book (Generation of Non-uniform random deviates) on random uniform sampling from a simplex, if I remeber correctly. Another approach is importance sampling, where the sampling points have a characterized distribution. I have seen a technique called polyEDA, based on Gibbs sampling and truncated multivariate normal distribution. I had previously emailed the authors of this approach for the code, but haven't received a reply yet. You can google "polyEDA" for more info. I am interested in various computational problems related to polyhedra (e.g. enumeration of vertices, locating extrema, random sampling). I would appreciate if you'd keep me posted on how you solved this problem. Best, Ravi. ----- Original Message ----- From: Robin Hankin <r.hankin at noc.soton.ac.uk> Date: Tuesday, July 10, 2007 6:58 am Subject: [R] integration over a simplex To: RHelp help <r-help at stat.math.ethz.ch>> 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. > ] > > > -- > Robin Hankin > Uncertainty Analyst > National Oceanography Centre, Southampton > European Way, Southampton SO14 3ZH, UK > tel 023-8059-7743 > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code.
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