Генерация случайных рациональных чисел
Рационалы перечислимы. Например, этот код находит k-е рациональное число в открытом интервале 0..1 с упорядочиванием, что {n1, d1} предшествует {n2, d2} if (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}
Или, может быть, можно эффективно сгенерировать рациональные числа с ограниченным знаменателем, имеющим другое «однородное по ощущениям» распределение?
РЕДАКТИРОВАТЬ Это код для системы 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
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]
