[R] sampling from the unoform distrubuton over a convex hull
(Ted Harding)
ted.harding at nessie.mcc.ac.uk
Sat Mar 24 20:26:21 CET 2007
On 24-Mar-07 14:00:44, Ted Harding wrote:
> [...]
> 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.
Well, I've written some rather crude code to implement the
above. Even allowing for possible "editorial" improvements,
the more I look at it the more I think there may be a better
way! (Of doing it directly, I mean, raher than a rejection
method).
Still, it seems to work (and quite slickly) ...
## Needed for rdiric() to sample from Dirichlet
## Better to write one's own rdiric() for this application!
library(VGAM)
## Make some random points, get their CH and centroid, and draw them
X<-cbind(rnorm(10),rnorm(10))
plot(X[,1],X[,2],pch="+",col="green")
H<-chull(X)
C<-colMeans(X[H,])
points(C[1],C[2],pch="+",col="red")
## Draw the CH and the triangulation
H<-c(H,H[1]) ## To close the contour
K<-length(H)-1 ## No of triangles
lines(X[H,1],X[H,2],col="blue")
for(i in (1:K)){lines(c(X[H[i],1],C[1]),c(X[H[i],2],C[2]),col="red")}
## Set up the sampling by triangles
As<-numeric(K)
Ns<-numeric(K)
## Get the areas of the triangles
for(i in (1:(K))){
V1<-X[H[i],]
V2<-X[H[i+1],]
V3<-C
As[i]<-abs(det(rbind(V1-V3,V2-V3)))
}
## As illustration, 1000 points over the CH
N<-2000
## How many points to go in eavh triangle
Cuts<-cumsum(As)/sum(As)
R<-runif(N)
Ns[1]<-sum(R<=Cuts[1])
for(i in (2:K)){Ns[i]<-sum(R<=Cuts[i])-sum(R<=Cuts[i-1])}
## Uniform distribution over each triangle
for(i in (1:(K))){
V1<-X[H[i],]
V2<-X[H[i+1],]
V3<-C
## The Dirichlet sample:
T<-rbind(X[H[i],],X[H[i+1],],C)
D<-rdiric(Ns[i],c(1,1,1))
Z<-D%*%T
points(Z[,1],Z[,2],pch="+",col="blue")
}
Any advance on the above??
Best wishes,
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: 19:26:18
------------------------------ XFMail ------------------------------
More information about the R-help
mailing list