# Sampling inhomogeneous spatial point processes

I am looking to generate a point process (random set of points) in a disk of radius $$R$$, but where the intensity measure decays, for $$r$$ the distance from the origin of the disk, as $$\exp(-r^2/2)$$. Is there a standard way to do this in Mathematica?

In the simpler case where the intensity measure is uniform I can sample by using a set of random reals, where the number of points is Poisson. To model the non-uniform intensity, I could just discretise the space, and pick a random number of points in each cell which is decaying as appropriate, but is there a more `built in' way?

• Isn't this just sampling from a truncated 2D normal distribution? For n points, radius R, Select[ RandomVariate[NormalDistribution[], {n, 2}], Norm[#] <= R & ] – Szabolcs Mar 31 at 6:08
• And you can just pick $n$ as a Poisson? – apkg Mar 31 at 14:43
• Yes, exactly. Of course, sampling the x and y coordinates independently works only for the normal distribution. This simplified method won't work for most other distributions. – Szabolcs Mar 31 at 16:52

Assuming you're using the new spatial statistics from v12.2:

{x0, y0} = {1.5, 2};
σ = .7;
R = 3.0;
reg = Disk[{x0, y0}, R];
pts = RandomPointConfiguration[
InhomogeneousPoissonPointProcess[
Function[10000*Exp[-((#1 - x0)^2 + (#2 - y0)^2)/(2*σ)]],
2], reg];
Show[RegionPlot[reg], ListPlot[pts]]

• Ok I will download 12.2! – apkg Mar 30 at 15:22

Here is a pedestrian way to achieve the same that also works with older MMA:

We first define the normalized ProbabilityDistribution of the radius:p

r0 = 1;
pdf[x_] = Exp[-x^2/2]/Integrate[  Exp[-y^2/2], {y, 0, r0}]
p = ProbabilityDistribution[pdf[x], {x, 0, r0}];

With this distribution we can create r-values: r. The angle is uniformly distributed, what we can achieve simply by RandomReal. From r and phi we can get x/y coordinates and draw the points:

n = 10^4;
pts = MapThread[#1 {Cos[#2], Sin[#2]} &, {RandomVariate[p, n],
RandomReal[{0, 2 Pi}, n]}];
Graphics[{Circle[{0, 0}, r0], PointSize[0.001], Point[pts]}]

• Thank you, that works well, particularly if you just want the points. Only thing is p=ProbabilityDistribution should be should be p1=ProbabilityDistribution. – apkg Apr 1 at 13:36
• Thanks, I corrected this. – Daniel Huber Apr 1 at 14:14