[abetusk]

Interesting. The assumption is that the Normal Distribution is two dimensional, that it's rotationally invariant and that the X and Y coordinate are statistically independent. From this, they find the polar equation, phi(r) is proportional to a Cartesian one, f(x) * f(y). Since phi(r) = phi(sqrt(x^2 + y^2)) which is then proportional to f(x) * f(y).

More succinctly:

    phi(r) = phi(sqrt(x^2 + y^2)) = f(x) * f(y)
    phi(sqrt(x^2 + 0)) = f(x) * f(0) = lamba * f(x)
    -> phi(x) = lamba * f(x)
    -> phi(r) = lamba f(sqrt(x^2 + y^2)) = f(x) * f(y)

With the last line because of rotational invariance and statistical independence of the two dimensional axies. I haven't followed the rest but I assume (maybe with some other minor assumptions?) that the equation, lambda f(sqrt(x^2 + y^2)) = f(x) * f(y), uniquely determines the Gaussian.

I've long struggled to find a clear explanation of where the Normal formula comes from. The best description I've seen is by deriving limiting distribution of sums of uniform distributions on a unit interval, say. The sum of identical and independently distributed random variables is a convolution which can be 'de-convolved' by taking the Fourier Transform. After the Fourier transform, the sum of the random variables turns into a product which can easily be approximated. There's extra work involved in proving the Fourier transform of a Gaussian is itself Gaussian and some other technicalities (not to mention this is only for uniform distributions that are identical) but this seems much more motivated to me than any other descriptions I've heard, including this one.

As a benefit, if I remember correctly, the same trick works to derive the basics of Levy stable distributions as well.