From: Eugene Salamin <gene_salamin@yahoo.com> To: math-fun <math-fun@mailman.xmission.com> Sent: Wed, March 23, 2011 6:03:33 PM Subject: [math-fun] Integral of Gaussian over sphere Can anyone find a closed form expression for the integral of exp[(-1/2) (a[1] x[1]^2 + a[2] x[2]^2 + ... + a[n] x[n]^2)] over the unit (n-1)-sphere? In the n=2 case, we have a x^2 + b y^2 = a cos(t)^2 + b sin(t)^2 = (1/2) ((a+b) + (a-b) cos(2t)), and integrating over t in [0, 2 pi] gives 2 pi exp((-1/4) (a+b)) BesselI(0, (1/4) (a-b)). -- Gene ______________________________________________ I can express this as a single integral along a contour. I will slightly change the problem. Let I be the integral of exp(-(a[1] x[1]^2 + a[2] x[2]^2 + ... + a[n] x[n]^2)) over the (n-1)sphere of radius r. Since the integration is over a compact region, the integral is well defined for all complex a[k]. I will assume a[k] > 0, allowing analytic continuation to do the rest. Let ρ be the radial coordinate. Introducing the delta function δ(ρ-r) into the integrand, we can then integrate over all n-space. We have δ(ρ^2-r^2) = δ(ρ-r)/(2r) + δ(ρ+r)/(2r) = δ(ρ-r)/(2r), since ρ+r > 0. It is well-known, at least to physicists, that one can represent the delta function as δ(t) = (1/(2π)) int(exp(iωt), ω=-inf..inf). Then δ(ρ-r) = 2r δ(ρ^2-r^2) = 2r δ(sum(x[k]^2, k=1..n)-r^2). Then I is r/π times the integral of exp(-(sum((a[k]-iω) x[k]^2, k=1..n)) exp(-iωr^2) over all x[1], ..., x[n], ω from -infinity to +infinity. First integrate over the x[k] the first exponential expression above. This step requires positive a[k], or at least Re a[k] > 0. This integral is sqrt(π^n / prod((a[k]-iω), k=1..n)), and we have I = (r/π) π^(n/2) int(exp(-iωr^2)/sqrt(prod((a[k]-iω), k=1..n)), ω=inf..inf). The integrand has poles or branch points at ω = -i a[k], which lie on the negative imaginary axis. Also, exp(-iωr^2) is analytic in the lower ω half-plane, and --> 0 for Im ω --> -inf. Hence the contour can be deformed to go up from -i inf on the left side of the negative imaginary axis, clockwise around the -i a[k] closest to the origin, and down to -i inf on the right side of the imaginary axis. If two of the a[k] are equal, they contribute a pole to the integrand, and their contribution to I can be evaluated by residues. A pair of unequal a[k] can be joined by a branch cut, with the contour going around the cut. If n is odd, at least one branch cut necessarily extends to (-i inf). Note that poles are contoured clockwise, so they contribute -2πi times the residue. When n=2, a[1]=a[2]=a, we get I = 2πr exp(-a r^2), as expected. But for n>2, it looks wrong; r should appear to the n-1 power. But consider n=4, a[1]=a[2]=a, a[3]=a[4]=b. Working out the residues, we get I = 2π^2 r (exp(-ar^2) - exp(-br^2))/(b-a). As b-->a, this becomes 2π^2 r^3 exp(-ar^2), as it should. The double pole at ω=-ia differentiates the exp(-iωr^2), providing the extra two powers of r. More generally, suppose we scale r-->sr, a[k]-->a[k]/s^2, ω-->ω/s^2, dω-->dω/s^2. Then exp(-iωr^2) is unchanged, 1/sqrt(prod((a[k]-iω), k=1..n)) is multiplied by s^n, the integral is multiplied by s^(n-2), and finally I is multiplied by s^(n-1), displaying the r^(n-1) dependence. When n=2, a[1] /= a[2], there is a pair of branch points joined by a cut. Its Fourier transform integrated along the cut is a Bessel function, but I'm too lazy to work it out to verify its correctness. When n=4 with no two a[k] equal, the formula asks us to take the Fourier transform along the cuts of the function, which if integrated just by itself, gives an elliptic integral. This, and for larger n, quickly exceeds my skill with special functions. Perhaps the reader is curious as to my interest in this problem. I've been working on the problem of fitting a line y = ax + b to data (x[i], y[i]) when both x[i] and y[i] are subject to error. The usual least squares solution assumes that the x[i] are free of error. A proper Bayesian solution must provide, not just a single "best fit", but rather a probability distribution p(a,b|x[i],y[i],σx,σy). The proposed problem, with n=2, arises in normalizing this probability to unit total integral. So this can be done with nothing worse than BesselI. More generally, fitting a plane y = sum(a[k] x[k], k=1..n) + b yields a probability distribution p(a[k],b) whose normalization requires the integral of an (n+1)-dimensional Gaussian over the n-sphere. -- Gene