I'm trying to compute the probability of the distances between points on a 2D domain that have complete spatial randomness (CSR). From this wikipedia page on CSR, the probability of locating the $N$th neighbour of a point within a radius $r$ is
$$ P(N, r) = \frac{D}{(N - 1)!} \lambda^N r^{DN - 1} e^{-\lambda r^D}, $$
where
$$ \lambda = \frac{\rho \pi^{D/2}}{\Gamma(D/2 + 1)}. $$
Here, $\rho$ is the event density and $D ( = 2)$ is the number of dimensions. However, what I ultimately want to calculate is, I guess, the integral of $P(N, r)$ wrt N, i.e.
$$ P(r) = \int_0^\infty \frac{D}{(N - 1)!} \lambda^N r^{DN - 1} e^{-\lambda r^D} dN = \frac{D}{r} e^{-\lambda r^D} \int_0^\infty \frac{\lambda^N r^{DN}}{(N - 1)!} dN. $$
Trouble is I don't know how to evaluate the integral. I've tried using Sage's symbolic integration but to no avail. Can anyone help?
Edit 1:
In light of joriki's comments, it might be useful if I added some more background to my problem. What I have is $n$ points which are supposedly randomly distributed in the $(x,y)$ plane, as is shown in 1 below:
I calculate the distribution of these points by 'binning' the distances between all pairs of points. Finally, I want to compare this distribution (2) to some theoretical distribution (hence my question on this forum).
Edit 2:
Following Ritz' comment, I tried taking the sum over $N$, such that
$$ P(r) = \frac{D}{r} \lambda r^D e^{-\lambda r^D} \sum_1^\infty \frac{\left(\lambda r^D\right)^N}{(N - 1)!}, $$
Or, using a change of variable, $x = \lambda r^D$,
$$ P(x) = \frac{D}{(x/\lambda)^{1/D}} x e^{-x} \sum_1^\infty \frac{x^N}{(N - 1)!} = \frac{D}{(x/\lambda)^{1/D}} x e^{-x} \sum_0^\infty \frac{x^N}{N!} = D x^{1 - 1/D} \lambda^{1/D}. $$
Clearly, this new function isn't going to describe the data in (2). What (to me) is curious is that the $P(N,r)$ function, if scaled and stretched, describes the data very well (imaging shifting the red curve over a bit and lowering its height). Is this because I've misunderstood the meaning of this distribution? Should it perhaps be the probability of finding $N$ neighbours within a radius $r$ and, actually, I'm only trying to find one neighbour at a time?
It seems like you're actually just calculating the probability distribution of $d(x,y)$, where $x$ and $y$ are placed at random within a rectangle (square?). Note that the pdf of the distance from a point to its $k$-th nearest neighbor doesn't depend on the size of the square, in the limit that the square becomes large, but what you're calculating certainly does. The result (for a square of side $1$) is known to be $$ f(d)=\begin{cases}2d(\pi+d^2-4d) && 0<d\le 1 \\ 2d\left(-2+4\sin^{-1}\left(\frac{1}{d}\right)+4\sqrt{d^2-1}-\pi+d^2\right) && 1<d\le \sqrt{2}, \end{cases} $$ if I haven't made any typos. (For instance, see J. Philip's paper here.)