理性是可以列举的。例如,此代码在打开间隔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 < 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]
最佳答案
我强烈建议您查看The "guess the number" game for arbitrary rational numbers?以获取有关您的潜在问题的一些启发。
如果您的目标是尽快实现统一,并且您不介意选择具有不同概率的不同有理数,则以下算法应该是有效的。
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)
请注意,可以通过将斯特恩-布罗科树向中间走去来计算这两个辅助函数。还要注意,通过一些小的修改,您可以轻松地将其转换为迭代算法,该算法吐出有理数序列,并最终在区间中的任何地方均等收敛。我认为该物业不错。
如果要获得最初指定的确切分布,并且
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倍。因此,在实践中这应该非常有效。关于math - 随机有理数生成,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/5720689/