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

Sun Feb 25 01:37:39 CET 2007

"My" method is for the surface, not for the interior. The constraint d*X/|| \Gamma^{-1/2}X ||ensures the constraint, no? The uniformity is ensured by the density restricted to satisfy the constraint which makes it a constant.

Ranjan

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

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