随机有理数生成

18

有理数是可枚举的。例如,这段代码可以找到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}

也许有可能高效地生成分母受限的有不同“感觉均匀”分布的有理数?


编辑这是运行btilly建议的接受-拒绝生成的Mathematica代码。

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:
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]}]];

以下编译函数计算最大公约数,假设输入为两个正整数。

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]

在此输入图片描述


2
问题在于费雷序列只有一个递归表达式... - Dr. belisarius
http://oeis.org/A005728/a005728.nb - Dr. belisarius
@belisarius 所以我应该把它称为RandomFarey。可以通过随机遍历二进制树从Farey序列中生成随机元素吗?这可以在不存储内存的情况下完成吗? - Sasha
不是我目前所知的。但我的数论知识有些生疏。你可以尝试在MathOverflow上发问,我猜。 - Dr. belisarius
3个回答

8

我强烈建议参考“猜数字”游戏,适用于任意有理数?来获得关于您底层问题的灵感。

如果您的目标是尽快实现近似均匀分布,并且不介意选择不同的有理数具有不同的概率,则以下算法应该是有效的。

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)会给您一个介于1n之间的随机整数,那么以下伪代码将适用于分母限制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次。因此在实践中,这应该相当有效。


谢谢您。这正是我想要的。后一个算法更适合在Mathematica中进行向量化操作。据我所理解,您的算法是简单地选择具有相等概率的 (i,j)对,然后拒绝它们。我唯一的担心是计算二次项时可能会出现整数溢出问题。不过我知道如何解决这个问题。我将用Mathematica实现您的算法,并更新我的问题。 - Sasha

3

由于分母有限制,有理数并不均匀分布(例如,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}

你需要我为你工作吗?


是的,但随着分母的增加,有理数越来越接近。我的目标是生成有理数,使得样本上的EmpiricalDistributionCDF在分母边界变大时趋于直线。考虑cdf = CDF[EmpiricalDistribution[RandomRational1[50, 10^4]], x]; Plot[{x, cdf}, {x, 0, 1}, Exclusions -> None, PlotPoints -> 250, PlotStyle -> {Black, Orange}]。此外,Cramer von-Mises测试表明分布接近均匀分布Histogram[ Table[CramerVonMisesTest[RandomRational1[100, 10^3], UniformDistribution[]], {10^3}]] - Sasha
很遗憾,无法在评论中使用图像。评估 cdf1 = CDF[EmpiricalDistribution[RandomRational1[25, 10^4]], x];cdf2 = CDF[EmpiricalDistribution[Rationalize[RandomReal[1, 10^4], 1/25]], x]; Row@{Plot[{cdf1, x}, {x, 0, 1}, ImageSize -> 250], Plot[{cdf2, x}, {x, 0, 1}, ImageSize -> 250]} 您可以看到,有理化的随机实数比从 RandomRational1 生成的实数要慢得多。 - Sasha
@ Sasha RandomRational1 生成的不同有理数比Rationalize方法多约5倍... - Brett Champion

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的数字生成器。(请不要问我如何做到这一点;我不知道。)随意决定Left=0,Right=1。

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

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

如果该分数的分母等于或小于k,则停止,那就是您要的数字。

如果分母大于k,则沿着下降的路径向上爬树,直到找到一个分母不大于k的分数。

我不知道数字生成是否真正随机。我甚至不知道如何判断。但就价值而言,我没有发现任何明显的偏差来源。


这个算法可以工作,但不会收敛到任何接近均匀分布的东西。 - btilly
@btilly 请解释一下。它会偏袒某些分数吗? - DavidC
@David,[L,L,R,R,L](或任何人)在树中的位置公式是什么? - Dr. belisarius
@belisarius 你会在这个演示的源代码中找到它。(下载源代码; 它不是太长)。http://demonstrations.wolfram.com/FractionTreeAndContinuedFractions/左右都是矩阵。Knuth的《具体数学》为我开发的代码提供了基本原理。我可以解释一下,但我认为你可以从演示中的代码更快地理解它。 - DavidC
1
@david-carraher:它会强烈偏向某些地区。例如,只有在前10个选择为“L”的情况下,才能到达从0到0.05的区域,这种情况只发生约0.1%的时间。 - btilly

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接