4
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ Isn't this just sampling from a truncated 2D normal distribution? For n points, radius R, Select[ RandomVariate[NormalDistribution[], {n, 2}], Norm[#] <= R & ] $\endgroup$
    – Szabolcs
    Commented Mar 31, 2021 at 6:08
  • $\begingroup$ And you can just pick $n$ as a Poisson? $\endgroup$
    – apg
    Commented Mar 31, 2021 at 14:43
  • 1
    $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Commented Mar 31, 2021 at 16:52

2 Answers 2

4
$\begingroup$

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]]

point process

$\endgroup$
1
  • $\begingroup$ Ok I will download 12.2! $\endgroup$
    – apg
    Commented Mar 30, 2021 at 15:22
3
$\begingroup$

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]}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thank you, that works well, particularly if you just want the points. Only thing is p=ProbabilityDistribution should be should be p1=ProbabilityDistribution. $\endgroup$
    – apg
    Commented Apr 1, 2021 at 13:36
  • 1
    $\begingroup$ Thanks, I corrected this. $\endgroup$ Commented Apr 1, 2021 at 14:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.