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.

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.

Surely there would still be a bias towards the corners?

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.

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.

Yes, the fact that i.i.d. Gaussians give a rotationally symmetric joint distribution is a defining property of the Gaussian.

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.

Perhaps coincidentally, someone just asked about this at math.stackexchange: http://math.stackexchange.com/questions/1255637/joint-pdf-of...

Yes, that was me :)

I wanted to make sure I got a proof, since I didn't really find this elsewhere.

You have to assume that the one-dimensional marginals are Gaussian for otherwise the statement is not correct.

Hmm I'm not sure I get your point. I'm trying to prove that the joint pdf of N iid RV's is isotropic if and only if the RV's are gaussian. If I assume the pdfs are gaussian in the first place the proof isn't valid?

Oh, I see. I thought you were trying to prove something different.

No, OP is correct; he spoke about an isotropic normal distribution (randn as opposed to rand).

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.

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.

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

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

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.

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

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. :-)

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 [1] http://stats.stackexchange.com/questions/7977/how-to-generat...

Yeah, was going to say, this would most readily be addressed by randomly plotting in Spherical Polar Coordinates and then converting.

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.

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.

Is the rejection method better at low dimensions? If so, where does the gaussian method catch up?

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.

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.