[R] random uniform sample of points on an ellipsoid (e.g. WG

[Ranjan Maitra]

Hi,

Sorry for being a late entrant to this thread, but let me see if I understand the problem.

The poster wants to sample from an ellipsoid. Let us call this ellipsoid X'\Gamma X - d^2= 0.

There is no loss in assuming that the center is zero, otherwise the same can be done.

Let us consider the case when Gamma = I first. 

Then, let X \sim N_p(0, I) (any radially symmetric distribution will
do here), then d*X/||X|| is uniform on the sphere of radius d. 

How about imitating the same?

Let X \sim N_p(0, \Sigma), where \Sigma = \Gamma^{-1} then X restricted to X'\Gamma X = d^2 gives the required uniform density on the ellipsoid.

How do we get this easily? I don't think rejection works or is even necessary.

Isn't d*X / ||\Gamma^{1/2}X|| enough? Here \Gamma^{1/2} is the square root matrix of \Gamma.

Note that any distribution of the kind f(X'\Gamma X) would work, but the multivariate Gaussian is a convenient tool, for which two-lines of R code should be enough.


Many thanks and best wishes,
Ranjan




On Sat, 24 Feb 2007 13:45:56 -0000 (GMT) (Ted Harding) <ted.harding at nessie.mcc.ac.uk> wrote:

> [Apologies if this is a repeated posting for you. Something seems
>  to have gone amiss with my previous attempts to post this reply,
>  as seen from my end]
> 
> On 22-Feb-07 Roger Bivand wrote:
> > On 21 Feb 2007, Russell Senior wrote:
> > 
> >> 
> >> I am interested in making a random sample from a uniform distribution
> >> of points over the surface of the earth, using the WGS84 ellipsoid as
> >> a model for the earth.  I know how to do this for a sphere, but would
> >> like to do better.  I can supply random numbers, want latitude
> >> longitude pairs out.
> >> 
> >> Can anyone point me at a solution?  Thanks very much.
> >> 
> > 
> > http://www.csit.fsu.edu/~burkardt/f_src/random_data/random_data.html
> > 
> > looks promising, untried.
> 
> Hmmm ... That page didn't seem to be directly useful, since
> on my understanding of the code (and comments) listed under
> "subroutine uniform_on_ellipsoid_map(dim_num, n, a, r, seed, x)"
> "UNIFORM_ON_ELLIPSOID_MAP maps uniform points onto an ellipsoid."
> in
> 
> http://www.csit.fsu.edu/~burkardt/f_src/random_data/random_data.f90
> 
> it takes points uniformly distributed on a sphere and then
> linearly transforms these onto an ellipsoid. This will not
> give unform density over the surface of the ellipsoid: indeed
> the example graph they show of points on an ellipse generated
> in this way clearly appear to be more dense at the "ends" of
> the ellipse, and less dense on its "sides". See:
> 
> http://www.csit.fsu.edu/~burkardt/f_src/random_data/
> uniform_on_ellipsoid_map.png
> [all one line]
> 
> Indeed, if I understand their method correctly, in the case
> of a horizontal ellipse it is equivalent (modulo rotating
> the result) to distributing the points uniformly over a circle,
> and then stretching the circle sideways. This will preserve
> the vertical distribution (so at the two ends of the major axis
> it has the same density as on the circle) but diluting the
> horizontal distribution (so that at the two ends of the minor
> axis the density isless than on the circle).
> 
> I did have a notion about this, but sat on it expecting that
> someone would come up with a slick solution -- which hasn't
> happened yet.
> 
> For the application you have in hand, uniform distribution
> over a sphere is a fairly close approximation to uniform
> distriobution over the ellipspoid -- but not quite.
> 
> But a rejection method, applied to points uniform on the sphere,
> can give you points uniform on the ellipsoid and, because of
> the close approximation of the sphere to the ellipsoid, you
> would not be rejecting many points.
> 
> The outline strategy I had in mind (I haven't worked out details)
> is based on the following.
> 
> Consider a point X0 on the sphere, at radial distance r0 from
> the centre of the sphere (same as the centre of the ellipsoid).
> Let the radius through that point meet the ellipsoid at a point
> X1, at radial distance R1.
> 
> Let dS0 be an element of area at X0 on the sphere, which projects
> radially onto an element of area dS1 on the ellipsoid. You want
> all elements dS1 of equal size to be equally likely to receive
> a random point.
> 
> Let the angle between the tangent plane to the ellipsoid at X1,
> and the tangent plane to the sphere at X0, be phi.
> 
> The the ratio of areas dS1/dS0 is R(X0), say, where
> 
>   R(X0) = dS1/dS0 = r1^2/(r0^2 * cos(phi))
> 
> and the smaller this ratio, the less likely you want a point
> u.d. on the sphere to give rise to a point on the ellipsoid.
> 
> Now define an acceptance probability P(X0) by
> 
>   P(X0) = R(X0)/sup[R(X)]
> 
> taking the supremum over X on the sphere. Then sample points X0
> unformly on the sphere, accepting each one with probability
> P(X0), and continue sampling until you have the number of
> points that you need.
> 
> Maybe someone has a better idea ... (or code for the above!)
> 
> Ted.
> 
> 
> --------------------------------------------------------------------
> E-Mail: (Ted Harding) <ted.harding at nessie.mcc.ac.uk>
> Fax-to-email: +44 (0)870 094 0861
> Date: 24-Feb-07                                       Time: 13:45:50
> ------------------------------ XFMail ------------------------------
> 
> ______________________________________________
> R-help at 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.
>