随机有理数生成

Question

理性是可以枚举的.例如,这个代码在开放区间0..1中找到第k个有理数,如果假设是互质的,{n1, d1}则排序在之前.{n2, d2}(d1<d2 || (d1==d2 && n1<n2)){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}

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

---

编辑这是Mathematica代码,它运行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你得Try约2.55倍.所以在实践中这应该非常有效.