Генерация случайных рациональных чисел

Question

Рациональные числа перечисляются.Например, этот код находит k-тое рациональное значение в открытом интервале 0..1 с порядком, в котором {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}

Или, возможно, можно эффективно генерировать рациональные числа с ограниченным знаменателем, имеющим другое распределение "по ощущениям, как единообразное"?

---

EDIT Это код 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)

Обратите внимание, что обе эти вспомогательные функции можно рассчитать, если подойти к Дереву Штерн-Броко к середине. Также обратите внимание, что с некоторой незначительной модификацией вы можете легко преобразовать это в итерационный алгоритм, который выделяет последовательность рациональных чисел и в конечном итоге будет сходиться с равной вероятностью в любом месте интервала. Я считаю, что это было бы неплохо.

---

Если вам нужно точное распределение, которое вы изначально указали, а 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 раза. Так что на практике это должно быть довольно эффективно.

Answer 2

С ограничением на знаменатель, рациональные числа распределены неравномерно (например, 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}

работаешь на тебя?

Answer 3

Вот несколько случайных мыслей о проблеме, которую вы поднимаете.Я не тщательно проверил математику, чтобы я мог быть на 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.Произвольно сгенерируйте число m между 1 и k.Затем сгенерируйте случайный список R и L.Пройдите (т.е. спуститесь) по дереву Штерна-Броко, следуя списку поворотов.Остановитесь, когда вы достигнете целевой фракции.

Если эта дробь имеет знаменатель, равный или меньший k, остановитесь, это ваш номер.

Если знаменатель больше k, поднимайтесь по дереву (по тому же пути, по которому вы спустились), пока не достигнете дроби со знаменателем не более k.

Я не знаю, что генерация чисел действительно случайна.Я даже не знаю, как сказать.Но что бы это ни значило, я не обнаружил очевидного источника предвзятости.