# Integral of Gaussian ring/shell

I would like to know the integral of the function

$$f(\mathbf{x}) = {1 \over \sqrt{2\pi \sigma^2}} \exp\left\{- {(|\mathbf{x}| - \mu)^2 \over 2\cdot \sigma^2}\right\}$$

over an $n$-dimensional vector space. It is a ring of radius $\mu$ around the origin, where the deviation from the radius yields a fall-off according to a Gaussian with standard deviation $\sigma$.

I had no success with WolframAlpha/Maxima/SymPy, probably due to the norm $|\mathbf{x}|$.

Is the right approach to switch to spherical coordinates, and then integrate along the ring first, and finally over the radius?

$$g(R, \phi) = {1 \over \sqrt{2\pi \sigma^2}} \exp\left\{- {(R - \mu)^2 \over 2\cdot \sigma^2}\right\}$$ so that $$\int_0^\infty dR \int_0^{2\pi} d\phi ~~ g(R, \phi) = 2\pi\cdot 1 = 6.28$$

I do not think this result is correct, for 2d it should be 0.087 according to a paper I have (for 5d 0.0017; using $r=2$, $w=0.1$).

## Update

I missed out to use the volume element dV. In the general integration over the n-sphere (from here):

$$d^nV = r^{n-1}\sin^{n-2}(\phi_1)\sin^{n-3}(\phi_2)\cdots \sin(\phi_{n-2})\, dr\,d\phi_1 \, d\phi_2\cdots d\phi_{n-1}$$

Now the integral $$\int g(R) R^{n-1} dR$$ is the (n-1)th moment of the gaussian, $M_{n-1} = \operatorname{E}\left[X^{n-1}\right]$ (see here, and scipy.stats.norm.moment(n-1, loc=r, scale=sigma))

So then we are only missing the rest of the volume elements. This is the "surface area" of a unit n-sphere (from here and here). $$S_{n}={(n+1) V_{n+1}}$$ with $$V_n = \frac{\pi^{n/2}}{\Gamma(\frac{n}{2} + 1)}$$

So I should be able to write the complete integral as

$$I = \int_{\cal{R}^n} g(R, \mathbf{\theta}) dV = S_n \cdot M_{n-1}$$

For n=2: The first moment with $\mu=2$, $\sigma=0.1$ is $M_1=2$. The surface $S_2=4\pi$. This gives $I_2 = 25.1$, which is still not what it ought to be.

• In 3D, what are you integrating over? A sphere? Are we essentially integrating over an $n$-ball? – Chris K Dec 8 '13 at 18:19
• Integrating over $\cal{R}^n$. A $n$-ball is fine. – j13r Dec 8 '13 at 18:25
• As for your calculation you are missing the Jacobian in the 2D calculation. – Chris K Dec 8 '13 at 18:51
• @ChrisK, could you have another look now – j13r Dec 8 '13 at 19:07
• The update seems to assume that $\mu=0$, is this what you want? – Did Dec 8 '13 at 19:16

Let $\varphi$ and $\Phi$ denote the standard normal PDF and CDF, that is, $$\varphi(t)=\frac1{\sqrt{2\pi}}\mathrm e^{-t^2/2},\qquad\Phi(t)=\int_{-\infty}^t\varphi=\int_{-t}^\infty\varphi.$$ Using spherical coordinates as you explained one is led to compute $\displaystyle\int g(R)R^{n-1}\mathrm dR$ the integral in your post, which is $$A_n(\mu,\sigma^2)=\int_0^\infty r^{n-1}\varphi\left(\frac{r-\mu}\sigma\right)\mathrm dr/\sigma.$$ The change of variable $r=\sigma t$ yields $$A_n(\mu,\sigma^2)=\sigma^{n-1}B_{n-1}(\mu/\sigma),\qquad B_k(u)=\int_0^\infty t^k\varphi(t-u)\mathrm dt.$$ Now, for $k=0$, $$B_0(u)=\int_{-u}^\infty \varphi=\Phi(u).$$ Using the fact that $\varphi'(t)=-t\varphi(t)$, one sees that $$B_k'(u)=\int_0^\infty t^k(t-u)\varphi(t-u)\mathrm dt=B_{k+1}(u)-uB_k(u).$$ One deduces from this that there exists some sequences $(P_k)$ and $(Q_k)$ of polynomials such that, for every $k\geqslant0$, $$B_k=P_k\varphi+Q_k\Phi.$$ These are uniquely determined by the initial condition $(P_0,Q_0)=(0,1)$ and by the recursion $$P_{k+1}(u)=P'_k(u)+Q_k(u),\qquad Q_{k+1}(u)=Q'_k(u)+uQ_k(u).$$ For example, $$(P_1(u),Q_1(u))=(1,u),\quad (P_2(u),Q_2(u))=(u,1+u^2),\quad (P_3(u),Q_3(u))=(2+u^2,3u+u^3).$$ One recognizes (or a direct computation shows that), for every $k$, $$Q_k(u)=\mathrm i^{-k}H_k(\mathrm iu),$$ where $H_k$ is the $k$th Hermite polynomial.
• I also computed the momentum numerically, and it is indeed $A_2=M_1=2$. But if I compute the integral over the Gaussian shell numerically, I get, like in the paper, 0.00691. Something is missing (sorry for the edit-confusion). – j13r Dec 8 '13 at 22:39