Методы получения случайных комбинаторных объектов

Задача:

Необходимо сгенерировать случайный комбинаторный объект размера [math] n [/math] с равномерным распределением вероятности, если в наличии есть функция для генерации случайного числа в заданном интервале.

Описание алгоритма

Пусть [math] B = \{b_1, b_2 \ldots, b_k\} [/math] — множество различных элементов, которые могут находиться в данном комбинаторном объекте.

Будем получать элементы по порядку: сначала определим, какой элемент будет стоять на первом месте, потом на втором и так далее. Считаем, что мы построили префикс длинны [math] i [/math] : [math] P = \{a_1, a_2, \ldots, a_i\} [/math]. Будем выбирать элемент [math] a_{i+1} [/math] из множества всех возможных так, чтобы вероятность выбора элемнта [math] b \in B [/math], была пропорциональна числу комбинторных обьектов размера [math] n [/math] с префиксом [math] P + b [/math]. Для этого разобъем отрезок натуральных чисел [math] [1, s] [/math]. где [math] s [/math] — число различных комбинаторных объектов с текущим префиксом, на [math] k [/math] диапазонов так, чтобы размер диапазаона [math] d_j [/math] был равен числу объектов с префиксом [math] P + b_j [/math]. С помощью функции для генерации случайного числа получим число [math] r [/math] в интервале [math] [1, s] [/math] и добавим к префиксу [math] P [/math] элемент [math] b_j [/math] соответствующий диапазону отрезка в которм находится полученное число.

object randomObject(n: int, k: int):  // [math] n [/math] — размер комбинторного объекта, [math] k [/math] — число различных элемнтов.
  for i = 1 to n                               
    s = number(prefix)  // число комбинаторных объектов с текущим префиксом. 
    r = random(1, s)
    for j = 1 to k  
      if number(prefix + B[j]) < r  // [math] B [/math] — множество всех возможных элементов. 
        r = r - number(prefix + B[j])  // если [math] r [/math] не попало в текщий диапазон — перейдем к следующему.
      else 
        prefix[i] = b[j]
        break
  return prefix

Сложность алгоритма — [math]O(nk) [/math]. Количества комбинаторных объектов с заданными префиксами считаются известными, и их подсчет в сложности не учитывается. Стоит отметить, что подсчет количества комбинаторных объектов с заданным префиксом зачастую является задачей с достаточно большой вычислительной сложностью.

Доказательство корректности

Лемма:

Вероятность получить в процессе работы алгоритма некоторый префикс [math] P [/math] равна [math] \dfrac{S(P)}{C(n)} [/math], где [math] C(n) [/math] — число различных комбинаторных объектов данного типа длины [math] n [/math], а [math] S(P) [/math] — число различных комбинаторных обьектов длины [math] n [/math] с таким префиксом.

Доказательство:

[math]\triangleright[/math]

Докажем по индукции: База: Заметим, что любой комбинаторный объект имеет пустой префикс ([math] \varnothing [/math]), следовательно [math] S(\varnothing)=C(n) [/math]. Вероятность получить некоторый префикс [math] P [/math] длины [math] 1 [/math] равна [math] \dfrac{S(P)}{S(\varnothing)} [/math], что равно [math] \dfrac{S(P)}{C(n)} [/math] . Переход: Пусть вероятность получить префикс [math] P [/math] длины [math] l [/math] равна [math] \dfrac{S(P)}{C(n)} [/math]. Вероятность получить из [math] P [/math] некоторый префикс [math] P' [/math] длины [math] l+1 [/math] равна [math] \dfrac{S(P')}{S(P)} [/math] , следовательно вероятность получить префикс [math] P' [/math] равна [math] \dfrac{S(P)}{C(n)} [/math][math]\cdot[/math][math] \dfrac{S(P')}{S(P)} [/math] что равно [math] \dfrac{S(P')}{C(n)} [/math]. Следовательно предположение индукции верно для всех возможных префиксов любой длины.

[math]\triangleleft[/math]

Утверждение:

Результатом работы данного алгоритма является случайный комбинаторный объект размера [math] n [/math], причем вероятность получения одинакова для каждого результата.

[math]\triangleright[/math]

Заметим, что результатом работы данного алгоритма является любой из возможных префиксов размера [math] n [/math]. Любой такой префикс является комбинаторным объектом размера [math] n [/math] и наоборот, следовательно число объектов с таким префиксом равно [math] 1 [/math]. Из чего по доказанной лемме следует, что вероятность получения одинакова для всех различных результатов и равна [math] \dfrac{1}{C(n)} [/math].

[math]\triangleleft[/math]

Примеры решения задач

Рассмотрим примеры использования данного алгоритма для генерации различных типов комбинаторных объектов.

Битовые вектора

Рассмотрим алгоритм получения случайного битового вектора. В битовом векторе может находиться только два типа элементов: [math] 0 [/math] и [math] 1 [/math], следовательно [math] k = 2 [/math]. Заметим что для любого префикса длины [math] l [/math] число возможных комбинаторных объектов длины [math] n [/math] одинаково и равно [math] 2^{n-l} [/math], следовательно на каждом шаге алгоритма небходмо выбирать с равной вероятностью [math] 0 [/math] или [math] 1 [/math]

vector<int> randomBitVector(n: int):  // [math] n [/math] — размер битового вектора.
  for i = 1 to n                               
    r = random(0, 1)
    v[i] = r
  return prefix

Сложность алгоритма — [math] O(n) [/math], так как в случае двоичных векторов [math] k [/math] постоянно и равно [math] 2 [/math].

Правильные скобочные последовательности

Рассмотрим алгоритм получения случайной правильной скобочной последовательности. Правильная скобочная пследовательность состоит из двух типов элементов: открывающей и закрывающей скобок, следовательно [math] k = 2 [/math].

Рассмотрим "полуправильные" скобочные последовательности т.е. такие что всякой закрывающей скобке соответствует парная открывающая, но не все открытые скобки закрыты. Такую последовтеьность можно охарактеризовать двумя числами: [math] l [/math] — длина скобочной последовательности и [math] b [/math] — баланс (т.е. разность между количеством открывающих и закрывающих скобок). Заметим что любой префикс правильной скобочной последовательности является "полупраильной" скобочной последовательностью, и что для любого префикса [math] P [/math] длины [math] l [/math] число различных правильных скобочных последовательностей длины [math] n [/math] равно числу "полуправильных" скобочных последовательностей длины [math] n-l [/math] с таким же балансом как у [math] P [/math].

Научимся считать [math] D[l][b] [/math] — число последовательностей длины [math] l [/math] и баланса [math] b [/math]). Если [math] l = 0 [/math], то ответ понятен сразу: [math] D[0][0] = 1 [/math], все остальные [math] D[0][b] = 0 [/math]. Пусть теперь [math] l \gt 0 [/math], тогда переберём, чему мог быть равен последний символ этой последовательности. Если он был равен [math] '(' [/math], то до этого символа мы находились в состоянии [math] (l-1,b-1) [/math]. Если он был равен [math]')'[/math], то предыдущим было состояние [math](l-1,b+1)[/math]. Таким образом, получаем формулу:

[math]D[l][b] = D[l-1][b-1] + D[l-1][b+1][/math]

(считается, что все значения [math] D[l][b] [/math] при отрицательном [math]j[/math] равны нулю). Этот преподсчет можно выполнить за [math]O(n^2)[/math].

Будем строить префикс следующим образом: на каждом шаге интервал случайных чисел [math] [0, s] [/math] (где [math] s = D[n - l][b] [/math]) , будет разбиваться на два диапазона размерами [math] D[n - l - 1][b + 1] [/math] и [math] D[n - l - 1][b - 1] [/math] , и к префиксу будет прибавляться [math]'('[/math] или [math]')'[/math] если случайное число попало в первый или второй диапазон соответственно.

string randomBracketSequence(n: int):  // [math] n [/math] — длина скобочной последовательности. 
  b = 0
  l = 0
  for i = 1 to n                               
    s = D[n - l][b] 
    r = random(1, s)
     if D[n - l - 1][b + 1] >= r
       l = l + 1
       b = b + 1
       result = result + '('
     else
       l = l + 1
       b = b - 1
       result = result + ')'
  return result

Итоговая сложность алгоритма — [math] O(n) + O(n^2) [/math] на преподсчет.

Разбиения на множества

Разбиение на [math] k [/math] подмножеств

Рассмотрим множество первых [math] n [/math] натуральных чисел: [math] N_n = \{1, 2, \ldots, n\} [/math]. Необходимо разбить его на [math] k [/math] непустых подмножеств [math] \{B_1, B_2, \ldots, B_k\} [/math] с равным распределением вероятности.

Будем строить разбиение таким образом, чтобы в результате подмножества [math] \{B_1, B_2, \ldots, B_k\} [/math] оказались отсортированы в лексикографическом порядке (т.е. чтобы для любых [math]i, j \mid 1 \leqslant i \lt j \leqslant k [/math] наименьший элемент [math] B_i [/math] был меньше наименьшего элемента [math] B_j [/math]). Для этого будем по очереди добавлять каждое число от [math] n [/math] до [math] 1 [/math] в одно из подмножеств и для каждого из подмножеств начиная с [math] B_n [/math] и заканчивая [math] B_1 [/math] будем выбирать какой элемент будет добавлен в него последним(т.е. будет минимальным).

На каждом шаге префиксом считаем текущее разбиение. Оно характеризуется двумя значениями: [math] l [/math] — число добавленных элементов и [math] m [/math] — число подмножеств для которых определен последний элемент. Заметим, что количество разбиений на подмножества с заднным префиксом равно числу способов разбить еще не добавленные элементы на еще не законченные подмножества так, чтобы они оказались лексикографически упорядочены, то есть равно числу разбиений [math] n-l [/math] элементов на [math] k-m [/math] непустых подмножеств, что равно [math]\lbrace{n-l\atop k-m}\rbrace[/math] (т.е числу Стирлинга второго рода). Таким образом из префикса [math] P [/math] мы можем получить следующий префикс [math] P' [/math] двумя способами:

• Добавить текущий элемент ([math] n-l [/math]) в одно из [math] k-m [/math] незаконченных подмножеств. В таком случае число обьектов с префиксом [math] P' [/math] будет равно [math]\lbrace{n-l-1\atop k-m}\rbrace[/math] .
• Сделать текущий элемент последним в подмножестве [math] B_{k-m} [/math] . В таком случае это подмножество станет законченым, следовательно число обьектов с префиксом [math] P' [/math] будет равно [math]\lbrace{n-l-1\atop k-m-1}\rbrace[/math].

Таким образом на каждом шаге интервал случайных чисел [math] [0, s] [/math] (где [math] s = [/math][math]\lbrace{n-l\atop k-m}\rbrace[/math]) , будет разбиваться на два диапазона размерами [math]\lbrace{n-l-1\atop k-m-1}\rbrace[/math] и [math] (k-m)\cdot [/math][math]\lbrace{n-l-1\atop k-m}\rbrace[/math]. Если случайно сгенерированное число попадет в первый диапазон, то сделаем [math] n-l [/math] последним элементом в подмножестве [math] B_{k-m} [/math] . Иначе добавим [math] n-l [/math] в случайно выбранное из незаконченных подмножеств ([math] \{B_1, B_2, \ldots, B_{k-m}\} [/math]).

int[] randomSetPartition(n: int k: int):  // [math] n [/math] — количество элементов в множестве, [math] k [/math] — число подмножеств на которые нужно разбить исходное множество. 
  l = 0
  m = 0
  downto i = n to 1                               
    s = stirling(n - l, k - m)  // [math] stirling(a, b) [/math] — функция возвращающая число Стирлинга второго рода для заданных аргументов. 
    r = random(1, s)
     if stirling(n - l - 1, k - m - 1) >= r
       result[i] = k - m  // [math] result[i] [/math] — номер подмножества в котором находится элемент [math] i [/math]. 
       l = l + 1
       m = m + 1
     else
       p = random(1, k - m)  // Случайным образом выбираем номер одного из незаконченных подмножеств. 
       result[i] = p
       l = l + 1
  return result

Так как на каждом шаге интервал случайных чисел разделяется только на на два диапазона, а всего шагов — [math] n [/math] то итоговая сложность алгоритма — [math] O(n) + O(n^2) [/math] на преподсчет чисел Стирлинга второго рода (если преподсчитать их динамически).

Разбиение на случайное число подмножеств

Описаный алгоритм можно применить для получения разбиения множества на случайное число подмножеств. Для этого достаточно случайным образом выбрать число подмножеств из интевала [math] [1, n] [/math] , так чтобы вероятность получить число [math] k [/math] была пропорциональна числу способов разбить [math] n [/math] элементов на [math] k [/math] подмножеств (что равно [math]\lbrace{n\atop k}\rbrace[/math]).

Разделим интервал случайных чисел [math] [1, s] [/math] (где [math] s = [/math] [math] \sum\limits_{i=1}^{n}\left\{{n\atop i}\right\}[/math]) на [math] n [/math] диапазонов, так чтобы размер диапазона [math] d_i [/math] был равен [math] \lbrace{n\atop i}\rbrace [/math]. С помощью функции для генерации случайного числа получим число [math] r [/math] в интервале [math] [1, s] [/math] и выберем количество подмножеств [math] k [/math] соответствующее диапазону отрезка в которм находится полученное число и по выбранному [math] k [/math] получим случайное разбиение множества размера [math] n [/math] на [math] k [/math] подмножеств.

См. также

• Метод генерации случайной перестановки
• Методы генерации случайного сочетания

Источники информации

• Комбинаторные алгоритмы: учебное пособие / Т. И. Федоряева. - Новосибирск: Новосибирский гос. ун-т, 2011. 118 с. - ISBN 978-5-4437-0019-9
• Non-Uniform Random Variate Generation / Luc Devroye. - Springer, New York, NY, 1986. 843 c. - ISBN 978-1-4613-8645-2