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

Sat Feb 24 23:49:25 CET 2007

On 24-Feb-07 Ranjan Maitra wrote:
> 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

As I understand Rusell Senior's original query, he wants to
generate a uniform distribution on the surface of the ellipsoid,
not over its interior:

  "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."

Your solution, and the solution given in Roger Bivand's reference

Section "UNIFORM_IN_ELLIPSOID_MAP maps uniform points into an ellipsoid"
in:
http://www.csit.fsu.edu/~burkardt/f_src/random_data/random_data.f90

is valid for uniform points in the interior of an ellipsoid, I think.

But, since it is a linear transformation, it is not valid for the
points on the surface, as I explain for the case of an ellipse
(in particular, it results in higher linear density along the
circumference of an ellipse near the ends of the major axis
than near the ends of the minor axis).

It is also the method adopted for points on the surface in the
Section "UNIFORM_ON_ELLIPSOID_MAP maps uniform points onto an ellipsoid",
and I think this is wrong, as I explained.

Ted.

> 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.
>>
> 
> ______________________________________________
> 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.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 24-Feb-07                                       Time: 22:49:09
------------------------------ XFMail ------------------------------