It is of course the question of defintion but the routine your suggest for random distribution of points in the sphere is not uniformally  random in my understanding since   density of the points you obtains  depends on the distance to the  pole (is bigger near poles and smaller near equator). 

I suggest theirfore another routine:

x1:= 2* random(n)-1; y1:= 2*random(n)-1;  z1:=2*random(n)-1; r:=sqrt(x1^2+y1^2+z1^2);

x:=x1/r; y:=y1/r; z:= z1/r;

The points (x,y,z) lie on the sphere and are randomly distributed.

Now, in order to obtain the points randomly distributed withing any shape on the sphere, you

can add the command

check_whether_the_point_xyz_is_within_shape(x,y,z)

, and through  out points such that they this command returns 0.

The subroutine check_whether_the_point_xyz_is_within_shape depends of course on how the shape looks like. I do not know what is an ellips on a spere for you, but in the case the ellips is say a connected component  interior of the intersection of the  elliptic cone $\{(x,y,z) : 2*x^2+ y^2 = z^2\}$ with the sphere, the subroutine could look like

if  2*x^2+ y^2 = z^2 & z>1 then return 1 else return 0:

The recommended approach would be using a Poisson distribution which is then culled by the space constraints (surface of sphere in an ellipse shaped window).

There are several optimizations that can be added on top of this, but the baseline algorithm is given by Bridson (in his SIGGRAPH 07 Sketch), see attachment.

http://people.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf