Отбраковочная выборка

В численном анализе и статистике вычислительной отбраковочная выборка является основным методом, используемым для получения наблюдений из распределения . Его также часто называют методом принятия-отклонения или «алгоритмом принятия-отклонения», и он представляет собой разновидность метода точного моделирования. Метод работает для любого дистрибутива в $\mathbb {R} ^{m}$ с плотностью .

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

Описание

Чтобы визуализировать мотивацию отбраковки выборки, представьте себе, что вы строите график функции плотности вероятности (PDF) случайной величины на большой прямоугольной доске и бросаете в нее дротики. Предположим, что дротики равномерно распределены по доске. Теперь удалите все вытачки, находящиеся за пределами области под кривой. Остальные дротики будут равномерно распределены внутри области под кривой, а $x$‑положения этих дротиков будут распределены в соответствии с плотностью случайной величины. Это связано с тем, что дротикам больше всего места для приземления там, где кривая самая высокая, и, следовательно, плотность вероятности наибольшая.

Только что описанная визуализация эквивалентна особой форме выборки отклонений, при которой «распределение предложений» является равномерным. Следовательно, его график представляет собой прямоугольник. Общая форма выборки отклонений предполагает, что доска не обязательно имеет прямоугольную форму, а имеет форму, соответствующую плотности некоторого распределения предложений (не обязательно нормализованного к $1$), из которого мы знаем, как производить выборку (например, с помощью инверсионной выборки ). Его форма должна быть в каждой точке по крайней мере такой же высокой, как распределение, из которого мы хотим выполнить выборку, чтобы первое полностью охватывало второе. В противном случае были бы части изогнутой области, из которых мы хотим получить образец, и которые никогда не были бы доступны.

Отбраковочная выборка работает следующим образом:

1. Образец точки на $x$— ось распределения предложений.
2. Нарисуйте здесь вертикальную линию $x$-положение, вплоть до максимального значения y функции плотности вероятности распределения предложений.
3. Выполняйте выборку равномерно вдоль этой линии от 0 до максимума функции плотности вероятности. Если выбранное значение больше значения желаемого распределения на этой вертикальной линии, отклоните $x$—оценить и вернуться к шагу 1; еще $x$‑value — это выборка из желаемого распределения.

Этот алгоритм можно использовать для выборки из области под любой кривой, независимо от того, интегрируется ли функция до 1. Фактически, масштабирование функции по константе не влияет на выборку. $x$-позиции . Таким образом, алгоритм можно использовать для выборки из распределения, нормализующая константа которого неизвестна, что часто встречается в вычислительной статистике .

Теория

Метод бракованной выборки генерирует значения выборки из целевого распределения. $X$ с произвольной функцией плотности вероятности $f(x)$ с помощью распределения предложений $Y$ с плотностью вероятности $g(x)$. Идея состоит в том, что можно сгенерировать выборочное значение из $X$ вместо этого осуществляя выборку из $Y$ и принятие образца от $Y$ с вероятностью $f(x)/(Mg(x))$, повторяя розыгрыши из $Y$ пока значение не будет принято. $M$ вот постоянная, конечная граница отношения правдоподобия $f(x)/g(x)$, удовлетворяя $M<\infty$ из- поддержки за $X$; другими словами, M должно удовлетворять $f(x)\leq Mg(x)$ для всех значений $x$. Обратите внимание, что для этого требуется поддержка $Y$ должна включать поддержку $X$-другими словами, $g(x)>0$ в любое время $f(x)>0$.

Валидацией этого метода является принцип конверта: при моделировании пары ${\textstyle (x,v=u\cdot Mg(x))}$, производится однородное моделирование по подграфу ${\textstyle Mg(x)}$. Принимаем только такие пары, что ${\textstyle u<f(x)/(Mg(x))}$ затем создает пары $(x,v)$ равномерно распределены по подграфу $f(x)$ и таким образом, в некоторой степени, симуляция из $f(x).$

Это означает, что при достаточном количестве повторов алгоритм генерирует выборку из желаемого распределения. $f(x)$. Существует ряд расширений этого алгоритма, например алгоритм Метрополиса .

Этот метод относится к общей области методов Монте-Карло , включая алгоритмы Монте-Карло с цепями Маркова , которые также используют прокси-распределение для достижения моделирования на основе целевого распределения. $f(x)$. Он формирует основу для таких алгоритмов, как алгоритм Метрополиса .

Вероятность безусловного принятия — это доля принятых предложенных образцов, которая равна $${\begin{aligned}\mathbb {P} \left(U\leq {\frac {f(Y)}{Mg(Y)}}\right)&=\operatorname {E} \mathbf {1} _{\left[U\leq {\frac {f(Y)}{Mg(Y)}}\right]}\\[6pt]&=E\left[\operatorname {E} [\mathbf {1} _{\left[U\leq {\frac {f(Y)}{Mg(Y)}}\right]}|Y]\right]&({\text{by tower property }})\\[6pt]&=\operatorname {E} \left[\mathbb {P} \left(U\leq {\frac {f(Y)}{Mg(Y)}}{\biggr |}Y\right)\right]\\[6pt]&=E\left[{\frac {f(Y)}{Mg(Y)}}\right]&({\text{because }}\Pr(U\leq u)=u,{\text{when }}U{\text{ is uniform on }}(0,1))\\[6pt]&=\int \limits _{y:g(y)>0}{\frac {f(y)}{Mg(y)}}g(y)\,dy\\[6pt]&={\frac {1}{M}}\int \limits _{y:g(y)>0}f(y)\,dy\\[6pt]&={\frac {1}{M}}&({\text{since support of }}Y{\text{ includes support of }}X)\end{aligned}}$$где $U\sim \mathrm {Unif} (0,1)$, и значение $y$ каждый раз генерируется под функцией плотности $g(\cdot )$ распределения предложения $Y$.

Количество образцов, требуемое от $Y$ таким образом, чтобы получить принятое значение, следует геометрическое распределение с вероятностью $1/M$, что имеет в виду $M$. Интуитивно, $M$ — ожидаемое количество необходимых итераций как мера вычислительной сложности алгоритма.

Перепишите приведенное выше уравнение, $$M={\frac {1}{\mathbb {P} \left(U\leq {\frac {f(Y)}{Mg(Y)}}\right)}}$$ Обратите внимание, что ${\textstyle 1\leq M<\infty }$, по приведенной выше формуле, где ${\textstyle \mathbb {P} \left(U\leq {\frac {f(Y)}{Mg(Y)}}\right)}$ вероятность, которая может принимать значения только в интервале $[0,1]$. Когда $M$ выбирается ближе к единице, то вероятность безусловного принятия тем выше, чем меньше изменяется это отношение, поскольку $M$ это верхняя граница отношения правдоподобия ${\textstyle f(x)/g(x)}$. На практике значение $M$ значение ближе к 1 является предпочтительным, поскольку оно предполагает в среднем меньшее количество отклоненных выборок и, следовательно, меньшее количество итераций алгоритма. В этом смысле предпочитают иметь $M$ как можно меньше (но при этом удовлетворяя $f(x)\leq Mg(x)$, что говорит о том, что $g(x)$ в целом должен напоминать $f(x)$ каким-то образом. Однако обратите внимание, что $M$ не может быть равно 1: это означало бы, что $f(x)=g(x)$, т. е. целевое и предлагаемое распределения на самом деле являются одним и тем же распределением.

Отбраковочная выборка чаще всего применяется в тех случаях, когда форма $f(x)$ затрудняет выборку. Одна итерация алгоритма отклонения требует выборки из распределения предложений, извлечения из равномерного распределения и оценки $f(x)/(Mg(x))$ выражение. Таким образом, отбраковочная выборка более эффективна, чем какой-либо другой метод, если стоимость этих операций, умноженная на M (которая представляет собой ожидаемую стоимость получения выборки с помощью браковочной выборки), ниже, чем стоимость получения выборки с использованием другого метода.

Алгоритм

Алгоритм, который использовал Джон фон Нейман и восходит к Бюффону и его игле , получает образец из распределения $X$ с плотностью $f$ используя образцы из раздачи $Y$ с плотностью $g$ следующее:

• Получить образец $y$ от распространения $Y$ и образец $u$ от $\mathrm {Unif} (0,1)$ (равномерное распределение на единичном интервале).
• Проверьте, если ${\textstyle u<f(y)/Mg(y)}$.
  • Если это так, примите $y$ как образец, взятый из $f$;
  • если нет, отклоните значение $y$ и вернитесь к этапу выборки.

Алгоритм будет принимать в среднем $M$ итерации для получения образца.

Преимущества перед выборкой с использованием наивных методов

В некоторых ситуациях отбраковочная выборка может быть гораздо более эффективной по сравнению с наивными методами. Например, при такой задаче, как выборка ${\textstyle X\sim F(\cdot )}$ условно на $X$ учитывая набор $A$, то есть, ${\textstyle X|X\in A}$, иногда ${\textstyle X}$ можно легко смоделировать, используя наивные методы (например, путем выборки обратного преобразования ):

• Образец ${\textstyle X\sim F(\cdot )}$ независимо и принимать те, которые удовлетворяют $\{n\geq 1:X_{n}\in A\}$
• Выход: $\{X_{1},X_{2},...,X_{N}:X_{i}\in A,i=1,...,N\}$ (см. также усечение (статистика) )

Проблема в том, что выборка может быть сложной и неэффективной, если ${\textstyle \mathbb {P} (X\in A)\approx 0}$. Ожидаемое количество итераций будет ${\frac {1}{\mathbb {P} (X\in A)}}$, которая может быть близка к бесконечности. Более того, даже если вы применяете метод выборки «Отбраковка», всегда сложно оптимизировать границу. $M$ для отношения правдоподобия. Чаще всего, $M$ велик, а процент отказов высок, алгоритм может быть очень неэффективным. Семейство естественных экспонент (если оно существует), также известное как экспоненциальный наклон, представляет собой класс распределений предложений, которые могут снизить сложность вычислений, ценность $M$ и ускорить вычисления (см. примеры: работа с натуральными экспоненциальными семействами).

Отбраковочная выборка с использованием экспоненциального наклона

Учитывая случайную величину $X\sim F(\cdot )$, $F(x)=\mathbb {P} (X\leq x)$ является целевым распределением. Предположим для простоты, что функцию плотности можно явно записать как $f(x)$. Выберите предложение как

$${\begin{aligned}F_{\theta }(x)&=\mathbb {E} \left[\exp(\theta X-\psi (\theta ))\mathbb {I} (X\leq x)\right]\\&=\int _{-\infty }^{x}e^{\theta y-\psi (\theta )}f(y)dy\\g_{\theta }(x)&=F'_{\theta }(x)=e^{\theta x-\psi (\theta )}f(x)\end{aligned}}$$

где $\psi (\theta )=\log \left(\mathbb {E} \exp(\theta X)\right)$ и $\Theta =\{\theta :\psi (\theta )<\infty \}$. Четко, $\{F_{\theta }(\cdot )\}_{\theta \in \Theta }$, принадлежит к естественному экспоненциальному семейству . Более того, отношение правдоподобия

$$Z(x)={\frac {f(x)}{g_{\theta }(x)}}={\frac {f(x)}{e^{\theta x-\psi (\theta )}f(x)}}=e^{-\theta x+\psi (\theta )}$$

Обратите внимание, что $\psi (\theta )<\infty$ подразумевает, что это действительно функция генерации кумулянта , то есть

$\psi (\theta )=\log \mathbb {E} {\exp(tX)}|_{t=\theta }=\log M_{X}(t)|_{t=\theta }$.

Легко вывести функцию генерации кумулянта предложения и, следовательно, кумулянты предложения.

${\begin{aligned}\psi _{\theta }(\eta )&=\log \left(\mathbb {E} _{\theta }\exp(\eta X)\right)=\psi (\theta +\eta )-\psi (\theta )<\infty \\\mathbb {E} _{\theta }(X)&=\left.{\frac {\partial \psi _{\theta }(\eta )}{\partial \eta }}\right|_{\eta =0}\\\mathrm {Var} _{\theta }(X)&=\left.{\frac {\partial ^{2}\psi _{\theta }(\eta )}{\partial ^{2}\eta }}\right|_{\eta =0}\end{aligned}}$

В качестве простого примера предположим, что под $F(\cdot )$, $X\sim \mathrm {N} (\mu ,\sigma ^{2})$, с ${\textstyle \psi (\theta )=\mu \theta +{\frac {\sigma ^{2}\theta ^{2}}{2}}}$. Цель – взять образец $X|X\in \left[b,\infty \right]$, где $b>\mu$. Анализ происходит следующим образом:

• Выберите форму распространения предложения $F_{\theta }(\cdot )$, с кумулянт-генерирующей функцией как

${\textstyle \psi _{\theta }(\eta )=\psi (\theta +\eta )-\psi (\theta )=(\mu +\theta \sigma ^{2})\eta +{\frac {\sigma ^{2}\eta ^{2}}{2}}}$,

что далее означает, что это нормальное распределение $\mathrm {N} (\mu +\theta \sigma ^{2},\sigma ^{2})$.

• Решите, что хорошо выбрано $\theta ^{*}$ для распространения предложения. В этой настройке интуитивный способ выбора $\theta ^{*}$ это установить

$\mathbb {E} _{\theta }(X)=\mu +\theta \sigma ^{2}=b$,

то есть $\theta ^{*}={\frac {b-\mu }{\sigma ^{2}}}.$ Распределение предложений таким образом $g_{\theta ^{*}}(x)=\mathrm {N} (b,\sigma ^{2})$.

• Четко запишите цель, предложение и отношение правдоподобия.

$${\begin{aligned}f_{X|X\geq b}(x)&={\frac {f(x)\mathbb {I} (x\geq b)}{\mathbb {P} (X\geq b)}}\\g_{\theta ^{*}}(x)&=f(x)\exp(\theta ^{*}x-\psi (\theta ^{*}))\\Z(x)&={\frac {f_{X|X\geq b}(x)}{g_{\theta ^{*}}(x)}}={\frac {\exp(-\theta ^{*}x+\psi (\theta ^{*}))\mathbb {I} (x\geq b)}{\mathbb {P} (X\geq b)}}\end{aligned}}$$

• Выведите границу $M$ для отношения правдоподобия $Z(x)$, что является убывающей функцией для $x\in [b,\infty ]$, поэтому

$$M=Z(b)={\frac {\exp(-\theta ^{*}b+\psi (\theta ^{*}))}{\mathbb {P} (X\geq b)}}={\frac {\exp \left(-{\frac {(b-\mu )^{2}}{2\sigma ^{2}}}\right)}{\mathbb {P} (X\geq b)}}={\frac {\exp \left(-{\frac {(b-\mu )^{2}}{2\sigma ^{2}}}\right)}{\mathbb {P} \left(\mathrm {N} (0,1)\geq {\frac {b-\mu }{\sigma }}\right)}}$$

• Критерий браковки: для $U\sim \mathrm {Unif} (0,1)$, если

$$U\leq {\frac {Z(x)}{M}}=e^{-\theta ^{*}(x-b)}\mathbb {I} (x\geq b)$$

выполняется, примите значение $X$; если нет, продолжайте пробовать новые ${\textstyle X\sim _{i.i.d.}\mathrm {N} (\mu +\theta ^{*}\sigma ^{2},\sigma ^{2})}$ и новый ${\textstyle U\sim \mathrm {Unif} (0,1)}$ до принятия.

Для приведенного выше примера в качестве измерения эффективности, ожидаемого числа итераций, имеет порядок отбраковочная выборка на основе естественного экспоненциального семейства. $b$, то есть $M(b)=O(b)$, а при наивном методе ожидаемое количество итераций равно ${\textstyle {\frac {1}{\mathbb {P} (X\geq b)}}=O(b\cdot e^{\frac {(b-\mu )^{2}}{2\sigma ^{2}}})}$, что гораздо более неэффективно.

В общем, экспоненциальный наклон параметрического класса распределения предложений удобно решает задачи оптимизации благодаря своим полезным свойствам, которые непосредственно характеризуют распределение предложений. Для такого типа задач, чтобы смоделировать $X$ условно на $X\in A$В классе простых распределений хитрость заключается в использовании натурального экспоненциального семейства, которое помогает получить некоторый контроль над сложностью и значительно ускорить вычисления. Действительно, существуют глубокие математические причины для использования натурального экспоненциального семейства.

Недостатки

Отбраковочная выборка требует знания целевого распределения (в частности, способности оценить целевой PDF-файл в любой момент).

Отбраковочная выборка может привести к взятию большого количества нежелательных выборок, если выборка функции сильно сконцентрирована в определенной области, например функция, имеющая пик в каком-то месте. Для многих распределений эту проблему можно решить с помощью адаптивного расширения (см. адаптивную браковочную выборку ), либо соответствующей заменой переменных методом отношения униформ . Кроме того, по мере увеличения размеров задачи отношение внедренного объема к «углам» внедренного объема стремится к нулю, поэтому до того, как будет сгенерирован полезный образец, может произойти множество отказов, что делает алгоритм неэффективно и непрактично. См. проклятие размерности . В больших размерностях необходимо использовать другой подход, обычно метод Монте-Карло с цепью Маркова, такой как выборка Метрополиса или выборка Гиббса . (Однако выборка Гиббса, которая разбивает задачу многомерной выборки на серию выборок низкой размерности, может использовать браковочную выборку в качестве одного из этапов.)

Адаптивная отбраковочная выборка

Для многих дистрибутивов сложно найти дистрибутив предложения, включающий данный дистрибутив и не занимающий много места. Расширение браковочной выборки, которое можно использовать для преодоления этой трудности и эффективной выборки из широкого спектра распределений (при условии, что они имеют логарифмически вогнутые функции плотности, что фактически справедливо для большинства распространенных распределений, даже тех, плотность которых функции сами по себе не вогнуты) известен как адаптивная режекционная выборка (ARS) .

В этой технике есть три основные идеи, которые в конечном итоге предложил Гилкс в 1992 году:

1. Если это поможет, вместо этого определите распределение конвертов в пространстве журнала (например, логарифмическую вероятность или логарифмическую плотность). То есть работать с $h\left(x\right)=\log g\left(x\right)$ вместо $g\left(x\right)$ напрямую.
   • Часто распределения, которые имеют алгебраически беспорядочные функции плотности, имеют достаточно простые функции логарифмической плотности (т.е. когда $f\left(x\right)$ это грязно, $\log f\left(x\right)$ может быть проще работать или, по крайней мере, ближе к кусочно-линейному).
2. Вместо единственной функции однородной плотности оболочки используйте в качестве оболочки кусочно-линейную функцию плотности.
   • Каждый раз, когда вам приходится отклонять образец, вы можете использовать значение $f\left(x\right)$ что вы оценили, чтобы улучшить кусочное приближение $h\left(x\right)$. Таким образом, это снижает вероятность того, что ваша следующая попытка будет отклонена. Асимптотически вероятность того, что вам придется отклонить вашу выборку, должна стремиться к нулю, и на практике часто очень быстро.
   • Как предложено, каждый раз, когда мы выбираем отвергнутую точку, мы сжимаем огибающую другим отрезком линии, касательным к кривой в точке с той же координатой x, что и выбранная точка.
   • Кусочно-линейная модель распределения журнала предложений приводит к набору кусочно -экспоненциальных распределений (т.е. сегментов одного или нескольких экспоненциальных распределений, соединенных встык). Экспоненциальные распределения хорошо себя ведут и хорошо понятны. Логарифм экспоненциального распределения представляет собой прямую линию, и, следовательно, этот метод по существу включает в себя заключение логарифма плотности в ряд отрезков линии. Это источник логарифмически вогнутого ограничения: если распределение логарифмически вогнутое, то его логарифм вогнутый (имеет форму перевернутой буквы U), а это означает, что отрезок прямой, касательный к кривой, всегда будет проходить над кривой.
   • Если вы не работаете в логарифмическом пространстве, кусочно-линейную функцию плотности также можно выбрать с помощью треугольных распределений.
3. Мы можем еще больше воспользоваться требованием (логарифмической) вогнутости, чтобы потенциально избежать затрат на оценку $f\left(x\right)$ Когда ваш образец будет принят.
   • Точно так же, как мы можем построить кусочно-линейную верхнюю границу (функцию «огибающая»), используя значения $h\left(x\right)$ которые нам пришлось оценить в текущей цепочке отклонений, мы также можем построить кусочно-линейную нижнюю границу («сжимающую» функцию), используя также эти значения.
   • Прежде чем оценивать (потенциально дорогое) $f\left(x\right)$ Чтобы узнать, будет ли принят ваш образец, мы, возможно, уже узнаем , будет ли он принят, путем сравнения с (в идеале более дешевым) $g_{l}\left(x\right)$ (или $h_{l}\left(x\right)$ в данном случае) функция сжатия, которая есть в наличии.
   • Этот этап сжатия не является обязательным, даже если он предложен Гилксом. В лучшем случае это избавит вас от одной дополнительной оценки вашей (беспорядочной и/или дорогостоящей) целевой плотности. Однако, по-видимому, для особенно дорогих функций плотности (и при условии быстрого стремления процента брака к нулю) это может существенно повлиять на конечное время выполнения.

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

К сожалению, ARS можно применять только для отбора проб из логарифмически вогнутых целевых плотностей. По этой причине в литературе было предложено несколько расширений ARS для решения задач нелогарифмически вогнутых целевых распределений. Кроме того, были разработаны различные комбинации ARS и метода Метрополиса-Гастингса, чтобы получить универсальный пробоотборник, который создает самонастраивающуюся плотность предложений (т. е. предложение автоматически создается и адаптируется к цели). Этот класс методов часто называют алгоритмами Adaptive Rejection Metropolis Sampling (ARMS) . Полученные адаптивные методы можно применять всегда, но в этом случае сгенерированные выборки коррелируют (хотя корреляция быстро исчезает до нуля по мере роста числа итераций).

См. также

• Выборка обратного преобразования
• Соотношение униформы
• Выборка псевдослучайных чисел
• Алгоритм зиккурата

Дальнейшее чтение

• Роберт, CP; Казелла, Г. (2004). Статистические методы Монте-Карло (второе изд.). Нью-Йорк: Springer-Verlag.