# Numerical Integration of a Gaussian Distribution in Polar Coordinates

I want to numerically evaluate a 2D-integral of a specific probability distribution over some given area (I use MATLAB so all the code below is MATLAB code). I broke down the problem so that it appears in the following integration: $$\int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} \frac{1}{2 \pi \sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}} dx dy = 1$$ Doing this numerically in MATLAB with cartesian coordinates looks something like this:

Dx = 0.05;
Dy = 0.05;

x = -1:Dx:1;
y = -1:Dy:1;

[XX YY] = meshgrid(x,y);
pdf1 = 1/(2*pi*sigma^2) * exp(- (XX.^2 + YY.^2)/(2*sigma^2));

area1 = sum(sum(pdf1)) * Dx * Dy % numerical integration
error1 = abs(1-area1) % almost no error: 1.1102e-15


However, when I try the same thing in a straightforward manner with polar coordinates, I get the following:

rsteps = 50;
tsteps = 20;

Dr = 1/rsteps;
Dt = 2*pi/tsteps;

r = [0:1:rsteps]*Dr;
t = [0:1:tsteps-1]*Dt;

[RR TT] = meshgrid(r,t);
pdf2 = RR/(2*pi*sigma^2) .* exp(- (RR.^2)/(2*sigma^2));

area2 = sum(sum(pdf2)) * Dr * Dt % area in radius-phase-plane
error2 = abs(1-area2) % substantial part is missing: 0.0033


At first I thought that I did some error with the conversion to polar coordinates, but letting the value rsteps go to very large number, the error becomes almost zero, but the convergence is pretty bad (e.g., for 50000 steps in $r$-direction, the error is of the order $10^{-9}$ which is way worse than in cartesian coordinates with a much cruder grid).

My intuition says that my straightforward numerical integration fails because I cannot properly evaluate the distribution at the origin ($x=0$, $y=0$) with the polar coordinates, so there is an area element missing (something that would look like a circle around the origin in cartesian coordinates). Is that correct or is there another problem? How can I fix this with a reasonable grid size anyway. Ideally, I would like to use the same grid, and achieve the same error as in cartesian coordinates.

-

This is to be expected; in fact we can calculate from the error you report that your $\sigma$ value must have been around $0.1$.

From the way you set up your $r$ grid, you're effectively using the trapezoidal rule. You've got $N+1$ points and you're normalizing by $N$. The two outer points should be weighted by $1/2$, but since the values there are $0$ that doesn't matter. The error of the trapezoidal rule can be estimated by expanding the function at the centre of the trapezoid and integrating the missing quadratic term, which yields

$$\int_{-h/2}^{h/2}\frac{f''(x_0)}2(x-x_0)^2\mathrm dx=\frac{f''(x_0)}{12}h^3\;,$$

where $h$ is the interval length. Approximating the sum over all intervals by an integral then yields

$$\int_a^b\frac{f''(x_0)}{12}h^3\frac1h\mathrm dx=\frac{h^2}{12}\left(f'(b)-f'(a)\right)=\frac{(b-a)^2}{12N^2}\left(f'(b)-f'(a)\right)\;,$$

where $N$ is the number of intervals. In your case $b-a=1$, $N=50$, $f'(b)\approx0$ and $f'(a)=1/\sigma^2$ (after cancelling $2\pi$ with the angular integration), so the error would be expected to be about $-(12\cdot50^2\sigma^2)^{-1}$, which is $-0.0033$ for $\sigma\approx0.1$.

The reason you didn't have this problem in the first integration is that you had $f'(a)\approx f'(b)\approx0$ in that case.

I guess the moral of the story is to use quadrature rules with higher error orders.

-
 Thanks a lot for your answer! Some comments: I thought I was using a simple rectangle-rule, but ok, its equivalent to the trapezoidal rule in my case. So for me the moral of the story is that with this rule one has to be careful if the the derivative of the function at the boundaries of the integration is not zero, otherwise the error term will be high. – Christian Aug 10 '12 at 12:00 @Christian: If you interpret this as using the rectangle rule, then you've got the wrong integration region: For the rectangle rule to be of third order like the trapezoidal rule, the values have to be taken in the middle of the interval, not at the lower end, so you'd be integrating from $-Dr/2$ to $1+Dr/2$. It so happens that it works out as the trapezoidal rule instead because the boundary values are both $0$. About the moral: Yes, a third-order error is quite high, and it's easy and computationally inexpensive to use e.g. Gaussian quadrature that can have any error order you want. – joriki Aug 10 '12 at 12:18