9
$\begingroup$

Given a positive integer n and a small positive real number s.

How to create distribution for n real numbers p1, p2, p3, ..., pn such that

0 ≤ p1 < p2 < p3 <... < pn ≤ 1
and
p1, p2-p1, p3-p2, ... 1-pn ≥ s

?

$\endgroup$
4
  • $\begingroup$ Do you want a symbolic distribution or a way to generate samples that satisfy your restrictions? $\endgroup$
    – JimB
    Mar 7 at 3:41
  • 1
    $\begingroup$ I's not exactly a solution to your problem, but related is HardcorePointProcess like this: points = Sort@Flatten@RandomPointConfiguration[HardcorePointProcess[50.0, 0.01, 1], Line[{{0}, {1}}]]["Points"]; . I don't think this has the boundary contraints you want though, and it also doesn't let you specify an exact number of points. $\endgroup$
    – flinty
    Mar 7 at 10:50
  • 4
    $\begingroup$ Hello. I am a lost mathematician who has wandered in from MSE. I don't know anything about Mathematica! However, if I understand your question correctly - "how can I uniformly generate real numbers $p_i$ according to these constraints" - I have something to say. This problem is amenable to a nice elementary trick: consider instead the numbers $q_i = p_i - is$ (and note that this transformation is measure-preserving). The corresponding constraint is simply $0 \le q_1 \le q_2 \le \dots \le q_n \le 1 - (n + 1)s$... $\endgroup$ Mar 8 at 2:00
  • 2
    $\begingroup$ So you can solve this problem without rejection sampling by using your favourite method to uniformly generate sorted numbers in $[0, 1 - (n + 1)s]$ (eg stackoverflow.com/a/66366347 or just generate random numbers & sort them), and then shifting them to get $p_i$. $\endgroup$ Mar 8 at 2:00

5 Answers 5

9
$\begingroup$

If you want to generate independent random samples under the stated restrictions and $s<1/(n+1)$, then consider the following:

n = 10;
s = 1/20;
p = ConstantArray[0, n];
count = 0;
SeedRandom[12345];
While[Min[p] < s || Max[p] > 1 - s || Min[Differences[p]] < s,
 count = count + 1;
 p = Sort[RandomVariate[UniformDistribution[{s, 1 - s}], n]]]

p
(* {0.16256, 0.296157, 0.368503, 0.446481, 0.510143, 0.56314, 0.675157,  0.726587, 0.836507, 0.93472} *)
s // N
(* 0.05 *)
Min[Differences[p]]
(* 0.0514309 *)
count
(* 359 - Number of samples that had to be rejected *)
$\endgroup$
4
  • $\begingroup$ neat, +1 of course :) $\endgroup$
    – ubpdqn
    Mar 7 at 11:14
  • $\begingroup$ Hello Jim. This might be a really silly question, so sorry about that - but I am trying to understand how your code works and I can't work out the purpose of the part Min[p] < s || Max[p] > 1 - s - isn't that (almost always*) guaranteed by sampling from a uniform distribution on $[s, 1 - s]$? Does it have to do with the starting condition of the loop, maybe? *: It also seems that the question only asks for a non-strict inequality between $p_1$ and $s$ and $p_n$ and $1 - s$ anyway. $\endgroup$ Mar 8 at 2:07
  • 1
    $\begingroup$ Not a silly question (and thank you for the diplomatic approach). I needed to include that restriction because I started out with p = ConstantArray[0, n]. But had I started out with p=RandomVariate[UniformDistribution[{s,1-s}], n, I would be able to remove those pieces of code. But certainly you are correct in that given the way I did it, Max[p] > 1 - s isn't necessary at all. $\endgroup$
    – JimB
    Mar 8 at 3:36
  • $\begingroup$ Also, note that this approach is extremely inefficient when $s$ is close to $1/(n+1)$. $\endgroup$
    – JimB
    Mar 8 at 6:18
8
$\begingroup$

Thanks to the comment from Izaak van Dongen, it can be solved with straightly :

umd[n_, s_] := Sort[RandomReal[{0, 1 - (n + 1) s}, n]] + s  Range[n];

Pick 0 < p1 < p2 < p3 < p4 < p5 < 1, such that difference of adjacent elements is at least 0.1

Table[ListPlot[umd[5, 1/10], ImageSize -> 50, AspectRatio -> 2,
   PlotRangePadding -> {0.5}, PlotStyle -> Red, Filling -> Axis,
   FillingStyle -> Red, PlotRange -> {{1, 5}, {0, 1}}], 10] // Row

enter image description here

I wonder if this is better than the previous answers.

$\endgroup$
2
  • $\begingroup$ This is far faster than any of the others. Also, the HardcorePointProcess seems to involve a "mystery parameter" that can affect the distribution of any of the ordered points. In other words, there are speed and potential distribution differences among the answers so "better" depends on what you're looking for. $\endgroup$
    – JimB
    Mar 8 at 21:23
  • $\begingroup$ You probably want some error handling and to raise a message if n*s > 1, because in something like umd[11, .1] the first and last elements are out of bounds. $\endgroup$
    – flinty
    Mar 20 at 13:01
6
$\begingroup$

Thanks for @Ulrich Neumann's advice for code improvement.


The idea of the below code is to iteratively generate points, ensuring that each new point maintains the minimum distance $s$ from all previously generated points.

(*Initialize variables*)n = 10;


(*Number of points*)s = 1/20;

(*Minimum distance*)
points = {};(*List to hold the points*)

(*Function to check if a new point is at least distance s from all existingpoints*)
isValidPoint[newPoint_, points_, s_] := 
 AllTrue[points, RealAbs[# - newPoint] >= s &];

(*Generate points*)

While[Length[points] < n, newPoint = RandomReal[{s, 1 - s}];
  If[isValidPoint[newPoint, points, s], AppendTo[points, newPoint];];
  (*Optional:
  include a check to break out if it's not possible to add more points*)];

(*Sort points for better visualization*)
points = Sort[points];


points
Min[Differences[points]] >= s  (*True*)
$\endgroup$
1
  • $\begingroup$ @Aspen Nice solution! Change Abs to RealAbs (and newPoint = RandomReal[{s, 1-s}]) which makes code run considerable faster for greater n. $\endgroup$ Mar 7 at 8:17
3
$\begingroup$

The most "don't make me think" approach in order to avoid scratching head regarding uniformity of sampling of such a distribution is to specify it as an $n$-dimensional implicit geometric region and sample points from it using RandomPoint. The problem with this approach is that RandomPoint performs some sort of an algebraic decomposition of the system of inequalities and the run time of this is dependent on $n$ and $s$, in practice increasing exponentially as $n$ increases, putting limit of practicality somewhere at $n \ll 15$. The good part is that once this decomposition is done, sampling many values is not horribly slow.

In numerical sampling (and when $s \gt 0$) you are never going to see $0$ or $1$ in results anyway, so I simplified $\le$ to $\lt$ in these edge cases.

samples =
  With[{n = 8, s = 1/10},
   With[{vals = {0, Sequence @@ Array[p, n], 1}},
    RandomPoint[
     ImplicitRegion[
      Less @@ vals &&
       And @@ (# >= s & /@ Differences[vals]),
      Evaluate[Array[p, n]]],
     25]]];

We can see that values are distributed "evenly":

NumberLinePlot[samples]

enter image description here

Also, differences between consecutive values are at least $\frac{1}{10}$:

NumberLinePlot[Differences /@ samples]

enter image description here

With $25000$ samples we can also see distributions of each $p_n$:

Histogram[Transpose[samples], {0.005}, "PDF"]

enter image description here


Similarly, symbolic CDFs for individual $p_n$ can be acquired using Probability and Conditioned:

cdfs =
  With[{n = 8, s = 1/10},
   With[{vals = {0, Sequence @@ Array[p, n], 1}},
    ParallelTable[
     Probability[
      Conditioned[
       p[i] <= x,
       Less @@ vals &&
        And @@ (# >= s & /@ Differences[vals])],
      Distributed[Array[p, n], UniformDistribution[Table[{0, 1}, n]]]],
     {i, n}]]];

Symbolic PDFs correspond to those seen on the histogram:

Plot[D[cdfs, x], {x, 1/10, 9/10}, Evaluated -> True]

enter image description here

$\endgroup$
2
$\begingroup$

You can do this by simulating a HardcorePointProcess till you've hit the right number of points:

n = 10; (* number of points *)
r = 0.05; (* exclusion radius *)
reg = Interval[{0, 1}];
mu =  30 (* You'll need to guess an appropriate value for this *)
proc = HardcorePointProcess[mu, r, 1];
ResourceFunction["Until"][
  points["PointCount"] === 10
  ,
  points = RandomPointConfiguration[proc, reg]
];
Sort @ Flatten @ points["Points"]
Differences[%] (* check the distances *)

You can make a rough estimate for mu as:

volume = RegionMeasure[reg];
mu =  n/(volume + r - n * r)

You can find this by estimating by assuming that the number of points is roughly equal to the intensity of the process times the effective volume of the region (i.e., the volume without the excluded hard cores):

Clear[mu, volume, n, r]
Solve[mu * (volume - (n - 1)  r) == n, mu]

If the formula give a negative number, you're probably trying to stuff too many points into the interval to fit, but you can still try a really high value of mu to see if you can get a result.

$\endgroup$
2
  • $\begingroup$ When I run this with SeedRandom[12345], I get p[1] < 0.05 and p[n] > 1-0.05 although all of the distances between points are at least 0.05. Maybe the OP doesn't really need the restrictions at the extremes. $\endgroup$
    – JimB
    Mar 8 at 4:37
  • $\begingroup$ Oh, I missed the boundary constraints. I guess you can fit those in by just shrinking the region a bit. $\endgroup$ Mar 8 at 8:31

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.