[R] sampling from the unoform distrubuton over a convex hull

[(Ted Harding)]

On 24-Mar-07 05:43:06, Ranjan Maitra wrote:
> Dear list,
> 
> Does anyone have a suggestion (or better still) code for sampling from
> the uniform distribution over the convex hull of a set of points?
> 
> Many thanks and best wishes,
> Ranjan

I was curious if anyone would come up with some ready-made efficient
code for this problem -- it cannot be the first time it has been
addressed! But if Duncan Murdoch doesn't, that reduces the probability
that such code is already available in R!

Duncan's suggestion of a rejection method for uniform points in an
enclosing rectangle will probably be efficient, provided the convex
hull occupies a good proportion of the rectangle. So I would suggest
for this method that a preliminary transformation be made, rotating
onto principal axes of the convex hull. This would avoid the situation
where the convex hull is a narrow cigar-shape at 45 degrees. At the
end, transform back to the original coordinates.

Another possible approach (again in two dimensions) could be based
on the following.

First, if A, B, C are the vertices of a triangle, let (w1,w2,w3)
be sampled from the 3-variate Dirichlet distribution with index
("shape") parameters (1,1,1). Then the weighted sum

   w1*A + w2*B + w3*C

is uniformly distributed over the triangle.[1] (This does not
generalise to planar convex hulls with more than three vertices).

Next, triangulate the convex hull (e.g. joining each vertex to the
centroid of the convex hull, getting K triangles, say), and calculate
the area of each triangle. Then, to sample N points, partition N at
random into K parts with probabilities proportional to the areas.
For instance, by cutting

sort(runif(N))

at the breakpoints given by

cumsum(A)/sum(A)

where A is a vector of areas

Then, for each component Nj of Ns, sample Nj points uniformly
over triangle j using the Dirichlet method above.

[1] This can be seen geometrically: (w1,w2,w3) is uniformly
distributed over the triangle T1 with vertices (1,0,0), (0,1,0)
and (0,0,1). Given any other triangle T2, by rotating T1 in
space and projecting orthogonally onto the plane containing T2,
you can match 2 (and therefore all 3) of the angles of T1 with
the angles of T2. Then dilate the projection of T1 until it
is congruent with T2. Since equal areas project orthogonally
onto equal areas, a uniform distribution on T1 projects into
a uniform on the projection of T1, therefore on T2.


PS I could envisage the above approach being generalised
to more than 2 dimensions. In 3 dimensions, for instance,
since the Dirchlet(1,1,1,1) is uniform on the simplex with
vertices (1,0,0,0), (0,1,0,0), (0,0,1,0), (0,0,0,1) we
similarly have that for a general simplex with vertices
A,B,C,D the point w1*A + w2*B + w3*C + w4*D is uniformly
distributed in the simplex.

But this requires "simplectification" of the convex hull,
(e.g. joining the vertices of its outer faces to its centroid)
so it gets more complicated, so no doubt Duncan's proposal
wins on simplicity, if not on efficiency (since the more
dimensions, the greater the proportion of points rejected).

Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 24-Mar-07                                       Time: 14:00:40
------------------------------ XFMail ------------------------------