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

Ziggurat algorithm

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

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

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

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

Теория Операции

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

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

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

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

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

Слой я простирается вертикально от уя к уя+1, и может быть разделен на две области по горизонтали: (как правило, большую) часть от 0 до Икся+1 который полностью содержится в желаемом распределении, а (небольшая) часть из Икся+1 к Икся, который содержится лишь частично.

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

  1. Выберите случайный слой 0 ≤ я < п.
  2. Позволять Икс = U0Икся.
  3. Если Икс < Икся+1, возвращаться Икс.
  4. Позволять у = уя + U1(уя+1уя).
  5. Вычислить ж(Икс). Если у < ж(Икс), возвращаться Икс.
  6. В противном случае выберите новые случайные числа и вернитесь к шагу 1.

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

В случае близко расположенных слоев алгоритм завершается на шаге 3 очень большую часть времени. Для верхнего слоя п - 1, однако этот тест всегда терпит неудачу, потому что Иксп = 0.

Слой 0 также можно разделить на центральную область и край, но край - это бесконечный хвост. Чтобы использовать тот же алгоритм для проверки, находится ли точка в центральной области, сгенерируйте фиктивный Икс0 = А/у1. Это будет генерировать точки с Икс < Икс1 с правильной частотой, и в редких случаях выбирается слой 0 и ИксИкс1используйте специальный резервный алгоритм для случайного выбора точки из хвоста. Поскольку резервный алгоритм используется реже одного раза из тысячи, скорость не имеет значения.

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

  1. Выберите случайный слой 0 ≤ я < п.
  2. Позволять Икс = U0Икся
  3. Если Икс < Икся+1, возвращаться Икс.
  4. Если я = 0, сгенерировать точку из хвоста, используя резервный алгоритм.
  5. Позволять у = уя + U1(уя+1уя).
  6. Вычислить ж(Икс). Если у < ж(Икс), возвращаться Икс.
  7. В противном случае выберите новые случайные числа и вернитесь к шагу 1.

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

Запасные алгоритмы для хвоста

Поскольку алгоритм зиккурата генерирует только наиболее выводится очень быстро и требует резервного алгоритма всякий раз, когда Икс > Икс1, это всегда сложнее, чем более прямая реализация. Алгоритм отката, конечно, зависит от распределения.

Для экспоненциального распределения хвост выглядит точно так же, как тело распределения. Один из способов - вернуться к простейшему алгоритму E = −ln (U1) и разреши Икс = Икс1 - ln (U1). Другой - вызвать алгоритм зиккурата рекурсивно и добавить Икс1 к результату.

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

  1. Позволять Икс = −ln (U1)/Икс1.
  2. Позволять у = −ln (U2).
  3. Если 2у > Икс2, возвращаться Икс + Икс1.
  4. В противном случае вернитесь к шагу 1.

С Икс1 ≈ 3,5 для типичных размеров таблиц, тест на шаге 3 почти всегда проходит успешно.

Оптимизация

Алгоритм может быть эффективно выполнен с предварительно вычисленными таблицами Икся и уя = ж(Икся), но есть некоторые модификации, чтобы сделать это еще быстрее:

  • В алгоритме зиккурата ничего не зависит от нормализуемой функции распределения вероятностей (интеграл под кривой равен 1), удаляя нормализующие константы может ускорить вычисление ж(Икс).
  • Большинство однородных генераторов случайных чисел основаны на генераторах целых случайных чисел, которые возвращают целое число в диапазоне [0, 232 - 1]. Стол из 2−32Икся позволяет использовать такие числа напрямую для U0.
  • При вычислении двусторонних распределений с использованием двустороннего U0 как описано ранее, случайное целое число можно интерпретировать как число со знаком в диапазоне [−231, 231 - 1] и масштабный коэффициент 2−31 может быть использован.
  • Вместо того, чтобы сравнивать U0Икся к Икся+1 на шаге 3 можно предварительно вычислить Икся+1/Икся и сравните U0 с этим напрямую. Если U0 является генератором целых случайных чисел, эти пределы можно умножить на 232 (или 231, в зависимости от ситуации), поэтому можно использовать целочисленное сравнение.
  • С двумя вышеупомянутыми изменениями таблица неизмененных Икся values ​​больше не нужны и могут быть удалены.
  • При создании IEEE 754 значения с плавающей запятой одинарной точности, которые имеют только 24-битную мантиссу (включая неявную начальную 1), младшие биты 32-битного целого случайного числа не используются. Эти биты могут использоваться для выбора номера слоя. (См. Ссылки ниже для подробного обсуждения этого.)
  • Первые три шага можно поместить в встроенная функция, который может вызвать автономную реализацию менее часто необходимых шагов.

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

Можно сохранить всю предварительно вычисленную таблицу или просто включить значения п, у1, А, и реализация ж −1(у) в исходном коде и вычислить оставшиеся значения при инициализации генератора случайных чисел.

Как описано ранее, вы можете найти Икся = ж −1(уя) и уя+1уя + А/Икся. Повторение п - 1 раз за слои зиккурата. В конце у вас должно быть уп = ж(0). Конечно, будут ошибка округления, но это полезный тест на вменяемость чтобы убедиться, что он достаточно мал.

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

обнаружение Икс1 и А

Учитывая начальное (предположение) Икс1, вам нужен способ вычисления площади т хвоста, для которого Икс > Икс1. Для экспоненциального распределения это просто еИкс1, а для нормального распределения, предполагая, что вы используете ненормализованное ж(Икс) = еИкс2/2, это π/2erfc(Икс/2). Для более неудобных распределений численное интегрирование может потребоваться.

С этим в руках от Икс1, ты можешь найти у1 = ж(Икс1), площадь т в хвосте, а область базового слоя А = Икс1у1 + т.

Затем вычислите серию уя и Икся как указано выше. Если уя > ж(0) для любого я < п, то начальная оценка Икс1 был слишком низким, что привело к слишком большой площади А. Если уп < ж(0), то начальная оценка Икс1 было слишком высоко.

Учитывая это, используйте алгоритм поиска корней (такой как метод деления пополам), чтобы найти значение Икс1 который производит уп−1 так близко к ж(0) насколько это возможно. В качестве альтернативы найдите значение, которое делает область самого верхнего слоя, Иксп−1(ж(0) − уп−1), как можно ближе к желаемому значению А насколько возможно. Это экономит одну оценку ж −1(Икс) и на самом деле представляет наибольший интерес.

Рекомендации

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