随机有理数生成 [英] Random rational numbers generation

问题描述

比率是可枚举的.例如，这段代码在打开间隔0 ... 1中找到第k个有理数，如果(d1<d2 || (d1==d2 && n1<n2))假设{n,d}是互素的，则命令{n1, d1}在{n2, d2}之前.

Rationals are enumerable. For example this code finds k-th rational in open interval 0..1, with ordering that {n1, d1} is before {n2, d2} if (d1<d2 || (d1==d2 && n1<n2)) assuming {n,d} is coprime.

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}

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

Now I would like to generate random rationals, given an upper bound on the denominator sort-of uniformly, so that for large enough denominator rationals will be uniformly distributed over the unit interval.

从直觉上讲，人们可以在所有具有相同权重的小分母的有理数中进行选择:

Intuitively, one could pick among all rationals with small denominators with equal weights:

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

在不构建所有有理数的情况下，能否更有效地生成具有这种分布的随机有理数?这张桌子变大不需要什么.

Can one generate random rationals with this distribution more efficiently, without constructing all of them ? It does not take much for this table to become huge.

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

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

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

Or maybe it is possible to efficiently generate rationals with bounded denominator having a different "feels-like-uniform" distribution ?

---

编辑:这是运行btilly建议的接受/拒绝生成的Mathematica代码.

---

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:

The following compiled function generations pairs of integers {i,j} such that 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.假定输入是一对正整数.

and the following compiled function computes gcd. It assumes the input is a pair of positive integers.

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]

推荐答案

我强烈建议您查看有关您潜在问题的一些启发.

I strongly suggest looking at The "guess the number" game for arbitrary rational numbers? for some inspiration about your underlying problem.

如果您的目标是尽快实现统一，并且您不介意选择具有不同概率的不同有理数，则以下算法应该是有效的.

If your goal is to be approximately uniform ASAP, and you don't mind picking different rationals with different probabilities, the following algorithm should be efficient.

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)

请注意，可以通过向中间行走斯特恩-布罗科树来计算这两个辅助函数.还要注意，通过一些小的修改，您可以轻松地将其转换为迭代算法，该算法吐出有理数序列，并最终在区间中的任何位置以相等的可能性收敛.我认为该属性不错.

Note that both of those two helper functions can be calculated by walking the Stern-Brocot Tree towards the mid. Also note that, with some minor modification, you can easily transform this into an iterative algorithm that spits out a sequence of rational numbers, and eventually will converge with equal likelihood anywhere in the interval. I consider that property to be kind of nice.

如果您希望获得最初指定的确切分布，并且rand(n)为您提供一个从1到n的随机整数，那么下面的伪代码将适用于分母绑定的n:

If you want the exact distribution that you originally specified, and rand(n) gives you a random integer from 1 to n, then the following pseudocode will work for denominator bound 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倍.因此，在实践中，这应该非常有效.

On average for large n you'll have to Try about 2.55 times. So in practice this should be pretty efficient.

这篇关于随机有理数生成的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！