Date: Thu, 30 Nov 95 11:07:49 CST
From: rusin (Dave Rusin)
To: chrisman@ucdmath.ucdavis.edu
Subject: Re: Random points on ellipsoid
Newsgroups: sci.math

It's not easy even to place points on an ellipse. There is a method, which
I will describe for ellipses -- it generalizes to ellipsoids, but I'll
let you decide first if this looks like a method you consider appropriate.

Begin with the parameterization  (x,y)= (a cos theta, b sin theta)  of
the ellipse  (x/a)^2+(y/b)^2=1. If you randomly select  theta  in [0,2 pi],
you'll get points all over the ellipse, but they will tend to be concentrated
on the long sides of the ellipse and sparse on the pointy ends. What we
do is to use just this process but select points in the transformed domain
_not_ uniformly but rather according to some other distribution. We just
have to figure out what that distribution is.

Your goal is to have it be that the likelihood a point is put on a
certain arc of the ellipse is proportional to the length of that
arc. Thus we need
	Pr(  (a cos theta, b sin theta) in arc ) = K * length(arc)
or, writing the arc as the set of points from  (x0,y0) to (x1, y1) (where I
will assume the entire arc lies above the x-axis, say) you need
	Pr(  x0 <= a cos theta <= x1 ) = K * integral(sqrt(1+(y')^2) dx)
where the integral is taken from  x0  to  x1. 

The left hand side may also be seen as an integral, as follows. Let 
f(theta) be the probability density function for our ultimate distribution.
That is, we will eventually choose  theta  in  [0, 2 pi]  not uniformly
but such that the probability that  theta  lies in an interval [u1, u2]
is  integral( f(t) dt, t=u1..u2). In terms of this, the left side is
clearly the integral   integral(f(t) dt, t=arccos(x0/a)..arccos(x1/a)).

Now, with a fixed  x0,  both integrals may be viewed as functions of  x1.
We want these functions to be equal. They are equal when x1=x0, so it suffices
to have their derivatives be equal. By the Fundamental Theorem of Calculus,
these derivatives are, respectively
    f(arccos(x/a))*(d/dx)(arccos(x/a)) = -f(arccos(x/a))/[a.sqrt(1-(x/a)^2)]
and
    K* sqrt(1+(dy/dx)^2) = K * sqrt([1-c.(x/a)^2]/[1-(x/a)^2])
where  c=1-(b/a)^2. Setting  u=arccos(x/a)  (i.e., x=a cos(u)) we then
express the equality of these integrals by writing
	f(u)=-a.sqrt(1-cos(u)^2)*K*sqrt([1-c.cos(u)^2]/[1-cos(u)^2])
	    =K'.sqrt(1-c.cos(u)^2)
where  K'  will the constant necessary to make the integral of  f  be  1.

So that's your probability density function: f(u)=K' sqrt(1-c.cos(u)^2),
where  c=1-(b/a)^2  and  K' = 1/ integral( sqrt( 1-c. cos(u)^2), u=0..2 pi).
This integral must be evaluated numerically -- there is no closed form
antiderivative, unless for example  c=0 (ellipse is a circle).
Once  u  is selected, you put your point  (a cos(u), b sin(u)) on the ellipse.
These points will be uniformly distributed.

(There are several methods of generating random numbers in an interval
according to some predetermined non-uniform distribution. I think  Knuth
has a section on this.)

The procedure for ellipsoids and indeed any manifold is to choose a 
parameterization and select points in the parameter domain according to
a pdf which is not uniform but rather adjusted by the jacobian (to
account for the change in area by the parameterization). 

dave