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

Алгоритм зиккурата — это алгоритм выборки псевдослучайных чисел . Принадлежащий к классу алгоритмов выборки отклонения , он опирается на базовый источник равномерно распределенных случайных чисел, обычно получаемых из генератора псевдослучайных чисел , а также на предварительно вычисленные таблицы. Алгоритм используется для генерации значений из монотонно убывающего распределения вероятностей . Его также можно применить к симметричным унимодальным распределениям , таким как нормальное распределение , выбирая значение из одной половины распределения, а затем случайным образом выбирая, какая половина значения считается полученной. Он был разработан Джорджем Марсальей и другими в 1960-х годах.

Типичное значение, полученное с помощью алгоритма, требует создания только одного случайного значения с плавающей запятой и одного случайного индекса таблицы, за которым следуют один поиск в таблице, одна операция умножения и одно сравнение. Иногда (2,5% случаев, в случае нормального или экспоненциального распределения при использовании типичных размеров таблиц) требуются дополнительные вычисления. Тем не менее, алгоритм вычислительно намного быстрее. чем два наиболее часто используемых метода генерации нормально распределенных случайных чисел, полярный метод Марсальи и преобразование Бокса-Мюллера , которые требуют по крайней мере одного логарифма и одного вычисления квадратного корня для каждой пары сгенерированных значений. Однако, поскольку алгоритм зиккурата более сложен в реализации, его лучше всего использовать, когда требуется большое количество случайных чисел.

Термин «алгоритм зиккурата» восходит к статье Марсальи и Вай Ван Цанга в 2000 году; он назван так потому, что концептуально основан на покрытии распределения вероятностей прямоугольными сегментами, сложенными в порядке убывания размера, в результате чего получается фигура, напоминающая зиккурат .

Алгоритм Зиккурата, используемый для генерации выборочных значений с нормальным распределением . (Для простоты показаны только положительные значения.) Розовые точки представляют собой изначально равномерно распределенные случайные числа. Искомая функция распределения сначала сегментируется на равные площади «А». Один слой i выбирается случайным образом однородным источником слева. Затем случайное значение из верхнего источника умножается на ширину выбранного слоя, и результат проверяется по x, чтобы увидеть, в какую область слоя оно попадает, с тремя возможными результатами: 1) (слева, сплошная черная область) образец явно находится под кривой и может быть выведен немедленно. 2) (справа, область с вертикальной полосой) значение выборки может лежать под кривой и должно быть проверено дальше. случайное значение y В этом случае генерируется внутри выбранного слоя и сравнивается с f(x) . Если меньше, точка находится под кривой и значение x выводится . В противном случае (третий случай) выбранная точка x отклоняется и алгоритм перезапускается с начала.

Теория работы

Алгоритм зиккурата представляет собой алгоритм отбраковки выборки; он случайным образом генерирует точку в распределении, немного превышающем желаемое распределение, а затем проверяет, находится ли сгенерированная точка внутри желаемого распределения. Если нет, он пытается еще раз. Учитывая случайную точку под кривой плотности вероятности, ее координата x представляет собой случайное число с желаемым распределением.

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

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

Этот слой (назовем его слоем 0) имеет A. площадь Поверх этого добавьте прямоугольный слой шириной x ₁ и высотой A / x ₁ также имел площадь A. , чтобы он Вершина этого слоя находится на высоте y ₂ = y ₁ + A / x ₁ и пересекает функцию плотности в точке ( x ₂ , y ₂ ), где y ₂ = f ( x ₂ ). Этот слой включает в себя каждую точку функции плотности между y ₁ и y ₂ , но (в отличие от базового слоя) также включает такие точки, как ( x ₁ , y ₂ ), которые не входят в желаемое распределение.

Затем сверху укладываются следующие слои. Чтобы использовать предварительно вычисленную таблицу размера n ( обычно n = 256), выбирают x ₁ так, что x _n = 0, что означает, что верхний блок, слой n - 1, достигает пика распределения в (0, f (0) ) точно.

Слой i простирается вертикально от y _i до y _{i +1} и может быть разделен на две области по горизонтали: (обычно большую) часть от 0 до xi _{, +1} которая полностью содержится в желаемом распределении, и (маленькую) часть. от x _{i +1} до x _i , который содержится лишь частично.

Игнорируя на мгновение проблему слоя 0 и учитывая однородные случайные величины U ₀ и U ₁ ∈ [0,1), алгоритм зиккурата можно описать как:

1. Выберите случайный слой 0 ≤ i < n .
2. Пусть Икс знак равно U ₀ Икс _я .
3. Если x < x _{i +1} , верните x .
4. Пусть y = y _i + U ₁ ( y _{i +1} - y _i ).
5. Вычислите f ( x ). Если y < f ( x ), верните x .
6. В противном случае выберите новые случайные числа и вернитесь к шагу 1.

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

При близко расположенных слоях алгоритм завершается на шаге 3 в очень большой части времени. для верхнего слоя n Однако − 1 этот тест всегда не пройден, поскольку x _n = 0.

Слой 0 также можно разделить на центральную область и край, но край представляет собой бесконечный хвост. Чтобы использовать тот же алгоритм для проверки того, находится ли точка в центральной области, сгенерируйте фиктивный x ₀ = A / y ₁ . Это будет генерировать точки с x < x ₁ с правильной частотой, и в том редком случае, когда выбран слой 0 и x ≥ x ₁ , используется специальный резервный алгоритм для случайного выбора точки из хвоста. Поскольку резервный алгоритм используется реже, чем один раз из тысячи, скорость не имеет существенного значения.

Таким образом, полный алгоритм зиккурата для односторонних распределений выглядит так:

1. Выберите случайный слой 0 ≤ i < n .
2. Пусть x = U ₀ x _i
3. Если x < x _{i +1} , верните x .
4. Если i = 0, сгенерируйте точку из хвоста, используя резервный алгоритм.
5. Пусть y = y _i + U ₁ ( y _{i +1} - y _i ).
6. Вычислите f ( x ). Если y < f ( x ), верните x .
7. В противном случае выберите новые случайные числа и вернитесь к шагу 1.

Для двустороннего распределения результат должен быть отрицательным в 50% случаев. Часто это можно сделать удобно, выбрав U ₀ ∈ (−1,1) и на шаге 3 проверив, | х | < Икс _{я +1} .

Алгоритмы отката для хвоста

Поскольку алгоритм зиккурата генерирует большую часть выходных данных очень быстро и требует резервного алгоритма всякий раз, когда x > x ₁ , он всегда более сложен, чем более прямая реализация. Конкретный алгоритм возврата зависит от распределения.

Для экспоненциального распределения хвост выглядит так же, как тело распределения. Один из способов — вернуться к самому элементарному алгоритму E = −ln( U ₁ ) и положить x = x ₁ − ln( U ₁ ). вызвать алгоритм зиккурата Другой способ — рекурсивно и добавить x ₁ к результату .

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

1. Пусть x = −ln( U ₁ )/ x ₁ .
2. Пусть y −ln( U2 ₌ ).
3. Если 2 у > х ² , верните х + х ₁ .
4. В противном случае вернитесь к шагу 1.

Поскольку x ₁ ≈ 3,5 для типичных размеров таблиц, тест на шаге 3 почти всегда оказывается успешным. Поскольку −ln( U ₁ ) является экспоненциально распределенной переменной, можно использовать реализацию экспоненциального распределения.

Оптимизации

Алгоритм может быть эффективно реализован с предварительно вычисленными таблицами x _i и y _i = f ( x _i ), но есть некоторые модификации, которые сделают его еще быстрее:

• Ничто в алгоритме зиккурата не зависит от нормируемой функции распределения вероятностей (интеграл под кривой, равный 1), удаление нормализующих констант может ускорить вычисление f ( x ).
• Большинство однородных генераторов случайных чисел основаны на генераторах целочисленных случайных чисел, которые возвращают целое число в диапазоне [0, 2 ³² − 1]. Стол из 2 человек ⁻³² x _i позволяет использовать такие числа непосредственно для U ₀ .
• При вычислении двусторонних распределений с использованием двустороннего U _0, как описано ранее, случайное целое число можно интерпретировать как число со знаком в диапазоне [−2 ³¹ , 2 ³¹ − 1] и масштабный коэффициент 2 ⁻³¹ можно использовать.
• Вместо того, чтобы сравнивать U ₀ x _i с x _{i +1} на этапе 3, можно предварительно вычислить x _{i +1} / x _i и сравнить U ₀ с этим напрямую. Если U ₀ — генератор целочисленных случайных чисел, эти пределы можно предварительно умножить на 2. ³² (или 2 ³¹ , при необходимости), поэтому можно использовать целочисленное сравнение.
• После двух вышеуказанных изменений таблица неизмененных xi _{значений} больше не нужна и может быть удалена.
• При генерации значений с плавающей запятой одинарной точности IEEE 754 , которые имеют только 24-битную мантиссу (включая неявную ведущую 1), младшие биты 32-битного целого случайного числа не используются. Эти биты могут использоваться для выбора номера слоя. (Подробное обсуждение этого вопроса см. в ссылках ниже.)
• Первые три шага можно поместить во встроенную функцию , которая может вызывать внешнюю реализацию менее часто необходимых шагов.

Создание таблиц

Можно сохранить всю предварительно вычисленную таблицу или просто включить значения n , y1 _. , A и реализацию f ⁻¹ ( y ) в исходном коде и вычислите оставшиеся значения при инициализации генератора случайных чисел.

Как описано ранее, вы можете найти x _i = f ⁻¹ ( y _я ) и y _{я +1} знак равно y _я + А / Икс _я . Повторите n − 1 раз для слоев зиккурата. В конце у вас должно получиться y _n = f (0). Будет некоторая ошибка округления , но это полезная проверка работоспособности, позволяющая убедиться, что она приемлемо мала.

При фактическом заполнении значений таблицы просто предположите, что x _n = 0 и y _n = f (0), и примите небольшую разницу в площади слоя n - 1 как ошибку округления.

Нахождение x ₁ и A

Учитывая начальное (угадайте) x ₁ , вам нужен способ вычислить площадь t хвоста, для которого x > x ₁ . Для экспоненциального распределения это просто e ^{− х ₁} , а для нормального распределения, предполагая, что вы используете ненормализованное f ( x ) = e ^{− х 2 /2} , это √ π /2   erfc ( x / √ 2 ). Для более неудобных распределений численное интегрирование может потребоваться .

Имея это в виду, из x ₁ вы можете найти y ₁ = f ( x ₁ ), площадь t в хвосте и площадь базового слоя A = x ₁ y ₁ + t .

Затем вычислите ряды y _i и xi _, как указано выше. Если y _i > f (0) для любого i < n оценка x ₁ была слишком низкой, что приводило к слишком большой площади A. , то первоначальная Если y _n < f (0), то первоначальная оценка x ₁ была слишком высокой.

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

Вариант Макфарланда

Кристофер Д. Макфарланд предложил еще более оптимизированную версию. При этом применяются три алгоритмических изменения за счет немного большего размера таблиц.

Во-первых, в общем случае рассматриваются только прямоугольные части от (0, y _{i −1} ) до ( xi _, . y _i ). Области нечетной формы справа от них (в основном почти треугольные, плюс хвост) обрабатываются отдельно . алгоритма Это упрощает и ускоряет работу .

Во-вторых, используется точная площадь областей нестандартной формы; они не округляются в большую сторону, чтобы включить весь прямоугольник до ( x _{i −1} , y _i ). Это увеличивает вероятность того, что будет использован быстрый путь.

Одним из основных последствий этого является то, что количество слоев немного меньше n . Несмотря на то, что площадь частей нестандартной формы взята точно, общая сумма составляет более одного слоя. Площадь каждого слоя регулируется таким образом, чтобы количество прямоугольных слоев было целым числом. Если исходное 0 ≤ i < n превышает количество прямоугольных слоев, продолжается этап 2.

Если искомое значение лежит в какой-либо из областей нечетной формы, метод псевдонима используется для выбора одной из них на основе ее истинной площади. Это небольшой объем дополнительной работы, требующий дополнительных таблиц псевдонимов, но при этом выбирается одна из правых сторон слоев.

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

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

Если функция выпуклая (поскольку экспоненциальное распределение есть везде, а нормальное распределение — для | x | > 1), то функция строго содержится внутри нижнего треугольника. два единичных однородных отклонения U ₁ и U ₂ Выбираются , и перед их масштабированием до прямоугольника, охватывающего область нечетной формы, проверяется их сумма. Если U ₁ + U ₂ > 1, точка находится в верхнем треугольнике и может быть отражена в (1 − U ₁ , 1 − U ₂ ). Тогда, если U ₁ + U ₂ < 1− ε , для некоторого подходящего допуска ε точка определенно находится ниже кривой и может быть немедленно принята. Только для точек, очень близких к диагонали, необходимо вычислить функцию распределения f ( x ), чтобы выполнить точный тест отбраковки. (Теоретически допуск ε должен зависеть от слоя, но одно максимальное значение можно использовать для всех слоев с небольшими потерями.)

Если функция вогнутая (как нормальное распределение для | x | < 1), она включает в себя небольшую часть верхнего треугольника, поэтому отражение невозможно, но точки, нормализованные координаты которых удовлетворяют U ₁ + U ₂ ≤ 1, могут быть немедленно приняты. , а точки, для которых U ₁ + U ₂ > 1+ ε, можно сразу отбросить.

В одном слое, охватывающем | х | = 1, нормальное распределение имеет точку перегиба, и критерий точного отклонения необходимо применять, если 1− ε < U ₁ + U ₂ < 1+ ε .

Хвост обрабатывается так же, как в оригинальном алгоритме Зиккурата, и его можно рассматривать как четвертый случай формы области нечетной формы справа.

Ссылки

• Марсалья, Джордж ; Цанг, Вай Ван (2 октября 2000 г.). «Метод Зиккурата для генерации случайных величин» . Журнал статистического программного обеспечения . 5 (8) . Проверено 20 июня 2007 г. В этой статье слои нумеруются с 1, начиная сверху, а слой 0 внизу рассматривается как особый случай, тогда как в приведенном выше объяснении слои нумеруются с 0 внизу.
• C реализация метода зиккурата для нормальной функции плотности и функции экспоненциальной плотности , которая по сути является копией кода из статьи. (Потенциальные пользователи должны знать, что этот код C предполагает 32-битные целые числа.)
• AC# реализация алгоритма зиккурата и обзор метода.
• Юрген А. Дорник (2005). «Улучшенный метод зиккурата для создания нормальных случайных выборок» (PDF) . Наффилд-колледж, Оксфорд . Проверено 20 июня 2007 г. Описывает опасности использования младших битов генератора целочисленных случайных чисел для выбора номера слоя.
• Нормальное поведение Клива Молера, MathWorks, описывающего алгоритм зиккурата, представленный в MATLAB версии 5, 2001 г.
• Блоги о случайном нормальном генераторе Зиккурата на сайте MathWorks, опубликованные Кливом Молером, 18 мая 2015 г.
• Дэвид Б. Томас; Филип Х.В. Леонг; Уэйн Люк; Джон Д. Вилласенор (октябрь 2007 г.). «Гауссовские генераторы случайных чисел» (PDF) . Обзоры вычислительной техники ACM . 39 (4): 11:1–38. дои : 10.1145/1287620.1287622 . ISSN   0360-0300 . S2CID   10948255 . Проверено 27 июля 2009 г. [Когда] поддержание чрезвычайно высокого статистического качества является главным приоритетом, и с учетом этого ограничения также желательна скорость, метод Зиккурата часто будет наиболее подходящим выбором. Сравнение нескольких алгоритмов генерации гауссовских случайных чисел.
• Надлер, Боаз (2006). «Ошибки проектирования в реализации методов Зиккурата и Монти Пайтона (и некоторые замечания по Matlab randn)». arXiv : математика/0603058 . . Иллюстрирует проблемы с базовыми генераторами унифицированных псевдослучайных чисел и то, как эти проблемы влияют на выходные данные алгоритма зиккурата.
• Эдрис, Хасан М.; Чунг, Брайан; Сандора, Маккаллен; Намми, Дэвид; Стефан, Деян (13–16 июля 2009 г.). Аппаратно-оптимизированный алгоритм зиккурата для высокоскоростных гауссовских генераторов случайных чисел (PDF) . 2009 Международная конференция по разработке реконфигурируемых систем и алгоритмов. Лас Вегас.
• Марсалья, Джордж (сентябрь 1963 г.). Генерация переменной из хвоста нормального распределения (технический отчет). Научно-исследовательские лаборатории Боинга. Математическое примечание № 322, номер доступа DTIC AD0423993. Архивировано из оригинала 10 сентября 2014 г. – через Центр технической информации Министерства обороны .