Re: [R] sampling from the unoform distrubuton over a convex hull
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) <[EMAIL PROTECTED]>
Fax-to-email: +44 (0)870 094 0861
Date: 25-Mar-07 Time: 11:26:58
-- 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.
Re: [R] sampling from the unoform distrubuton over a convex hull
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.
Re: [R] sampling from the unoform distrubuton over a convex hull
Duncan Murdoch wrote: > > Suggestion: Find a rectangular region containing the hull, and sample > uniformly there. Accept points that don't expand the hull of the > original points. > This is a feasible solution only in lower dimensions; the area of the convex hull can become exponentially small relative to its external rectangle. For example, if the convex hull is approximately spherical, the relation of the (hyper)volumes are (time to summon Wikipedia)... http://en.wikipedia.org/wiki/Hypersphere#Hyperspherical_volume ... Volume.ratio <- pi^(n/2) * 2^(-n) / gamma(n/2 + 1) Even for n = 5 this would be 0.16 The hypervolume of a _hyperpyramidal_ convex hull would be even worse, as it would be 1/n! (= 1/gamma(n+1)) of the hypervolume of the hypercube. Alberto Monteiro __ 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.
Re: [R] sampling from the unoform distrubuton over a convex hull
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) <[EMAIL PROTECTED]>
Fax-to-email: +44 (0)870 094 0861
Date: 24-Mar-07 Time: 14:00:40
-- 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.
Re: [R] sampling from the unoform distrubuton over a convex hull
The convhulln() function in the "geometry" package provides an R interface to
qhull for computing n-dimensional convex hull.
Andy
From: [EMAIL PROTECTED] on behalf of Duncan Murdoch
Sent: Sat 3/24/2007 7:52 AM
To: Ranjan Maitra
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] sampling from the unoform distrubuton over a convex hull
[Broadcast]
On 3/24/2007 1:43 AM, 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?
Are you talking about two dimensional points, or higher dimensions? The
suggestion below works for any dimension, but the actual code is
2-dimensional. I don't know if there's an equivalent of chull available
for higher dimensions.
Suggestion: Find a rectangular region containing the hull, and sample
uniformly there. Accept points that don't expand the hull of the
original points.
For example:
rhull <- function(n,x) {
boundary <- x[chull(x),]
xlim <- range(boundary[,1])
ylim <- range(boundary[,2])
boundary <- rbind(c(NA,NA), boundary) # add space for new test point
result <- matrix(NA, n, 2)
for (i in 1:n) {
repeat {
boundary[1,1] <- runif(1, xlim[1], xlim[2])
boundary[1,2] <- runif(1, ylim[1], ylim[2])
if ( !(1 %in% chull(boundary)) ) {
result[i,] <- boundary[1,]
break
}
}
}
result
}
x <- matrix(rnorm(20), ncol=2)
plot(x, cex=2, col="red")
sample <- rhull(1000, x)
points(sample)
Duncan Murdoch
__
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.
--
Notice: This e-mail message, together with any attachments,...{{dropped}}
__
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.
Re: [R] sampling from the unoform distrubuton over a convex hull
On 3/24/2007 1:43 AM, 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?
Are you talking about two dimensional points, or higher dimensions? The
suggestion below works for any dimension, but the actual code is
2-dimensional. I don't know if there's an equivalent of chull available
for higher dimensions.
Suggestion: Find a rectangular region containing the hull, and sample
uniformly there. Accept points that don't expand the hull of the
original points.
For example:
rhull <- function(n,x) {
boundary <- x[chull(x),]
xlim <- range(boundary[,1])
ylim <- range(boundary[,2])
boundary <- rbind(c(NA,NA), boundary) # add space for new test point
result <- matrix(NA, n, 2)
for (i in 1:n) {
repeat {
boundary[1,1] <- runif(1, xlim[1], xlim[2])
boundary[1,2] <- runif(1, ylim[1], ylim[2])
if ( !(1 %in% chull(boundary)) ) {
result[i,] <- boundary[1,]
break
}
}
}
result
}
x <- matrix(rnorm(20), ncol=2)
plot(x, cex=2, col="red")
sample <- rhull(1000, x)
points(sample)
Duncan Murdoch
__
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.
[R] sampling from the unoform distrubuton over a convex hull
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 __ 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.
