
Однажды я задал на Stack Overflow вопрос о структуре данных для шулерских игральных костей. В частности, меня интересовал ответ на такой вопрос: «Если у нас есть n-гранная кость, у грани которой i есть вероятность выпадения pi. Какова наиболее эффективная структура данных для симуляции бросков такой кости?»
Такую структуру данных можно использовать для многих задач. Например, можно применять её для симуляции бросков честной шестигранной кости, присвоив вероятность
Если выйти за пределы игр, то можно применить эту структуру данных в симуляции роботов, датчики которых имеют известные уровни отказа. Например, если датчик дальности имеет 95-процентную вероятность возврата правильного значения, 4-процентную вероятность слишком маленького значения, и 1-процентную вероятность слишком большого значения, то можно использовать эту структуру данных для симуляции считывания показаний датчика генерацией случайного результата и симуляцией считывания датчиком этого результата.
Полученный мной в Stack Overflow ответ впечатлил меня по двум причинам. Во-первых, в решении мне посоветовали использовать мощную технику под названием alias-метод, которая при определённых разумных предположениях о модели машины способна после простого этапа предварительной подготовки симулировать броски кости за время
Эта статья является моей попыткой сделать краткий обзор различных подходов к симуляции шулерской кости, от простых и очень непрактичных техник, до очень оптимизированного и эффективного alias-метода. Я надеюсь, что мне удастся передать различные способы интуитивного понимания задачи и то, насколько каждый из них подчёркивает какой-то новый аспект симуляции шулерской кости. Моя цель для каждого подхода заключается в исследовании мотивирующей идеи, базового алгоритма, доказательства верности и анализа времени выполнения (с точки зрения требуемых времени, памяти и случайности).
Вступление
Прежде чем переходить к конкретным подробностям различных техник, давайте сначала стандартизируем терминологию и систему обозначений.
Во введении к статье я использовал термин «шулерская кость» для описания обобщённого сценария, в котором существует конечное множество результатов, с каждым из которых связана вероятность, Формально это называется дискретным распределением вероятностей, а задача симуляции шулерской кости называется выборкой из дискретного распределения.
Для описания нашего дискретного распределения вероятностей (шулерской кости) мы будем считать, что у нас есть множество из n вероятностей
Работа с вещественными числами на компьютере — это «серая область» вычислений. Существует множество быстрых алгоритмов, скорость которых обеспечивается исключительно способностью за постоянное время вычислять функцию floor произвольного вещественного числа, и численные неточности в представлении чисел с плавающей запятой могут полностью разрушить некоторые алгоритмы. Следовательно, прежде чем приступать к каким-либо обсуждениям алгоритмов, работающих с вероятностями, то есть вступать в тёмный мир вещественных чисел, я должен уточнить, что может и чего не может компьютер.
Здесь и далее я буду предполагать, что все указанные ниже операции могут выполняться за постоянное время:
- Сложение. вычитание, умножение, деление и сравнение произвольных вещественных чисел. Нам необходимо будет это делать для манипуляций с вероятностями. Может показаться, что это слишком смелое предположение, но если считать, что точность любого вещественного числа ограничена неким многчленом размера слова машины (например, 64-битное double на 32-битной машине), но я не думаю, что это слишком неразумно.
- Генерирование равномерного вещественного числа в интервале [0, 1). Для симуляции случайности нам нужен некий источник случайных значений. Я предполагаю, что мы можем за постоянное время генерировать вещественное число произвольной точности. Это намного превышает возможности реального компьютера, но мне кажется, что в целях этого обсуждения такое допустимо. Если мы согласимся пожертвовать долей точности, сказав, что произвольное IEEE-754 double находится в интервале [0, 1], то мы и в самом деле потеряем точность, но результат, вероятно, будет достаточно точным для большинства применений.
- Вычисление целочисленного floor (округления в меньшую сторону) вещественного числа. Это допустимо, если мы предположим, что работаем с IEEE-754 double, но в общем случае такое требование к компьютеру невыполнимо.
Стоит задаться вопросом — разумно ли считать, что мы можем выполнять все эти операции эффективно. На практике мы редко используем вероятности, указываемые до такого уровня точности, при котором свойственная IEEE-754 double ошибка округления может вызвать серьёзные проблемы, поэтому выполнить все представленные выше требования мы можем, просто работая исключительно с IEEE double. Однако если мы находимся в среде, где вероятности указываются точно как рациональные числа высокой точности, то подобные ограничения могут оказаться неразумными.
Симуляция честной кости
Прежде чем мы перейдём к более общему случаю бросания произвольной шулерской кости, давайте начнём с более простого алгоритма, который послужит строительным блоком для последующих алгоритмов: с симуляции честной n-гранной кости. Например, нам могут пригодиться броски честной 6-гранной кости при игре в Monopoly или Risk, или бросание честной монетки (двусторонней кости) и т.д.
Для этого конкретного случая есть простой, элегантный и эффективный алгоритм симуляции результата. В основе алгоритма лежит следующая идея: предположим, что мы можем генерировать действительно случайные, равномерно распределённые вещественные числа в интервале

Теперь если мы хотим бросить

Далее мы генерируем случайно выбираемое вещественное число в интервале

то мы можем сказать, что на кости выпало 2 (если считать, что грани кости проиндексованы с нуля).
Графически легко увидеть, в какую из областей попало случайное значение, но как нам закодировать это в алгоритме? И здесь мы воспользуемся тем фактом, что это честная кость. Поскольку все интервалы имеют равный размер, а именно
Алгоритм: симуляция честной кости
- Генерируем равномерно распределённое случайное значение
в интервале
.
- Возвращаем
.
Учитывая наши сделанные выше допущения о вычислениях, этот алгоритм выполняется за время.
Из этого раздела можно сделать два вывода. Во-первых, мы можем разделить интервал
Симуляция шулерской кости при помощи честной кости
Имея алгоритм симуляции честной кости, можем ли мы адаптировать его для симуляции шулерской кости? Интересно, что ответ положительный, но решение потребует больше пространства.
Из предыдущего раздела интуитивно понятно, что для симуляции броска шулерской кости нам достаточно разделить интервал

К сожалению, на этом шаге мы застрянем. Даже если бы мы знали случайное число в интервале
Однако мы можем адаптировать технику, используемую для честной кости, чтобы она работала и в этом случае. Давайте возьмём для примера рассмотренную выше кость. Вероятность выпадения граней равна

Тогда мы назначим их различным результатам следующим образом:

Теперь симулировать бросок кости будет очень просто — мы просто бросаем эту новую честную кость, а затем смотрим, какая грань выпала и считываем её значение. Этот первый шаг можно выполнить представленным выше алгоритмом, что даст нам целое числов в интервале

Чтобы формализировать это в виде алгоритма, мы опишем и этап инициализации (получение таблицы), и этап генерации (симуляцию броска случайной кости). Оба эти этапа важно учитывать в этом и последующих алгоритмах, потому что время подготовки должно быть отличным.
На этапе инициализации мы начинаем с поиска наименьшего общего кратного всех вероятностей, заданных для граней кости (в нашем примере НОК равно 12). НОК полезно здесь, потому что оно соответствует наименьшему общему делителю, который мы можем использовать для всех дробей, а значит, и количеству граней новой честной кости, которую мы будем бросать. Получив это НОК (обозначим его L), мы должны определить, сколько граней новой кости будет распределено на каждую из граней исходной шулерской кости. В нашем примере грань с вероятностью
Вот псевдокод представленного выше алгоритма:
Алгоритм: симуляция шулерской кости при помощи честной кости
- Инициализация:
- Находим НОК знаменателей вероятностей
; обозначаем его
- Выделяем массив
размера
для сопоставления результатов бросков честной кости с бросками исходной кости.
- Для каждой грани
исходной кости выполняем в любом порядке следующее:
- Присваиваем следующим
элементам
значение
.
- Генерирование:
- Генерируем бросок честной кости для
-гранной кости; называем грань
.
- Возвращаем
.
Возможно, этот алгоритм и прост, но насколько он эффективен? Само генерирование бросков кости достаточно быстрое — каждый бросок кости требует
Однако этап инициализации может быть чрезвычайно затратным. Чтобы этот алгоритм заработал, нам нужно выделить пространство под массив размером с НОК знаменателей всех входных дробей. В нашем примере (
К сожалению, всё может быть ещё хуже. В предыдущем примере мы по крайней мере можем «ожидать», чтоб алгоритм займёт много памяти, поскольку оба знаменателя дробей равны одному миллиону. Однако у нас может появиться множество вероятностей, для которых НОК значительно больше, чем каждый отдельный знаменатель. Для примера давайте рассмотрим вероятности
Иными словами, во многих случаях этот алгоритм ведёт себя хорошо. Если все вероятности одинаковы, то все получаемые на входе вероятности равны

Это даёт нам следующую информацию об алгоритме:
Алгоритм | Время инициализации | Время генерации | Занимаемая память | |||
---|---|---|---|---|---|---|
Наилучшее | Наихудшее | Наилучшее | Наихудшее | Наилучшая | Наихудшая | |
Шулерская кость из честной кости |
Ещё одна важная подробность об этом алгоритме: он предполагает, что мы будем получать удобные вероятности в виде дробей с хорошими знаменателями. Если вероятности задаются как IEEE-754 double, то такой подход скорее всего ждёт катастрофа из-а небольших ошибок округления; представьте, что у нас есть вероятности 0,25 и 0,250000000001! Поэтому такой подход, вероятно, лучше не использовать, за исключением особых случаев, когда вероятности ведут себя хорошо и указаны в формате, соответствующем операциям с рациональными числами.
Симуляция несимметричной монеты
Наше объяснение простого случайного примитива (честной кости) привело к простому, но потенциально ужасно неэффективному алгоритму симуляции шулерской кости. Возможно, изучение других простых случайных примитивов прольёт немного света на другие подходы к решению этой задачи.
Простая, но на удивление полезная задача заключается в симуляции с помощью генератора случайных чисел несимметричной монеты. Если у нас есть монета с вероятностью выпадения орла
Ранее мы выработали один интуитивно понятный подход: разбиение интервала

А потом сгенерировать равномерно распределённое случайное значение в интервале
Алгоритм: симуляция несимметричной монеты
- Генерируем равномерно распределённое случайное значение в интервале
.
- Если
, возвращаем «орёл».
- Если
, возвращаем «решка».
Поскольку мы можем сгенерировать равномерно распределённое случайное значение в интервале
Симуляция честной кости с помощью несимметричных монет
Из предыдущего обсуждения мы знаем, что можно симулировать шулерскую кость с помощью честной кости, если предположить, что мы готовы потратить лишнее место в памяти. Поскольку мы можем воспринимать несимметричную монету как шулерскую двустороннюю кость, то это значит, что мы можем симулировать несимметричную монету с помощью честной кости. Интересно, что можно сделать и обратное — симулировать честную кость с помощью несимметричной монеты. Конструкция получается простой, элегантной и её можно легко обобщить для симуляции шулерской кости с помощью множества несимметричных монет.
Конструкция для симуляции несимметричной монеты разбивает интервал

Теперь предположим, что нас интересует симуляция броска этой честной кости с помощью набора несимметричных монет. Один из способов решения заключается в следующем: представьте, что мы обходим эти области слева направо, каждый раз спрашивая, хотим ли мы останавливаться в текущей области, или будем двигаться дальше. Например, давайте предположим, что мы хотим выбрать случайным образом одну из этих областей. Начиная с самой левой области, мы будем подбрасывать несимметричную монету, которая сообщает нам, должны ли мы остановиться в этой области, или продолжить движение дальше. Поскольку нам нужно выбирать из всех этих областей равномерно с вероятностью
Если монеты падают решкой вверх, то мы оказывается во второй области и снова спрашиваем, должны ли мы снова выбрать эту область или продолжать движение. Вы можете подумать, что для этого мы должны бросить ещё одну монету с вероятностью выпадения орла
Чтобы определить вероятность, с которой наша несимметричная монета должна выбрасывать орла после пропуска первой области, нам нужно заметить, что после пропуска первой области их осталось всего три. Поскольку мы бросаем честную кость, нам нужно, чтобы каждая из этих трёх областей была выбрана с вероятностью
Такое интуитивное понимание приводит нас к следующему алгоритму. Заметьте, что мы не обсуждали правильность или ошибочность этого алгоритма; вскоре мы этим займёмся.
Алгоритм: симуляция честной кости с помощью несимметричных монет
- For
to
:
- Бросаем несимметричную монету с вероятностью выбрасывания орла
.
- Если выпадает орёл, то возвращаем
.
Этот алгоритм прост и в наихудшем случае выполняется за время
Теорема: приведённый выше алгоритм возвращает сторонус вероятностью
для любого выбранного
.
Доказательство: рассмотрим любое постоянное. Мы доказываем с помощью сильной индукции, что каждая из
граней имеет вероятность выбора
.
Для нашего примера мы показываем, что гранькости имеет вероятность выбора
. Но это непосредственно следует из самого алгоритма — мы выбираем грань 0, если на несимметричной монете с вероятностью выпадения орла
выпадает орёл, то есть мы выбрали её с вероятностью
.
Для этапа индукции предположим для граней, что эти грани выбираются с вероятностью
и рассмотрим вероятность выбора грани
. Грань
будет выбрана тогда и только тогда, когда не выбраны первые
граней, а на монете с вероятностью выпадения орла
выпал орёл. Поскольку каждая из первых
граней имеет вероятность выбора
, и поскольку выбирается только одна грань, то вероятность выбора одной из первых граней
задаётся как
. Это значит, что вероятность того, что алгоритм не выберет одну из первых
граней равен
. То есть вероятность выбора грани
задаётся как
, что и требуется показать. Таким образом, каждая грань кости выбирается равномерно и случайно.
Разумеется, алгоритм довольно неэффективен — с помощью первой техники мы можем симулировать бросок честной кости за время
Симуляция шулерской кости с помощью несимметричных монет
Представленный выше алгоритм интересен тем, что даёт нам простой каркас для симуляции кости с помощью набора монет. Мы начинаем с броска монеты для определения того, нужно ли выбрать первую грань кости или двигаться к оставшимся. При этом процессе нам нужно внимательно обращаться с изменением масштаба оставшихся вероятностей.
Давайте посмотрим, как можно использовать эту технику для симуляции броска шулерской кости. Воспользуемся нашим примером с вероятностями

Теперь давайте подумаем о том, как можно симулировать такую шулерскую кость с помощью несимметричных монет. Мы можем начать с броска монеты с вероятностью выпадения орла
Подведём итог — вероятности монет будут следующими:
- Первый бросок:
- Второй бросок:
- Третий бросок:
- Четвёртый бросок:
Может быть интуитивно понятно, откуда берутся эти числа, но чтобы превратить подбор в алгоритм, нам придётся создать формальную конструкцию выбора вероятностей. Идея будет следующей — на каждом этапе мы запоминаем оставшуюся часть массы вероятностей. В начале, до броска первой монеты, она равна
Алгоритм: шулерская кость из несимметричных монет
- Инициализация:
- Сохраняем вероятности
для дальнейшего использования.
- Генерирование:
Задаём
Forto
:
- Бросаем несимметричную монету с вероятностью выпадения орла
.
- Если выпал орёл, то возвращаем
.
- В противном случае задаём
С интуитивной точки зрения это логично, но верно ли это математически? К счастью, ответ положительный благодаря обобщению вышеприведённого доказательства:
Теорема: показанный выше алгоритм возвращает граньс вероятностью
при любом выбранном
.
Доказательство: рассмотрим любое постоянное. С помощью сильной индукции мы доказываем, что каждая из
граней имеет вероятность выбора
.
В качестве базового случая докажем, что гранькости имеет вероятность выбора
. Мы выбираем грань
, если при броске самой первой монеты получаем орла, что происходит с вероятностью
. Поскольку
изначально равна
, то эта вероятность равна
, то есть грань 0 выбирается с вероятностью
, как и требуется.
Для этапа индукции предположим, чтоб гранивыбираются с вероятностью
и рассмотрим вероятность выбора ребра
. Грань
будет выбрана тогда и только тогда, когда не выбраны первые
граней, после чего на монете с вероятностью выпадения орла
выпадает орёл. Поскольку каждая из первых
граней была выбрана с правильной вероятностью, а выбирается всегда только одна грань, то вероятность выбора одной из первых
задаётся как
. Это означает, что вероятность не выбора алгоритмом одной из первых
граней задаётся как
. Вероятность того, что монета для грани
выпадет орлом, задаётся как
и после
итераций мы можем при помощи быстрой индукции увидеть, что
. Это означает, что общая вероятность выбора грани
задаётся как
, как и требовалось.
Теперь давайте оценим временную сложность этого алгоритма. Мы знаем, что время инициализации может быть
Однако после размышлений становится понятно, что на количество необходимых бросков монет сильно влияет входящее распределение. В самом наилучшем случае у нас будет распределение вероятностей, при котором вся масса вероятностей сосредоточена на первой грани кости, а все остальные вероятности равны нулю. В таком случае нам достаточно одного броска монеты. В самом наихудшем случае вся масса вероятностей сосредоточена в самой последней грани кости, а на всех остальных гранях равна нулю. В таком случае нам для получения результата придётся бросать
Мы можем чётким и математическим образом характеризовать ожидаемое количество бросков монеты в этом алгоритме. Давайте представим случайную переменную
Каково же значение
Заметьте, что в последнем упрощении первый член эквивалентен
Алгоритм | Время инициализации | Время генерации | Занятая память | |||
---|---|---|---|---|---|---|
Наилучшее | Наихудшее | Наилучшее | Наихудшее | Наилучшая | Наихудшая | |
Шулерская кость из честной кости | ||||||
Шулерская кость из несимметричных монет |
Обобщение несимметричных монет: симуляция шулерской кости
В показанном выше примере мы смогли эффективно симулировать несимметричную монету, поскольку нам необходимо было учитывать только одну точку разбиения. Как нам эффективно обобщить эту идею до шулерской кости, в которой количество граней может быть произвольным?
Как можно заметить, несимметричная монета — это шулерская кость, только с двумя гранями. Следовательно, мы можем воспринимать несимметричную монету просто как особый случай более общей задачи, которую хотим решить. При решении задачи несимметричной монеты мы разделяем интервал

Заметьте, где расположены эти области. Первая область начинается с
Заметьте, что разность между этими двумя значениями равна
Теперь мы знаем, где находятся области. Если мы хотим выбрать равномерно распределённое случайное значение

Из чего мы можем заключить, что на кости с индексацией с нуля выпала грань 3.
Такой алгоритм линейного сканирования даст нам алгоритм времени

Это даёт нам алгоритм со временем
Этот алгоритм иногда называют алгоритмом выбора колесом рулетки, потому что он выбирает случайную область с помощью техники, похожей на колесо рулетки — броском шарика в интервал и наблюдением за тем, где он остановится. В псевдокоде алгоритм выглядит так:
Алгоритм: выбор колесом рулетки
- Инициализация:
- Выделяем массив
размера
- Задаём
.
- Для каждой вероятности
от
до
:
- Задаём
- Генерация:
- Генерируем равномерно распределённое случайное значение
в интервале
- С помощью двоичного поиска находим индекс
наименьшего элемента
, который меньше
.
- Возвращаем
.
Сравнение между этим алгоритмом и приведённым ранее выглядит довольно впечатляюще:
Алгоритм | Время инициализации | Время генерации | Занятая память | |||
---|---|---|---|---|---|---|
Наилучшее | Наихудшее | Наилучшее | Наихудшее | Наилучшая | Наихудшая | |
Шулерская кость из честной кости | ||||||
Шулерская кость из несимметричных монет | ||||||
Выбор колесом рулетки |
Очевидно, что теперь у нас есть гораздо лучший алгоритм, чем изначальный. Дискретизация вероятностей только поначалу казалась многообещающей, но этот новый подход, основанный на непрерывном значении и двоичном поиске, выглядит гораздо лучше. Однако всё ещё возможно улучшить эти показатели с помощью умного использования набора гибридных техник, которые мы рассмотрим чуть ниже.
Интересная деталь этого алгоритма заключается в том, что хотя использование двоичного поиска гарантирует время наихудшего времени генерации случайных чисел
Допустим мы перейдём от двоичного поиска по списку накапливаемых вероятностей к использованию двоичного дерева поиска. Например, имея приведённое выше множество вероятностей, мы можем построить для их накапливаемого распределения следующее двоичное дерево поиска:

Теперь, если мы захотим симулировать бросок кости, то сможем сгенерировать равномерно распределённое число в интервале
Однако если предположить, что мы знаем больше о распределении вероятностей, то можно сделать гораздо лучше. Например, предположим, что наши вероятности равны

Хотя это двоичное дерево поиска идеально сбалансировано, для наших задач оно не очень подходит. Поскольку мы знаем, что в 99 из 100 случаев случайное значение будет находиться в интервале

Теперь мы скорее всего будем завершать поиск, сразу же найдя нужную область после первой попытки. В очень маловероятном случае, когда нужная область находится в оставшихся
В обобщённом виде мы хотим решить следующую задачу:
При заданном множестве вероятностей найти двоичное дерево поиска этих вероятностей, минимизирующее ожидаемое время поиска.
К счастью, эта задача очень хорошо изучена и называется задачей оптимального двоичного дерева поиска. Существует множество алгоритмов для решения этой задачи; известно, что точное решение можно найти за время
Интересно что наилучший случай поведения таких оптимизированных двоичных деревьев поиска возникает, когда распределения вероятностей чрезвычайно перекошены, поскольку мы просто можем переместить в корень дерева узлы, содержащие подавляющую часть массы вероятностей, а наихудший случай — это когда распределение сбалансировано, потому что тогда дерево должно быть широким и неглубоким. Это противоположно поведению предыдущего алгоритма, в котором для симуляции шулерской кости использовалась честная!
В наилучшем случае у нас есть шулерская кость, в которой всегда выпадает одна грань (то есть она имеет вероятность 1, а все остальные грани имеют вероятность 0). Это чрезвычайное преувеличение нашего предыдущего примера, но в таком случае поиск будет всегда завершаться после первой попытки. В наихудшем случае все вероятности равны, и у нас получается стандартный поиск по BST. Мы приходим к следующему:
Алгоритм | Время инициализации | Время генерации | Занятая память | |||
---|---|---|---|---|---|---|
Наилучшее | Наихудшее | Наилучшее | Наихудшее | Наилучшая | Наихудшая | |
Шулерская кость из честной кости | ||||||
Шулерская кость из несимметричных монет | ||||||
Выбор колесом рулетки | ||||||
Оптимальный выбор колесом рулетки |
Бросание дротиков
Пока мы рассматривали два примитива, которые помогли нам построить алгоритмы для симуляции шулерской кости: честную кость и несимметричную монету. При использовании только честной кости мы приходим к (увы, непрактичному) алгоритму симуляции шулерской кости, а начав с несимметричными монетами, смогли изобрести быстрый алгоритм симуляции шулерской кости. Можно ли скомбинировать эти два подхода, чтобы создать алгоритм на основе и честных костей, и несимметричных монет? Оказывается, что да, и на самом деле получающийся алгоритм лучше обоих этих подходов.
До этого момента мы визуализировали интервал
Давайте начнём с возврата к нашему предыдущему примеру с вероятностями

Заметьте, что общая площадь этих прямоугольников равна
Теперь предположим, что мы нарисуем вокруг этих прямоугольников ограничивающий прямоугольник, ширина которого равна

Мы можем представить, что этот прямоугольник разделён на пять областей — четыре области соответствуют различным вероятностям, а одной областью обозначено неиспользуемое пространство. Взяв это разбиение, мы можем подумать об алгоритме симуляции случайных бросков кости как об игре в дартс. Предположим, что мы бросаем в эту цель (идеально равномерно распределённый) дротик. Если он попадает в неиспользуемое пространство, то мы достаём дротик и бросаем его снова, повторяя процесс, пока не попадём в один из прямоугольников. Поскольку чем больше вероятность, тем больше прямоугольник, то чем больше вероятность выбрасывания грани кости, тем выше вероятность попадания в её прямоугольник. На самом деле если мы поставим условие, что уже попали в какой-то прямоугольник, то у нас получится следующее:
Другими словами, когда мы наконец попадаем в какой-то прямоугольник своим равномерно распределённым дротиком, мы выбираем прямоугольник грани
Один из способов имитации бросков дротика в этот прямоугольник заключается в выборе двух равномерно распределённых значений в интервале
Первое наблюдение: мы показали, что ширину этих прямоугольников можно выбирать произвольно, потому что все они имеют равную ширину. Высоты, разумеется, зависят от вероятностей граней костей. Однако если мы равномерно отмасштабируем все высоты на какое-то положительное вещественное число
Теперь рассмотрим вероятность выбора любого отдельного прямоугольника, ограничившись условием, что мы точно попали в какой-то прямоугольник. Используя те же расчёты, мы получаем следующее:
То есть на самом деле вероятность выбора любого отдельного прямоугольника не изменяется, если мы масштабируем их линейно и равномерно.
Поскольку мы можем выбрать любой подходящий коэффициент масштабирования, то почему бы нам не отмасштабировать эти прямоугольники так, чтобы высота ограничивающего прямоугольника всегда равнялась 1? Так как высота ограничивающего прямоугольника определяется максимальным значением

Теперь мы готовы подумать о том, как мы можем бросить случайный дротик в прямоугольник и определить, во что он попал. Самое важное заключается в том, что мы можем разбить прямоугольник так, чтобы он не состоял из нескольких более мелких прямоугольников и пустого пространства странной формы. Вместо этого область разрезается на набор из

Заметьте, как образуется этот прямоугольник. Для каждой грани шулерской кости у нас есть один столбец шириной 1 и высотой 1, разделённая на два пространства — полупространство «да», соответствующее прямоугольнику этого размера, и полупространство «нет», соответствующее оставшейся части столбца.
Теперь давайте подумаем, как мы можем бросать дротик. Идеально равномерный дротик, брошенный в этот прямоугольник, будет иметь компоненты
Выбор случайной точки в этом интервале аналогичен бросанию честной кости и броску несимметричной монеты.
На самом деле этот результат можно воспринимать как гораздо более мощную возможность. Для симуляции шулерской кости мы строим набор несимметричных монет, по одной для каждой грани кости, а затем бросаем честную кость для определения того, какую монету нам бросать. На основании броска кости, если на соответствующей монете выпал орёл, то мы выбираем соответствующую грань, а если выпала решка, то бросаем кость снова и повторяем процесс.
Давайте повторим эти важные моменты. Во-первых, эти прямоугольники имеют следующие размеры — для каждой грани высота прямоугольника «да» задаётся как
Алгоритм: шулерская кость из честных костей/несимметричных монет
- Инициализация:
- Находим максимальное значение
; называем его
.
- Выделяем массив
длиной
, соответствующий высотам прямоугольников «да» в каждой строке.
- Для каждой вероятности
от
до
:
- Задаём
- Генерация:
- Пока не найдено значение:
- Бросаем честную n-гранную кость и получаем индекс
в интервале
.
- Бросаем несимметричную монету, на которой выпадает орёл с вероятностью
.
- Если на монете выпадает орёл, то возвращаем
.
Давайте проанализируем сложность этого алгоритма. На этапе инициализации поиск максимальной вероятности занимает время
Если это вероятность выбора какой-то монеты в любой одной итерации, то ожидаемое количество итераций, которое может произойти, задаётся величиной, обратной этой дроби, то есть
Интересно, что ожидаемое количество бросков монеты зависит исключительно от значения
Сравнение этого алгоритма с другими даёт нам следующую информацию:
Алгоритм | Время инициализации | Время генерации | Занятая память | |||
---|---|---|---|---|---|---|
Наилучшее | Наихудшее | Наилучшее | Наихудшее | Наилучшая | Наихудшая | |
Шулерская кость из честной кости | ||||||
Шулерская кость из несимметричных монет | ||||||
Выбор колесом рулетки | ||||||
Оптимальный выбор колесом рулетки | ||||||
Шулерская кость из честных костей/несимметричных монет |
В наилучшем случае этот алгоритм лучше показанного выше алгоритма двоичного поиска, потому что требует всего одного броска монеты. Однако наихудший случай экспоненциально хуже. Можно ли устранить этот наихудший случай?
Alias-метод
Предыдущая техника отлично показывает себя в наилучшем случае, генерируя случайный бросок кости с помощью всего лишь одной честной кости и одной монеты. Но ожидаемое наихудшее поведение намного хуже, и потенциально требует линейного количества бросков кости и монеты. Причина заключается в том, что в отличие от предыдущих техник, алгоритм может «промахиваться» и вынужден повторно выполнять итерации, пока не выберет решение. Графически это можно представить так, что мы бросаем дротики в цель, в которой есть большое количество пустого пространства, не привязанного ни к какому результату. Если бы был какой-то способ устранить всё это пустое пространство таким образом, чтобы каждая часть цели была покрыта прямоугольником, соответствующим какой-то грани шулерской кости, то мы бы просто бросили единственный дротик и считали бы результат.
Хитрый ход мы можем сделать, подстроив высоту прямоугольника таким образом, чтобы она соответствовала не наибольшей вероятности, а средней. Давайте вернёмся к нашему примеру. Вероятности равны

Теперь мы отмасштабируем все эти прямоугольники так, чтобы вероятность

На этом этапе давайте нарисуем поверх этого изображения прямоугольник

Как видите, это не совсем правильно, потому что прямоугольники для

Итак, у нас остались лишние полосы, но их уже не очень много. Мы могли бы полностью устранить лишнее, переместив выделяющиеся части из полос

Наконец, мы перетащим лишнее из второго столбца в первый, получив готовый прямоугольник:

Получившаяся схема обладает замечательными свойствами. Во-первых, общие площади прямоугольников. обозначающих каждую грань шулерской кости, не изменились по сравнению с исходной схемой; мы просто разрезали эти прямоугольники на части и переместили их. Это означает, что пока площади исходных прямоугольников пропорционально распределены согласно исходному распределению вероятностей, то общая площадь каждой грани кости остаётся той же. Во-вторых, заметьте, что в этом новом прямоугольнике нет свободного места, то есть при любом броске дротика мы гарантировано куда-нибудь попадём и получим окончательный ответ, а не пустое пространство, требующее ещё одного броска. Это значит, что единственного броска дротика будет достаточно для генерирования нашего случайного значения. Последнее и самое важное — заметьте, что в каждом столбце содержится не более двух разных прямоугольников. Это значит, что можно воспользоваться предыдущим изобретением — мы бросаем кость, чтобы выбрать несимметричную монету, а затем бросаем монету. Отличие на этот раз заключается в том, что означает бросок монеты. Выпавший орёл означает, что мы выбираем одну грань кости, а выпавшая решка теперь означает, что мы должны выбрать какую-то другую грань кости (а не бросать кость заново).
На высоком уровне alias-метод работает следующим образом. Во-первых, мы создаём прямоугольники, обозначающие каждую из разных вероятностей граней кости. Далее мы разрезаем эти прямоугольники на части и перераспределяем их таким образом, чтобы они полностью заполняли прямоугольную цель, чтобы при этом каждый столбец имел постоянную ширину и содержал прямоугольники максимум двух разных сторон шулерской кости. Наконец, мы симулируем броски кости, случайным образом метая дротик в цель, что можно реализовать сочетанием честной кости и несимметричных монет.
Но как нам узнать, возможно ли вообще разрезать исходные прямоугольники таким образом, чтобы каждый столбец содержал не более двух разных вероятностей? Это не кажется очевидным, но удивительно, что этого действительно можно добиться всегда. Более того, мы не только можем разрезать прямоугольники на части таким образом, что каждый столбец содержит не более двух разных прямоугольников, но и сделать так, чтобы один из прямоугольников в каждом столбце был прямоугольником, изначально расположенным в этом столбце. Можно заметить, что в показанной выше схеме прямоугольников мы всегда вырезали часть прямоугольника и перемещали его в другой столбец и никогда не убирали прямоугольник их исходного столбца полностью. Это означает, что каждый столбец в готовой схеме будет состоять из какого-то прямоугольника, соответствующего вероятности, изначально назначенной этому столбцу, плюс (не обязательно) из второго прямоугольника, перетащенного из какого-то другого столбца. Такой второй прямоугольник часто называют псевдонимом (alias) столбца, потому что оставшаяся вероятность столбца используется в качестве «псевдонима» какой-то другой вероятности. Из-за использования термина «alias» этот метод получил название «alias-метод».
Прежде чем мы перейдём к доказательству того, что можно всегда распределить вероятности таким образом, давайте вкратце опишем суть работы алгоритма. Поскольку каждый столбец готовой схемы всегда содержит какую-то часть исходного прямоугольника из этого столбца (потенциально с нулевой высотой!), то для хранения (потенциально) двух различных прямоугольников, занимающих столбец, реализации alias-метода обычно хранят две разные таблицы: таблицу вероятностей

Доказательство существования таблиц Alias
Теперь нам нужно формально доказать, что возможность создания таблиц
- Создать прямоугольники
для каждой из вероятностей
,
- разрезать их горизонтально на меньшие части и
- распределить их на
столбцов
- так, что высота каждого столбца равна
,
- ни в одном столбце нет больше двух прямоугольников и
- прямоугольник, соответствующий грани
, расположен в столбце
.
- так, что высота каждого столбца равна
Прежде чем переходить к доказательствам того, что это возможно всегда, давайте рассмотрим пример. Допустим, у нас есть четыре вероятности

Заметьте, что два из четырёх прямоугольников, которые нам нужно распределить (

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

Теперь посмотрим, на что похожа наша схема. У нас есть три прямоугольника, общая площадь которых равна

Теперь нам нужно заполнить столбец до конца, поэтому мы возьмём какую-нибудь вероятность, которая не меньше 1, и используем её для дополнения оставшейся части столбца. Здесь у нас есть только один вариант (использовать

Теперь всё свелось к двум прямоугольникам, общая площадь которых равна двум. Мы повторяем этот процесс, найдя какой-нибудь прямоугольник, высота которого не больше 1 (здесь это

Теперь мы найдём какой-нибудь прямоугольник высотой не менее

Теперь у нас остался только один прямоугольник, и его площадь равна 1. Поэтому мы можем просто завершить создание, поместив этот прямоугольник в его собственный столбец:

И вуаля! Мы заполнили таблицу.
Заметим, что в основе конструирования этой таблицы лежит общий паттерн:
- Находим какой-то прямоугольник, высота которого не больше 1, и помещаем его в его собственный столбец, задавая в таблице
высоту этого прямоугольника.
- Находим какой-то прямоугольник, высота которого не меньше 1, используем его для заполнения столбца, задавая в таблице
соответствие грани кости, которую представляет этот прямоугольник.
Можно ли доказать, что такое построение в общем случае всегда возможно? То есть что мы не «застрянем», распределяя вероятности таким образом? К счастью, доказательство есть. Можно объяснить это так: мы отмасштабировали все вероятности так, чтобы среднее значение новых вероятностей равнялось 1 (поскольку изначально оно было равно
Мы можем формализировать это рассуждение в следующей теореме:
Теорема: при наличиипрямоугольников единичной ширины с высотами
,
, ...,
, такие, что
, всегда есть способ разделения прямоугольников и распределения их в
столбцов, каждый из которых имеет высоту 1, таким образом, что каждый столбец будет содержать не более двух разных прямоугольников, и в
-том столбце будет содержаться хотя бы одна часть
-того прямоугольника.
Доказательство: по индукции. В базовом случае, если, то у нас есть только один прямоугольник и его высота должна быть равна 1. Поэтому мы можем занести его в
-ой столбец. Следовательно, каждый столбец имеет высоту 1, содержит не более двух прямоугольников, а в
-том столбце содержится не менее одной части
-того прямоугольника.
Для этапа индукции предположим, что для какого-то натурального числатеорема справедлива и рассмотрим любые
прямоугольников шириной
и высотами
,
, ...,
, такие, что
. Сначала мы утверждаем, что существует некая высота
, такая, что
, и какая-то другая высота
(такая, что
), такая, что
. Чтобы убедиться в этом, пойдём от обратного и допустим, что нет такой
с
; это будет значит, что
для всех натуральных чисел
в интервале
. Но тогда у нас получится, что
, что очевидно невозможно. Следовательно, существует какой-то индекс
, такой, что
. Теперь снова пойдём от обратного и допустим, что нет никакой другой высоты
(с
), такой, что
. Тогда у нас получится, что каждая другая
, то есть (по аналогичной логике)
и мы пришли к противоречию. Следовательно,
и
.
Теперь рассмотрим следующую конструкцию. Поместимв столбец
и заполним оставшееся пространство
в
-том столбце частью прямоугольника
(такое пространство должно существовать, поскольку
и
). Так мы полностью заполним столбец. Теперь у нас остался набор из
разных частей прямоугольников, чья общая сумма равна
, поскольку мы убрали из прямоугольников общую площадь
, а изначальная общая сумма была равна
. Более того, мы полностью заполнили столбец
, то есть нам больше не понадобится размещать там другие части прямоугольников. Следовательно, по индуктивной гипотезе, мы можем распределить оставшиеся
в
столбцов так, чтобы это удовлетворяло приведённым выше условиям. В сочетании с тем фактом, что мы теперь заполнили столбец
, это означает, что у нас есть способ заполнения всех столбцов, чтобы это удовлетворяло ограничениям. На этом индукция завершена.
Это доказательство возможности построения, гласящее, что мы не только всегда можем построить таблицу alias, но и то, что приведённый выше алгоритм поиска высоты не больше единицы и соединение её в пару с прямоугольником высотой не меньше единицы всегда оказывается успешным. Из этого мы можем начать выводить всё более и более сложные алгоритмы вычисления таблиц alias.
Генерирование таблиц Alias
Используя сказанное выше, мы можем получить достаточно хороший алгоритм для симуляции бросков шулерской кости с помощью alias-метода. Инициализация заключается в повторяющемся сканировании входящих вероятностей для нахождения значения не больше 1 и значения не меньше 1, с последующим объединением их в пару для заполнения столбца:
Алгоритм: наивный Alias-метод
- Инициализация:
- Умножаем каждую вероятность
на
.
- Создаём массивы
и
, каждый из которых имеет размер
.
- For
:
- Находим вероятность
, удовлетворяющую условию
.
- Находим вероятность
(с
), удовлетворяющую условию
- Задаём
.
- Задаём
.
- Удаляем
из списка изначальных вероятностей.
- Задаём
.
- Пусть
будет последней оставшейся вероятностью, которая должна иметь высоту 1.
- Задаём
.
- Генерация:
- Генерируем бросок честной кости из
-гранной кости; называем грань
.
- Бросаем несимметричную монету, выпадающую орлом с вероятностью
.
- Если на монете выпадает орёл, возвращаем
.
- В противном случае возвращаем
.
Этап генерации этого алгоритма точной такой же, как описанный выше метод, и выполняется за время
Алгоритм | Время инициализации | Время генерации | Занятая память | |||
---|---|---|---|---|---|---|
Наилучшее | Наихудшее | Наилучшее | Наихудшее | Наилучшая | Наихудшая | |
Шулерская кость из честной кости | ||||||
Шулерская кость из несимметричных монет | ||||||
Выбор колесом рулетки | ||||||
Оптимальный выбор колесом рулетки | ||||||
Шулерская кость из честных костей/несимметричных монет | ||||||
Наивный Alias-метод |
По сравнению с другими эффективными методами симуляции этот наивный alias-метод имеет высокие затраты на инициализацию, но затем может чрезвычайно эффективно симулировать броски кости. Если бы мы могли каким-то образом снизить затраты на инициализацию (допустим, до
Один простой способ снижения затрат на инициализацию заключается в использовании более подходящей структуры данных для сохранения высот в процессе выполнения. В наивной версии мы используем для хранения всех вероятностей неотсортированный массив, то есть для нахождения двух нужных вероятностей требуется время
Алгоритм: Alias-метод
- Инициализация:
- Создаём массивы
и
, каждый размером
.
- Создаём сбалансированное двоичное дерево поиска
.
- Вставляем
в
для каждой вероятности
.
- For
:
- Находим и удаляем наименьшее значение в
; назовём его
.
- Находим и удаляем наибольшее значение в
; назовём его
.
- Задаём
.
- Задаём
.
- Задаём
.
- Добавляем
к
.
- Пусть
будет последней оставшейся вероятностью, которая должна иметь вес 1.
- Задаём
.
- Генерация:
- Генерируем бросок честной кости из
-гранной кости; называем грань
.
- Бросаем несимметричную монету, на которой выпадает орёл с вероятностью
.
- Если на монете выпадает орёл, возвращаем
.
- В противном случае возвращаем
.
Теперь инициализация нашего алгоритма происходит гораздо быстрее. Создание
Алгоритм | Время инициализации | Время генерации | Занятая память | |||
---|---|---|---|---|---|---|
Наилучшее | Наихудшее | Наилучшее | Наихудшее | Наилучшая | Наихудшая | |
Шулерская кость из честной кости | ||||||
Шулерская кость из несимметричных монет | ||||||
Выбор колесом рулетки | ||||||
Оптимальный выбор колесом рулетки | ||||||
Шулерская кость из честных костей/несимметричных монет | ||||||
Наивный Alias-метод | ||||||
Alias-метод |
Однако существует алгоритм, работающий ещё быстрее. Он на удивление прост, и, вероятно, является наиболее чистым из алгоритмов реализации alias-метода. Этот алгоритм впервые описан в статье «A Linear Algorithm For Generating Random Numbers With a Given Distribution» Майкла Воуза, и стал стандартным алгоритмом реализации alias-метода.
В основе алгоритма Воуза лежит использование двух рабочих списков: один содержит элементы с высотой меньше 1, другой содержит элементы с высотой не меньше 1. Алгоритм постоянно соединяет парами первые элементы каждого рабочего списка. На каждой итерации мы потребляем элемент из списка «малых» значений, и потенциально перемещаем остаток элемента из списка «больших» значений в список «малых». В этом алгоритме используется несколько инвариантов:
- Все элементы «малого» рабочего списка меньше 1.
- Все элементы «большого» рабочего списка не меньше 1.
- Сумма элементов в рабочих списка всегда равна общему количеству элементов.
Для упрощения каждый список хранит не саму вероятность, а некий указатель на список исходных вероятностей, сообщающий, на какую грань шулерской кости выполняется ссылка. Имея эти инварианты, мы получаем следующий алгоритм:
Алгоритм: (нестабильный) Alias-метод Воуза
Предупреждение: этот алгоритм подвержен численным неточностям. Более численно стабильный алгоритм представлен ниже.
- Инициализация:
- Создаём массивы
и
, каждый размером
.
- Создаём два рабочих списка,
и
.
- Умножаем каждую вероятность на
.
- Для каждой отмасштабированной вероятности
:
- Если
, добавляем
к
.
- В противном случае (
) добавляем
к
.
- Пока список
не пуст:
- Удаляем первый элемент из
; называем его
.
- Удаляем первый элемент из
; называем его
.
- Задаём
.
- Задаём
.
- Задаём
.
- Если
, добавляем
в
.
- В противном случае ($p_g \ge 1$) добавляем
в
.
- Пока список
не пуст:
- Удаляем первый элемент из
; называем его
.
- Задаём
.
- Генерация:
- Генерируем бросок честной кости из
-гранной кости; называем грань
.
- Бросаем несимметричную монету, на которой выпадает орёл с вероятностью
.
- Если на монете выпадает орёл, возвращаем
.
- В противном случае возвращаем
.
С учётом трёх представленных выше инвариантов первая часть алгоритма (всё, за исключением последнего цикла) должна быть достаточно понятна: мы постоянно соединяем в пары какой-то малый элемент из
В этом алгоритме тип рабочего списка не важен. В первоначальной статье Воуза в качестве рабочего списка используются стеки, потому что их можно эффективно реализовать с помощью массивов, но если мы захотим, то можем использовать очередь. Однако для простоты мы воспользуемся стеком.
Прежде чем переходить к анализу алгоритма, давайте разберём его работу на примере. Рассмотрим пример с использованием семи вероятностей

Теперь мы разместим вершину стека

Теперь мы используем вершину стека

Затем мы повторяем этот процесс. Перемещаем прямоугольник из вершины стека

И ещё раз:

Когда мы повторим этот процесс в следующий раз, то обнаружим, что хотя мы и можем использовать бирюзовый блок для заполнения пустого пространства в таблице, но в результате этого высота бирюзового блока станет меньше единицы. Следовательно, нам нужно будет переместить бирюзовый блок на вершину стека малых значений:

Теперь при обработке списка

Затем мы обрабатываем стек

И наконец, поскольку стек

И у нас получилась хорошо сформированная таблица alias для этих вероятностей.
Практичная версия алгоритма Воуза
К сожалению, приведённый выше алгоритм в таком виде численно нестабилен. Он вполне точен на идеализированной машине, способной выполнять вычисления с вещественными числами произвольной точности, но если попробовать выполнить его с помощью IEEE-754 double, то он в результате может оказаться совершенно ошибочным. Вот два источника неточностей, с которыми нам нужно разобраться, прежде чем двигаться дальше:
- Вычисление того, принадлежит ли вероятность к рабочему списку
или к
, может быть неточным. В частности, может возникнуть ситуация, в которой масштабирование вероятностей на коэффициент
приводит к том, что вероятности, равные
, в результате оказываются немного меньше
(то есть попадают в список
, а не в
).
- Вычисление, вычитающее соответствующую массу вероятностей из большей вероятности, численно нестабильно и может вносить значительные ошибки округления. Это может приводить к том, что вероятность, которая должна быть в списке
, оказывается в списке
.
Сочетание этих двух факторов означает, что в результате алгоритм может случайно поместить все вероятности в список
К счастью, исправить это оказывается не так сложно. Мы изменим внутренний цикл алгоритма так, чтобы он завершался, когда любой из списков оказывается пустым, чтобы нам случайно не пришлось искать несуществующие элементы в списке
Алгоритм: Alias-метод Воуза
- Инициализация:
- Создаём массивы
и
, каждый размером
.
- Создаём два рабочих списка,
и
.
- Умножаем каждую вероятность на
.
- Для каждой отмасштабированной вероятности
:
- Если
, добавляем
в
.
- В противном случае (
) добавляем
в
.
- Пока
и
не пусты: (
может быть сначала опустошён)
- Удаляем первый элемент из
; называем его
.
- Удаляем первый элемент из
; называем его
.
- Задаём
.
- Задаём
.
- Задаём
. (Это более численно стабильный вариант.)
- Если
, добавляем
в
.
- В противном случае (
) добавляем
в
.
- Пока
не пуст:
- Удаляем первый элемент из
; называем его
.
- Задаём
.
- Пока
не пуст: Это возможно только из-за численной нестабильности.
- Удаляем первый элемент из
; называем его
.
- Задаём
.
- Генерация:
- Генерируем бросок честной кости из
-гранной кости; называем грань
.
- Бросаем несимметричную монету, на которой выпадает орёл с вероятностью
.
- Если на монете выпадает орёл, возвращаем
.
- В противном случае возвращаем
.
Единственное, что нам осталось — проанализировать сложность алгоритма. Создание рабочих списков в общем занимает время
Алгоритм | Время инициализации | Время генерации | Занятая память | |||
---|---|---|---|---|---|---|
Наилучшее | Наихудшее | Наилучшее | Наихудшее | Наилучшая | Наихудшая | |
Шулерская кость из честной кости | ||||||
Шулерская кость из несимметричных монет | ||||||
Выбор колесом рулетки | ||||||
Оптимальный выбор колесом рулетки | ||||||
Шулерская кость из честных костей/несимметричных монет | ||||||
Наивный Alias-метод | ||||||
Alias-метод | ||||||
Alias-метод Воуза |
Мысли в заключение
Ого! Мы проделали большую работу! Исследовали различные методы симуляции шулерской кости, начиная с самого простого набора техник и заканчивая чрезвычайно быстрыми и эффективными алгоритмами. Каждый метод демонстрирует собственный набор техник, и я считаю окончательную версию (alias-метод Воуза) одним из самых интересных и элегантных алгоритмов, с которыми мне когда-либо приходилось встречаться.
Если вам интересно посмотреть код alias-метод Воуза, в том числе краткое описание сложностей, возникающих на практике из-за численной неточности, то у меня есть реализация alias-метода на Java, выложенная в архиве интересного кода.
Надеюсь, статья была вам полезна!