6 ms·

The easiest way to generate random points on a (hyper)sphere, at least with respect to most sampling APIs, is to sample from an isotropic normal distribution of

• josephg 11y ago
  This is not quite correct. Generating a random vector using three calls to rand() is the equivalent of picking a random point inside a cube. If you normalize the points, you're effectively squishing the corners of the cube in toward the center to make a sphere. Points near the corners of the containing cube will appear more frequently than the top, bottom and 'sides'. @ginko's comment is correct - you can fix the algorithm by throwing out any points that lie outside the sphere before normalizing.
  • mtinkerhess 11y ago
    I don't totally understand the math, but davmre is suggesting calling randn(), not rand() -- so sampling from a normal distribution instead of a uniform distribution.
    • darkmighty 11y ago
      The math is fairly simple. The probability distribution of a multivariate standard gaussian has a simple form that is f=a * exp(-x1^2+...+xN^2)=b * exp(-R^2), where a and b are some normalizing constants and R is the norm of the position , that is, it is obviously rotationally symmetric [1]. But that pdf is also the joint pdf of N i.i.d. gaussians, evident by decomposing f=a * exp(-x1^2) * ... * exp(-xN^2) [2], which is the joint pdf x1,...xN s.t. fx1=c * exp(-x1^2), ..., fxn=c * exp(-xN^2). [1] Since exp(-R^2) does not depend on direction but only on distance from the origin [2] The fact that f(x1,...xN)=f1 * ... * fN if x1,...xN are independent follows directly from the fact that P(A & B) = P(A)*P(B) if A and B are independent events.
  • T-hawk 11y ago
    Great point, you are correct. It's not enough for the distribution of the random vectors to be rotationally symmetric. It must be isotropic, looking the same in every direction from the origin, which a cube does not. (Nor does any other method deriving from a polyhedral solid, as I was thinking might have been possible, through something like selecting a random face from an icosahedron then a random point on that face.)
  • kasdkashdkas 11y ago
    No, OP is correct; he spoke about an isotropic normal distribution (randn as opposed to rand).
  • davmre 11y ago
    As others have pointed out, the distinction between randn() and rand() is crucial here - the former gives you points from a (spherically symmetric) normal distribution, the latter gives you uniform points from the unit cube. One of the advantages of the normal-sampling route over @ginko's rejection-based method is that in high dimensions almost all of the volume of the unit cube is situated outside of the unit sphere (the volume of the unit cube is always 1, whereas the volume of the unit sphere decreases exponentially with dimension). So the rejection method becomes exponentially slow in high dimensions, while the Gaussian method still works just fine.
    • Retr0spectrum 11y ago
      Surely there would still be a bias towards the corners?
      • darkmighty 11y ago
        Interestingly, there isn't for Gaussians (see my comment below). Other simple distributions tend to give biases towards the corners or axis. Perhaps the Gaussian is unique in this regard? I'm not sure.
        • madcaptenor 11y ago
          In three dimensions (to keep the notation reasonably simple) you need a distribution with density f such that f(x)f(y) f(z) is a function of x^2 + y^2 + z^2 (i. e. the distance from the origin). It looks like the Gaussian (up to scaling) is the only one.
        • OscarCunningham 11y ago
          Yes, the fact that i.i.d. Gaussians give a rotationally symmetric joint distribution is a defining property of the Gaussian.
          • darkmighty 11y ago
            Ah I suspected as much. Can you provide proof? I think I found one, but I'm not sure: Without loss of generality, take f(xi) = k1 * exp(-g(xi)) [1], for some g. Then we need the joint pdf to satisfy f(x1,...,xn) = k2 * exp(-h(R^2)), R=sum(xi^2)^1/2 (the R^2 and h(.) is w.l.g. too). So we get g(x1)+g(x2)=h(x1^2+x2^2). Then assuming the functions g and h analytic we end up needing g(x)= k * x^2, otherwise we get cross terms in the Taylor expansion that can't be cancelled out for all xi. Sounds good? [1] The function f trivially needs to be symmetric, justifying no loss of generality.
  • deleted 11y ago
    [deleted]
• Qantourisc 11y ago
  The easiest way I know how to do it, call rand 3 times: you get a cube -> Calculate the length of the vertex if outside of the sphere reject the sample -> Normalize. IIRC this should be well distributed, yet not fast.
  • ColinWright 11y ago
    And in higher dimensions, extremely slow, as the volume of the sphere fairly rapidly approaches zero, so nearly all your points are outside. You need to work harder than that, as most of the other comments are pointing out.
    • im3w1l 11y ago
      Is the rejection method better at low dimensions? If so, where does the gaussian method catch up?
      • ColinWright 11y ago
        In even dimensions, volume of a sphere of radius 1 is: V_{2k} = pi^k/k! For 10 dimensions, the hypersphere is of volume 2.55, but the enclosing box is of volume 2^10 = 1024. So roughly one in 400 trials gets you a point inside the sphere. Not good odds.
      • Qantourisc 11y ago
        The rejection method is never better I think, work smarter not harder. But it's an easy way to remember and do it. Don't recommend using it when you any actual performance. In terms of CPU power, calculating the length of the vector and doing a branch, I think, is going to be slower then the proper way to do it.
• carlob 11y ago
  Not necessarily the most efficient though: the Box-Muller transform that is probably used in the call to randn is related to the transform you'd need to start from uniformly distributed reals. So your solution ends up pretty much using the same math (while probably calling cos a few more times) and throwing away at least one random number and possibly half of them (depending on the implementation of randn). http://en.wikipedia.org/wiki/Box–Muller_transform http://en.wikipedia.org/wiki/Box–Muller_transform
  • davmre 11y ago
    Yeah, the Gaussian trick is like programming in a high-level language: using existing abstractions lets you write clean and concise code, at the price of some efficiency versus the optimal engineered-from-scratch solution.
    • carlob 11y ago
      I'm not sure I agree, I think the really high-level option is to use a library and not roll your own. For example: https://www.gnu.org/software/gsl/manual/html_node/Spherical-Vector-Distributions.html https://www.gnu.org/software/gsl/manual/html_node/Spherical-...
      • davmre 11y ago
        Fair enough. Though it looks like gsl_ran_dir_nd is implemented using normalized Gaussians, so in this case we're back where we started. :-)
• new1234567 11y ago
  Here's a nice article with pictures on sphere point picking[0]. Here is another [1]. The surface of a sphere in 3-dimensional space has only two dimensions. You can specify any point on the sphere with two coordinates (spherical coordinates). The mathematical expression for this sphere is S-2. If you're using three coordinates to generate the point, you're using more than you need. It's straight algebra to convert from spherical to linear coordinates if you need those for your work. [0] http://mathworld.wolfram.com/SpherePointPicking.html http://mathworld.wolfram.com/SpherePointPicking.html [1] http://stats.stackexchange.com/questions/7977/how-to-generate-uniformly-distributed-points-on-the-surface-of-the-3-d-unit-sphe#comment13055_7984 http://stats.stackexchange.com/questions/7977/how-to-generat...
  • madaxe_again 11y ago
    Yeah, was going to say, this would most readily be addressed by randomly plotting in Spherical Polar Coordinates and then converting.