随机有理数生成

Question

有理数是可数的。例如，此代码在开区间 0..1 中查找第 k 个有理数，如果 ，则 {n1, d1} 位于 {n2, d2} 之前(d1 假设 {n,d} 是互质的。

RankedRational[i_Integer?Positive] := 
 Module[{sum = 0, eph = 1, den = 1},
  While[sum < i, sum += (eph = EulerPhi[++den])];
  Select[Range[den - 1], CoprimeQ[#, den] &][[i - (sum - eph)]]/den
  ]

In[118]:= Table[RankedRational[i], {i, 1, 11}]

Out[118]= {1/2, 1/3, 2/3, 1/4, 3/4, 1/5, 2/5, 3/5, 4/5, 1/6, 5/6}

现在我想生成随机有理数，给定分母上界，均匀地排序，以便对于足够大的分母有理数将在单位区间上均匀分布。

直观上，我们可以在具有相同权重的小分母的所有有理数中进行选择：

RandomRational1[maxden_, len_] := 
 RandomChoice[(Table[
     i/j, {j, 2, maxden}, {i, 
      Select[Range[j - 1], CoprimeQ[#, j] &]}] // Flatten), len]

是否可以更有效地生成具有这种分布的随机有理数，而不需要构造所有有理数？这张桌子不需要太多就能变得很大。

In[197]:= Table[RankedRational[10^k] // Denominator, {k, 2, 10}]

Out[197]= {18, 58, 181, 573, 1814, 5736, 18138, 57357, 181380}

或者也许可以有效地生成具有不同“感觉均匀”分布的有界分母的有理数？

---

EDIT This is Mathematica code which runs acceptance-rejection generation suggested by btilly.

Clear[RandomFarey];
RandomFarey[n_, len_] := Module[{pairs, dim = 0, res, gcds},
  Join @@ Reap[While[dim < len,
      gcds = cfGCD[pairs = cfPairs[n, len - dim]];
      pairs = Pick[pairs, gcds, 1];
      If[pairs =!= {}, 
       dim += Length@Sow[res = pairs[[All, 1]]/pairs[[All, 2]]]];
      ]][[2, -1]]
  ]

以下编译函数生成整数对 {i,j} ，使得 1<=i j<=n：

cfPairs = 
  Compile[{{n, _Integer}, {len, _Integer}}, 
   Table[{i, RandomInteger[{i + 1, n}]}, {i, 
     RandomChoice[2 (n - Range[n - 1])/(n (n - 1.0)) -> Range[n - 1], 
      len]}]];

以下编译函数计算 gcd。它假设输入是一对正整数。

cfGCD = Compile[{{prs, _Integer, 1}}, Module[{a, b, p, q, mod},
    a = prs[[1]]; b = prs[[2]]; p = Max[a, b]; q = Min[a, b]; 
    While[q > 0, mod = Mod[p, q]; p = q; q = mod]; p], 
   RuntimeAttributes -> Listable];

然后

In[151]:= data = RandomFarey[12, 10^6]; // AbsoluteTiming

Out[151]= {1.5423084, Null}

In[152]:= cdf = CDF[EmpiricalDistribution[data], x];

In[153]:= Plot[{cdf, x}, {x, 0, 1}, ImageSize -> 300]

Answer 1

我强烈建议查看任意的“猜数字”游戏有理数？ 以获得有关您的根本问题的一些灵感。

如果您的目标是尽快达到大致均匀，并且您不介意选择具有不同概率的不同有理数，则以下算法应该是有效的。

lower = fractions.Fraction(0)
upper = fractions.Fraction(1)

while lower < upper:
    mid = (upper + lower)/2
    if 0 == random_bit():
        upper = largest_rational_under(mid, denominator_bound)
    else:
        lower = smallest_rational_over_or_equal(mid, denominator_bound)

请注意，这两个辅助函数都可以通过将 Stern-Brocot 树向中间移动来计算。另请注意，通过一些细微的修改，您可以轻松地将其转换为迭代算法，该算法输出一系列有理数，并最终以相等的可能性在区间内的任何位置收敛。我认为那处房产还不错。

---

如果您想要最初指定的精确分布，并且 rand(n) 为您提供从 1 到 n 的随机整数，则以下代码伪代码适用于分母绑定的 n：

Try:
    k = rand(n * (n+1) / 2)
    do binary search for largest j with j * (j-1) / 2 < k
    i = k - (j * (j-1) / 2)
    if (i, j) are not relatively prime:
        redo Try
answer = i/j

平均而言，对于较大的 n，您需要尝试大约 2.55 次。所以在实践中这应该是非常有效的。

Answer 2

由于分母有界，有理数不是均匀分布的（例如，1/2 与其他所有值之间有一个很好的间隙）。

也就是说，类似的东西

In[300]:= Rationalize[RandomReal[1, 10], 0.001]

Out[300]= {17/59, 45/68, 11/31, 9/16, 1/17, 13/22, 7/10, 1/17, 5/21, 8/39}

对你有用吗？

Answer 3

以下是对您提出的问题的一些随机想法。我没有仔细检查数学，所以我可能会在这里或那里落后 1。但这代表了我会遵循的推理。

我们只考虑区间 (0,1) 中的分数。这样就容易多了。稍后我们可以处理 1/1 和假分数。

Stern - Brocot Tree 唯一地列出了每个约简正公分式（以及每个正有理数）小于或等于一）一次，按顺序并以简化形式作为树中的节点。在这个二叉树中，任何节点以及任何分数都可以通过从最上层（为方便起见，我们将其称为层 -1）开始的有限序列的左右转弯来到达，其中包含 0/1 和 1/0。 [是的，1/0。这不是印刷错误！]

给定分母 k，您最多需要转 k 圈才能达到任何约简分数 j/k，其中 j 小于 k。例如，如果分母为 101，则分母为 101 或以下的所有可能分数将位于树中第 1 级（包含 1/1）和第 101 级（最左边位置包含 1/101）之间的某个位置。

假设我们有一个生成 0 和 1 的数字生成器。 （请不要问我该怎么做；我不知道。）Lef 任意决定 Left=0 和 Right=1。

假设我们有另一个数字生成器，可以随机生成 1 到 n 之间的整数。进一步假设生成的第一个数字是0，即。向左转：这保证分数落在区间 (0,1) 内。

选择最大分母 k。随机生成一个介于 1 和 k 之间的数字 m。然后生成 R 和 L 的随机列表。按照转弯列表遍历（即下降）Stern-Brocot 树。当您到达目的地部分时停止。

如果该分数的分母等于或小于 k，停止，这就是您的数字。

如果分母大于 k，则沿树上升（沿着下降的同一路径），直到到达分母不大于 k 的分数。

我不知道数字的生成是否真的是随机的。我什至不知道该怎么说。但就其价值而言，我没有发现任何明显的偏见来源。