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

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

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

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

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

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

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

У заданной монотонно убывающей вероятностной функции плотности ${\displaystyle f(x)}$, определенной для всех ${\displaystyle x\geqslant 0}$, основание зиккурата определяется как все точки внутри распределения и ниже некоторого ${\displaystyle y_{1}=f(x_{1})}$. Оно состоит из прямоугольной части от ${\displaystyle (0,0)}$ до ${\displaystyle (x_{1},y_{1})}$, и (обычно бесконечного) остатка (хвоста) распределения, где ${\displaystyle x>x_{1}}$ (и ${\displaystyle y<y_{1}}$).

Этот уровень (назовем его уровень 0) имеет площадь ${\displaystyle A}$. Добавим на его вершину новый прямоугольный уровень ширины ${\displaystyle x_{1}}$ и высоты ${\displaystyle A/x_{1}}$, так что у него тоже площадь будет равна ${\displaystyle A}$. Вершина этого уровня находится на высоте ${\displaystyle y_{2}=y_{1}+A/x_{1}}$, и пересекает функцию плотности в точке ${\displaystyle (x_{2},y_{2})}$, где ${\displaystyle y_{2}=f(x_{2})}$. Этот уровень включает в себя все точки функции плотности между ${\displaystyle y_{1}}$ и ${\displaystyle y_{2}}$, но (в отличие от базового уровня) также включает прочие точки, такие, как ${\displaystyle (x_{1},y_{2})}$, которые не принадлежат нужному распределению.

Все последующие уровни накладываются друг на друга аналогичным образом. Для использования предварительно вычисленной таблицы размера ${\displaystyle n}$ (${\displaystyle n=256}$ используется очень часто), следует выбрать ${\displaystyle x_{1}}$ так, что ${\displaystyle x_{n}=0}$, таким образом верхний прямоугольный уровень с номером ${\displaystyle n-1}$ достигнет пика распределения в точности в точке ${\displaystyle (0,f(0))}$.

Уровень с номером ${\displaystyle i}$ в высоту занимает место от ${\displaystyle y_{i}}$ до ${\displaystyle y_{i+1}}$, и по ширине может быть разделен на две области: часть от ${\displaystyle 0}$ до ${\displaystyle x_{i+1}}$ (обыкновенно бо́льшую), которая целиком содержится внутри заданного распределения, и часть от ${\displaystyle x_{i+1}}$ до ${\displaystyle x_{i}}$ (меньшую), которая содержится внутри лишь частично.

Забывая ненадолго о вопросе особого случая с уровнем 0, и имея равномерно распределенные числа ${\displaystyle U_{0}}$ и ${\displaystyle U_{1}}$ ${\displaystyle \in [0,1)}$, алгоритм может быть описан следующим образом:

1. Выбрать случайным образом уровень ${\displaystyle 0\leqslant i<n}$.
2. Положить ${\displaystyle x=U_{0}x_{i}}$.
3. Если ${\displaystyle x<x_{i+1}}$, вернуть ${\displaystyle x}$.
4. Положить ${\displaystyle y=y_{i}+U_{1}(y_{i+1}-y_{i})}$.
5. Вычислить ${\displaystyle f(x)}$. Если ${\displaystyle y<f(x)}$, вернуть ${\displaystyle x}$.
6. В противном случае выбрать новые случайные числа и вернуться к шагу 1.

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

Если число ${\displaystyle n}$ уровней достаточно велико и они имеют малую высоту, то та самая «зона риска», проверка которой производится после шага 3, очень мала, и алгоритм останавливается на шаге 3 существенную часть времени. Стоит обратить внимание на то, что верхний уровень ${\displaystyle n-1}$, однако, этот тест всегда проваливает, так как ${\displaystyle x_{n}=0}$.

Уровень 0 также может быть разделен на центральную и граничную области, но граничная будет содержать бесконечный остаток функции. Для использования того же алгоритма для проверки принадлежности точки центральной области, стоит сгенерировать фиктивную ${\displaystyle x_{0}=A/y_{1}}$. Точки с координатой ${\displaystyle x<x_{1}}$ будут обрабатываться просто, а для того редкого случая, когда был выбран уровень 0 и ${\displaystyle x\geqslant x_{1}}$, придется использовать особый резервный алгоритм для случайной выборки точки из «хвоста» функции. Поскольку такой запасной алгоритм будет задействован чрезвычайно редко (редкость относительна и зависит от разбиения на уровни), то его скорость не окажет существенного влияния на производительность в целом.

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

1. Выбрать случайный уровень ${\displaystyle 0\leqslant i<n}$.
2. Положить ${\displaystyle x=U_{0}x_{i}}$.
3. Если ${\displaystyle x<x_{i+1}}$, вернуть ${\displaystyle x}$.
4. Если ${\displaystyle i=0}$, сгенерировать точку из «хвоста» с использованием запасного алгоритма.
5. Положить ${\displaystyle y=y_{i}+U_{1}(y_{i+1}-y_{i})}$.
6. Вычислить ${\displaystyle f(x)}$. Если ${\displaystyle y<f(x)}$, вернуть ${\displaystyle x}$.
7. В противном случае выбрать новые случайные числа и вернуться к шагу 1.

Для симметричного распределения результат, конечно же, может просто становиться противоположного знака в 50 % случаев. Часто может быть удобно сгенерировать ${\displaystyle U_{0}\in (-1,1)}$ и на шаге 3 проверить ${\displaystyle |x|<x_{i+1}}$.

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

Так как алгоритм «Зиккурат» очень быстро генерирует только лишь большую часть значений и требует наличия запасного алгоритма в случаях ${\displaystyle x>x_{1}}$, дела обстоят сложнее непосредственной реализации из 6 шагов. Запасной алгоритм зависит от заданного распределения.

В случае показательного распределения, хвостовая часть имеет вид тела распределения. Один из способов — вернуться к самому элементарному алгоритму ${\displaystyle E=-\ln(U_{1})}$ и положить ${\displaystyle x=x_{1}-\ln(U_{1})}$. Другой способ состоит в рекурсивном вызове алгоритма «Зиккурат» и прибавлении ${\displaystyle x_{1}}$ к результату.

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

1. Положить ${\displaystyle x=-\ln(U_{1})/x_{1}}$.
2. Положить ${\displaystyle y=-\ln(U_{2})}$.
3. Если ${\displaystyle 2y>x^{2}}$, вернуть ${\displaystyle x+x_{1}}$.
4. В противном случае вернуться к шагу 1.

Так как ${\displaystyle x_{1}\approx 3.5}$ для таблиц более-менее типичных размеров, тест на шаге 3 почти всегда успешен.

Оптимизации

Алгоритм может быть выполнен эффективно с использованием заранее вычисленных таблиц ${\displaystyle x_{i}}$ и ${\displaystyle y_{i}=f(x_{i})}$, но есть несколько модификаций, чтобы ускорить его еще больше:

• В алгоритме ничего не зависит от того, нормализована ли вероятностная функция распределения (значение интеграла равняется 1), так что удаление нормализующей константы может ускорить вычисление ${\displaystyle f(x)}$.
• Большинство генераторов равномерно распределенных случайных чисел основаны на генераторах случайных целых чисел, которые возвращают целое число из диапазона ${\displaystyle [0,2^{32}-1]}$. Таблица, содержащая ${\displaystyle 2^{-32}x_{i}}$ позволит использовать такие числа напрямую в качестве ${\displaystyle U_{0}}$.
• В случае работы с симметричными распределениями при использовании симметричной ${\displaystyle U_{0}}$ как было описано выше, случайное целое число может быть интерпретировано как знаковое число в диапазоне ${\displaystyle [-2^{31},2^{31}-1]}$, и может использоваться масштабирующий коэффициент ${\displaystyle 2^{-31}}$.
• Вместо сравнения ${\displaystyle U_{0}x_{i}}$ с ${\displaystyle x_{i+1}}$ на шаге 3, возможно вычислить заранее ${\displaystyle x_{i+1}/x_{i}}$ и сравнивать ${\displaystyle U_{0}}$ с этим значением напрямую. Если ${\displaystyle U_{0}}$ — это генератор целых случайных чисел, значения могут быть заранее домножены на ${\displaystyle 2^{32}}$ (или ${\displaystyle 2^{31}}$, в соответствующем случае) так что будет проводиться целочисленное сравнение.
• С учетом двух изменений выше, таблица исходных значений ${\displaystyle x_{i}}$ более не нужна и может быть удалена.
• В случае генерации чисел с плавающей точкой одинарной точности согласно стандарту IEEE 754, где используется 24-битная мантисса (включая неявно заданную 1), наименее значимые биты 32-битного целого случайного числа не используются. Эти биты могут быть использованы при выборе уровня. (тут^[1] подробно описана суть вопроса).

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

Возможно или хранить таблицу с заранее вычисленными ${\displaystyle x_{i}}$ и ${\displaystyle y_{i}}$ целиком, или всего лишь включить значения ${\displaystyle n}$, ${\displaystyle y_{1}}$, ${\displaystyle A}$, и реализацию ${\displaystyle f^{-1}(y)}$ в исходный код, и вычислить оставшиеся значения при инициализации генератора случайных чисел (зависит от того, что нам дороже: вычислительное время или память).

Можно находить ${\displaystyle x_{i}=f^{-1}(y_{i})}$ и ${\displaystyle y_{i+1}=y_{i}+A/x_{i}}$. Повторить ${\displaystyle n-1}$ для всех уровней зиккурата. В конце должно получиться ${\displaystyle y_{n}=f(0)}$.

При итоговом заполнении таблицы нужно положить ${\displaystyle x_{n}=0}$ и ${\displaystyle y_{n}=f(0)}$, приняв небольшие несостыковки (если они и правда вышли небольшими) как ошибки округления.

Поиск ${\displaystyle x_{1}}$ и ${\displaystyle A}$

Если имеется начальное значение ${\displaystyle x_{1}}$ (вычисленное если и не точно, то приближенно), останется лишь вычислить площадь ${\displaystyle t}$ хвостовой части функции, для которой выполнено ${\displaystyle x>x_{1}}$. Вычислить можно методами численного интегрирования.

Далее, из ${\displaystyle x_{1}}$ можно найти ${\displaystyle y_{1}=f(x_{1})}$, из площади ${\displaystyle t}$ хвостовой части найдется площадь базового уровня: ${\displaystyle A=x_{1}y_{1}+t}$.

Затем вычисляется серия ${\displaystyle y_{i}}$ и ${\displaystyle x_{i}}$ как показано выше. Если ${\displaystyle y_{i}>f(0)}$ для любых ${\displaystyle i<n}$, тогда начальное значение ${\displaystyle x_{1}}$ было слишком малым, что привело к большой площади ${\displaystyle A}$. Если ${\displaystyle y_{n}<f(0)}$, тогда начальное значение ${\displaystyle x_{1}}$ было слишком большим.

Учитывая сказанное, можно использовать численное решение уравнений (например, метод бисекции) для нахождения значения ${\displaystyle x_{1}}$_, при котором значение ${\displaystyle y_{n-1}}$ так близко к ${\displaystyle f(0)}$ как только возможно. В качестве альтернативы можно рассматривать и находить значения площади верхнего уровня, ${\displaystyle x_{n-1}(f(0)-y_{n-1})}$, настолько близкие к нужному значению ${\displaystyle A}$ как только возможно.

Примечания

Литература

• 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: сайт

Ссылки

• Реализация алгоритма на языке C для нормальной и показательной функций плотности, по существу, копия кода из статьи.
• Реализация на языке C# и обзор самого алгоритма.
• The Ziggurat Random Normal Generator Blogs of MathWorks, posted by Cleve Moler, May 18, 2015.

Эта страница в последний раз была отредактирована 28 мая 2022 в 17:20.