WikiSort.ru - Не сортированное

Алгоритм «Зиккурат» (англ. 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}}$ .

Примечания

Литература

• 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.