Алгоритм Зиккурат

Материал из Википедии — свободной энциклопедии
Перейти к навигации Перейти к поиску

Алгоритм «Зиккурат» (англ. Ziggurat Algorithm, Ziggurat Method) — это алгоритм выборки псевдослучайных чисел. Будучи представителем класса алгоритмов выборки с отклонением, он в работе своей опирается на источник равномерно распределённых случайных чисел — обыкновенно это генератор псевдослучайных чисел, либо же предварительно вычисленная таблица. Алгоритм используется для генерации значений на основе монотонно убывающего вероятностного распределения. Также может быть применён по отношению к симметричному унимодальному распределению, такому как нормальное, с помощью выбора значений из одной его половины, а затем, при необходимости, перехода к симметричному значению с помощью операции арифметического отрицания. Одним из авторов алгоритма, разработанного в 1960-е, является Джордж Марсалья.

В простейшем случае для вычисления значения, возвращаемого алгоритмом, требуется только генерация одного числа с плавающей точкой и одного случайного табличного индекса, за которой следует один табличный поиск, одна операция умножения и одно сравнение. Иногда (в гораздо меньшем количестве случаев) требуются более сложные вычисления. Тем не менее, данный алгоритм гораздо быстрее с вычислительной точки зрения, чем два наиболее часто используемых метода генерации нормально распределённых случайных чисел: полярного метода Марсальи и преобразования Бокса — Мюллера, где требуется вычисление по меньшей мере одного логарифма и одного квадратного корня для каждой пары генерируемых значений. Однако, так как алгоритм «Зиккурат» более сложен в реализации, наиболее часто он используется в случаях, где требуется большое количество случайных чисел.

Сам термин «алгоритм „Зиккурат“» (Ziggurat Algorithm) фигурирует в совместной работе Марсальи и Ваи Ван Тсанга 2000 года и назван так потому, что концептуально основан на покрытии распределения вероятностей прямоугольными сегментами, сложенными друг на друге в порядке уменьшения размера (если рассматривать снизу вверх), что приводит к появлению фигуры, напоминающей зиккурат.

Визуализация алгоритма
Пример работы с одной частью нормального распределения. Розовые точки — равномерно распределённые случайные числа. Заданная функция распределения делится на области равных площадей . Уровень выбирается случайным образом (источник — равномерно распределённые числа слева). Затем случайное значение из верхнего источника домножается на ширину выбранного уровня, и результирующий проверяется на принадлежность одной из 3 возможных областей: 1) (левая, черная область) выборка точно под кривой, тут же возвращается, 2) (заштрихованная область) значение может как принадлежать кривой, так и нет. В этом случае генерируется случайный , принадлежащий выбранному уровню и сравнивается с . Если меньше, то точка под кривой, и возвращается. Если же нет, 3) выбранная точка отклоняется алгоритмом и всё сначала.

Теоретическая база

[править | править код]

Алгоритм «Зиккурат» — это алгоритм выборки с отклонением. Он случайным образом генерирует точку, незначительно отклоненную от нужного распределения, а затем проверяет попала ли сгенерированная точка точно внутрь такового. Если нет, алгоритм пробует заново. Если точка лежит под кривой вероятностной функции плотности, то ее x-координата и будет искомым случайным числом с нужным распределением.

Распределение, из которого алгоритм производит выборку, состоит из областей равной площади; прямоугольник покрывает основную часть нужного распределения и располагается «пирамидкой» на не-прямоугольном основании, которое включает в себя остаточную часть или «хвост» распределения.

У заданной монотонно убывающей вероятностной функции плотности , определенной для всех , основание зиккурата определяется как все точки внутри распределения и ниже некоторого . Оно состоит из прямоугольной части от до , и (обычно бесконечного) остатка (хвоста) распределения, где ).

Этот уровень (назовем его уровень 0) имеет площадь . Добавим на его вершину новый прямоугольный уровень ширины и высоты , так что у него тоже площадь будет равна . Вершина этого уровня находится на высоте , и пересекает функцию плотности в точке , где . Этот уровень включает в себя все точки функции плотности между и , но (в отличие от базового уровня) также включает прочие точки, такие, как , которые не принадлежат нужному распределению.

Все последующие уровни накладываются друг на друга аналогичным образом. Для использования предварительно вычисленной таблицы размера ( используется очень часто), следует выбрать так, что , таким образом верхний прямоугольный уровень с номером достигнет пика распределения в точности в точке .

Уровень с номером в высоту занимает место от до , и по ширине может быть разделен на две области: часть от до (обыкновенно бо́льшую), которая целиком содержится внутри заданного распределения, и часть от до (меньшую), которая содержится внутри лишь частично.

Забывая ненадолго о вопросе особого случая с уровнем 0, и имея равномерно распределенные числа и , алгоритм может быть описан следующим образом:

  1. Выбрать случайным образом уровень .
  2. Положить .
  3. Если , вернуть .
  4. Положить .
  5. Вычислить . Если , вернуть .
  6. В противном случае выбрать новые случайные числа и вернуться к шагу 1.

Шаг 1 является случайной выборкой уровня. Шаг 3 проверяет, лежит ли координата чётко внутри заданной функции плотности даже без какой-либо информации о координате . Если не лежит, шаг 4 производит вычисление координаты , и шаг 5 производит проверку на попадание внутрь нужной области.

Если число уровней достаточно велико и они имеют малую высоту, то та самая «зона риска», проверка которой производится после шага 3, очень мала, и алгоритм останавливается на шаге 3 существенную часть времени. Стоит обратить внимание на то, что верхний уровень , однако, этот тест всегда проваливает, так как .

Уровень 0 также может быть разделен на центральную и граничную области, но граничная будет содержать бесконечный остаток функции. Для использования того же алгоритма для проверки принадлежности точки центральной области, стоит сгенерировать фиктивную . Точки с координатой будут обрабатываться просто, а для того редкого случая, когда был выбран уровень 0 и , придется использовать особый резервный алгоритм для случайной выборки точки из «хвоста» функции. Поскольку такой запасной алгоритм будет задействован чрезвычайно редко (редкость относительна и зависит от разбиения на уровни), то его скорость не окажет существенного влияния на производительность в целом.

Таким образом, полный алгоритм «Зиккурат» для несимметричного распределения выглядит следующим образом:

  1. Выбрать случайный уровень .
  2. Положить .
  3. Если , вернуть .
  4. Если , сгенерировать точку из «хвоста» с использованием запасного алгоритма.
  5. Положить .
  6. Вычислить . Если , вернуть .
  7. В противном случае выбрать новые случайные числа и вернуться к шагу 1.

Для симметричного распределения результат, конечно же, может просто становиться противоположного знака в 50 % случаев. Часто может быть удобно сгенерировать и на шаге 3 проверить .

Запасные алгоритмы для хвостовой части функции

[править | править код]

Так как алгоритм «Зиккурат» очень быстро генерирует только лишь большую часть значений и требует наличия запасного алгоритма в случаях , дела обстоят сложнее непосредственной реализации из 6 шагов. Запасной алгоритм зависит от заданного распределения.

В случае показательного распределения, хвостовая часть имеет вид тела распределения. Один из способов — вернуться к самому элементарному алгоритму и положить . Другой способ состоит в рекурсивном вызове алгоритма «Зиккурат» и прибавлении к результату.

В случае нормального распределения, Марсалья предлагает компактный алгоритм:

  1. Положить .
  2. Положить .
  3. Если , вернуть .
  4. В противном случае вернуться к шагу 1.

Так как для таблиц более-менее типичных размеров, тест на шаге 3 почти всегда успешен.

Оптимизации

[править | править код]

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

  • В алгоритме ничего не зависит от того, нормализована ли вероятностная функция распределения (значение интеграла равняется 1), так что удаление нормализующей константы может ускорить вычисление .
  • Большинство генераторов равномерно распределенных случайных чисел основаны на генераторах случайных целых чисел, которые возвращают целое число из диапазона . Таблица, содержащая позволит использовать такие числа напрямую в качестве .
  • В случае работы с симметричными распределениями при использовании симметричной как было описано выше, случайное целое число может быть интерпретировано как знаковое число в диапазоне , и может использоваться масштабирующий коэффициент .
  • Вместо сравнения с на шаге 3, возможно вычислить заранее и сравнивать с этим значением напрямую. Если  — это генератор целых случайных чисел, значения могут быть заранее домножены на (или , в соответствующем случае) так что будет проводиться целочисленное сравнение.
  • С учетом двух изменений выше, таблица исходных значений более не нужна и может быть удалена.
  • В случае генерации чисел с плавающей точкой одинарной точности согласно стандарту IEEE 754, где используется 24-битная мантисса (включая неявно заданную 1), наименее значимые биты 32-битного целого случайного числа не используются. Эти биты могут быть использованы при выборе уровня. (тут[1] подробно описана суть вопроса).

Генерация таблиц

[править | править код]

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

Можно находить и . Повторить для всех уровней зиккурата. В конце должно получиться .

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

Поиск и

[править | править код]

Если имеется начальное значение (вычисленное если и не точно, то приближенно), останется лишь вычислить площадь хвостовой части функции, для которой выполнено . Вычислить можно методами численного интегрирования.

Далее, из можно найти , из площади хвостовой части найдется площадь базового уровня: .

Затем вычисляется серия и как показано выше. Если для любых , тогда начальное значение было слишком малым, что привело к большой площади . Если , тогда начальное значение было слишком большим.

Учитывая сказанное, можно использовать численное решение уравнений (например, метод бисекции) для нахождения значения , при котором значение так близко к как только возможно. В качестве альтернативы можно рассматривать и находить значения площади верхнего уровня, , настолько близкие к нужному значению как только возможно.

Примечания

[править | править код]
  1. Jurgen A. Doornik. "An Improved Ziggurat Method to Generate Normal Random Samples" (англ.) // Nuffield College, Oxford. — 2005. Архивировано 7 марта 2016 года.

Литература

[править | править код]
  • George Marsaglia, Wai Wan Tsang. The Ziggurat Method for Generating Random Variables // Journal of Statistical Software. — 2000. — 7 с. — URL: сайт
  • Jurgen A. Doornik. An Improved Ziggurat Method to Generate Normal Random Samples. — Nuffield College, Oxford: 2005. — 9 с. — URL: работа
  • David B. Thomas, Philip H.W. Leong, Wayne Luk, John D. Villasenor. Gaussian Random Number Generators // ACM Computing Surveys. — 2007. — 38 с. — URL: работа
  • Boaz Nadler. Design Flaws in the Implementation of the Ziggurat and Monty Python methods (and some remarks on Matlab randn) // The Journal of Business. — 2006. — 16 с. — URL: работа
  • Edrees, Hassan M.; Cheung, Brian; Sandora, McCullen; Nummey, David; Stefan, Deian. Hardware-Optimized Ziggurat Algorithm for High-Speed Gaussian Random Number Generators // 2009 International Conference on Engineering of Reconfigurable Systems & Algorithms. Las Vegas. — URL: сайт
  • Marsaglia, George. Generating a Variable from the Tail of the Normal Distribution // Technometrics. — 1964. — Т. 6, № 1. — С 101—102. — URL: сайт