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

[(Ted Harding)]

On 24-Mar-07 19:26:21, Ted Harding wrote:
> On 24-Mar-07 14:00:44, Ted Harding wrote:
>> [...]
> [...]
> 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) ...
> 
> 
And now I've done what I should have donein the first place:
the code for rdiric() in VGAM is "stand-alone" and can simply
be copied, so:

## library(VGAM) gives the following code for function rdiric()
rdiric<-function(n, shape, dimension = NULL) 
{
    dimension = if (!is.numeric(dimension)) length(shape)
    shape = rep(shape, len = dimension)
    ans = if (is.R()) 
        rgamma(n * dimension, rep(shape, rep(n, dimension)))
    else rgamma(n * dimension, rep(shape, each = n))
    dim(ans) = c(n, dimension)
    ans = ans/apply(ans, 1, sum)
    ans
}

## 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<-1000
 
## 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")
}

Ted. 

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