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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Оптимизации

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

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

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

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

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

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

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

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

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

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

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

Примечания

Jurgen A. Doornik. (англ.) // 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: