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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Оптимизации

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

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

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

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

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

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

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

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

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

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

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