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) <[EMAIL PROTECTED]> Fax-to-email: +44 (0)870 094 0861 Date: 24-Mar-07 Time: 19:26:18 ------------------------------ XFMail ------------------------------ ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.