# 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 & ] Commented Mar 31, 2021 at 6:08
• And you can just pick $n$ as a Poisson?
– apg
Commented Mar 31, 2021 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. Commented Mar 31, 2021 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]]


– apg
Commented Mar 30, 2021 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.
– apg
Commented Apr 1, 2021 at 13:36
• Thanks, I corrected this. Commented Apr 1, 2021 at 14:14