# Sampling from density specified in terms of spacing ratios?

Given a sorted vector of positive reals $$\{\lambda_1,\lambda_2,\ldots,\lambda_P\}$$ and corresponding vector of spacing ratios $$r_i=\frac{\lambda_i-\lambda_{i-1}}{\lambda_{i-1}-\lambda_{i-2}}$$, I know that as $$P\to \infty$$, the density of $$r$$ approaches the following

$$P(r)=\frac{27 \left(r^2+r\right)}{8 \left(r^2+r+1\right)^{5/2}}$$

How would I get a sample of $$\lambda$$'s?

density = 27 (r + r^2)/(8 (1 + r + r^2)^(5/2));


This comes from formula 7 in this paper

• So, your question is really just how to get a sorted list of reals? Commented Jul 8, 2023 at 10:04
• Sort[RandomReal[{0, 10}, 20]] gives you a sorted list of 20 reals randomly selected from the range {0,10}. Commented Jul 8, 2023 at 10:05
• If you know $r_i$, you can determine $\lambda_i$. Just notice that the map is not unique. The ratio is translationally invariant, therefore without loosing generality we can set $\lambda_1=0$. Furthermore, the ratio is scale-invariant. Therefore, we can set $\lambda_2=1$ without loosing generality. Thus, generate random $r_i$ according to your distribution and perform the transform. Commented Jul 8, 2023 at 11:25
• Or is your goal to find the distribution analytically? Commented Jul 8, 2023 at 11:26
• Is there reason to believe that the $R_i$’s are not serially correlated?
– JimB
Commented Jul 8, 2023 at 18:27

This follows @yarchik 's comment.

First, generate a random sample of values from the random variable $$R$$ with density

density = 27 (r + r^2)/(8 (1 + r + r^2)^(5/2));


Using the inverse transform sampling method find the cdf of $$R$$ and set that to a value $$U$$ from a standard uniform distribution:

cdf = (Integrate[density, {r, 0, r0}, Assumptions -> r0 > 0]) /. r0 -> r
(* (-2 - 3 r + 3 r^2 + 2 r^3 + 2 (1 + r + r^2)^(3/2))/(4 (1 + r + r^2)^(3/2)) *)
sol = r /. Solve[cdf == u, r, Reals]


It turns out that the 5-th and 6-th roots are what we want when $$0 and $$1/2, respectively. Writing a piecewise function for obtaining a random sample from $$R$$ follows:

randomR[uu_] := Piecewise[{{sol[[5, 1]] /. u -> uu, 0 < uu < 1/2},
{sol[[6, 1]] /. u -> uu, 1/2 < uu < 1}, {1, uu == 1/2}, {0, uu == 0}}, ∞]


As a partial test of whether this works, take random sample and compare the associated histogram with the density:

n = 10000;
SeedRandom[12345];
rr = randomR[#] & /@ RandomReal[{0, 1}, n];
Show[Histogram[rr, "FreedmanDiaconis", "PDF"],
Plot[density, {r, 0, 10}]]


Update 1: @kglr in a comment below showed that a random sample can be used with the simpler Mathematica code:

pd = ProbabilityDistribution[density, {r, 0, Infinity}];
SeedRandom[12345];
rr = RandomVariate[pd, n]


Now generate values of $$\lambda$$. From the relationship given we have:

Solve[ri == (λi - λi1)/(λi1 - λi2), λi][[1, 1, 2]]
(* λi1 + ri λi1 - ri λi2 *)


$$\lambda_i=\lambda_{i-1}+r_i (\lambda_{i-1}- \lambda_{i-2})$$

Putting this all together:

(* Random sample of R values *)
n = 1000;
SeedRandom[12345];
rr = randomR[#] & /@ RandomReal[{0, 1}, n];

(* Generate λ values *)
λλ = ConstantArray[0, n];
λλ[[1]] = 0;
λλ[[2]] = 1;
Do[λλ[[i]] = λλ[[i - 1]] + rr[[i]](λλ[[i - 1]] - λλ[[i - 2]]), {i, 3, n}]

(* Plot results *)
ListPlot[λλ, Frame -> True, FrameLabel -> (Style[#, Italic, 18] &) /@ {"i", "λ[i]"}]


I don't know the subject matter so I have no idea if the resulting figure makes any sense or if the cited paper makes any sense or if I've made an error. All are possible.

Update 2: Finding a nice closed-form formula for the density of $$\lambda_i$$ does not seem likely as the random variables $$\lambda$$ become more and more difficult to deal with. We have for the first few $$\lambda$$ random variables the following pattern:

$$\lambda_3=1+R_3$$ $$\lambda_4=1+R_3+R_3 R_4$$ $$\lambda_5=1+R_3+R_3 R_4+R_3 R_4 R_5$$ $$\lambda_6=1+R_3+R_3 R_4+R_3 R_4 R_5+R_3 R_4 R_5 R_6$$

@kglr also showed how to produce this using the following:

solλ = RSolve[{λ[m] == λ[m - 1] + r[m] (λ[m - 1] - λ[m - 2]), λ[1] == 0, λ[2] == 1},
λ[m], m][[1]] /. {K[1] -> j, K[2] -> i}


The above can be simplified for each value of $$m$$:

solλ /. m -> 3 // Expand
(* {λ[3] -> 1 + r[3]} *)

solλ /. m -> 6 // Expand
(* {λ[6] -> 1 + r[3] + r[3] r[4] + r[3] r[4] r[5] + r[3] r[4] r[5] r[6]}


However, one can produce the density for $$\lambda_4$$ numerically (as opposed to random sampling) by integrating the product of the densities times the Jacobian with r4 replaced by (-1 - r3 + \[Lambda]4)/r3 over r3 ranging from 0 to \[Lambda]4 - 1:

(* Jacobian *)
J = 1/r3;
(* Integrand *)
integrand = FullSimplify[(density /. r -> r3) (density /. r -> (-1 - r3 + λ)/r3) J,
Assumptions -> {r3 > 0 && λ > r3 + 1}];
(* Probability density function *)
pdfλ4[λ4_] := NIntegrate[integrand /. λ -> λ4, {r3, 0, λ4 - 1}]


Now to compare with a random sample:

n = 10000;
density = 27 (r + r^2)/(8 (1 + r + r^2)^(5/2));
dist = ProbabilityDistribution[density, {r, 0, ∞}];
SeedRandom[12345];
rr3 = RandomVariate[dist, n];
rr4 = RandomVariate[dist, n];
λλ4 = 1 + rr3 + rr3*rr4;

Show[Histogram[λλ4, "FreedmanDiaconis", "PDF"],
Plot[pdfλ4[t], {t, 1, 20}]]


• we can simplify the first step using pd = ProbabilityDistribution[density, {r, 0, Infinity}]; SeedRandom[12345]; rr = RandomVariate[pd, n]
– kglr
Commented Jul 10, 2023 at 22:41
• @kglr Thanks! I'll incorporate that. I foolishly assumed that Mathematica wouldn't be able to handle getting random numbers directly from that density.
– JimB
Commented Jul 10, 2023 at 22:52
• @kglr Very good. I did just find the same result in a brute force manner. It simplifies to $1+r_3$, $1+r_3+r_3 r_4$, $1+r_3+r_3 r_4+r_3 r_4 r_5$, $1+r_3+r_3 r_4+r_3 r_4 r_5+r_3 r_4 r_5 r_6$, etc. The mean continues to increase with increasing $m$ and that is the only integer moment that exists: i.e., the variance doesn't exist even for $m=3$. And I'm more that a bit skeptical that there's a closed form for the density.
– JimB
Commented Jul 10, 2023 at 23:22
• Nice investigation! The CDF looks a bit suspicious however. One way to get spacings with $P(r)$ close to above is to let lambda be the eigenvalues of a sample from GaussianOrthogonalMatrixDistribution[10000], CDF will be a lot less peaked. Commented Jul 11, 2023 at 7:49
• Not that it isn't related but I have no idea how your comment relates to determining the distribution of the $\lambda$'s from the distribution of the spacings.
– JimB
Commented Jul 11, 2023 at 15:41