Muestreo de rechazo

En análisis numéricos y estadísticas computacionales, muestreo de rechazo es una técnica básica utilizada para generar observaciones de una distribución. También se llama comúnmente método de aceptación-rechacción o " algoritmo de aceptación-reyecto" y es un tipo de método de simulación exacta. El método funciona para cualquier distribución en ${displaystyle mathbb {R} {} {m}$ con densidad.

El muestreo de rechazo se basa en la observación de que para muestrear una variable aleatoria en una dimensión, se puede realizar un muestreo aleatorio uniforme del gráfico cartesiano bidimensional y mantener las muestras en la región bajo el gráfico de su función de densidad. Tenga en cuenta que esta propiedad se puede extender a funciones de dimensión N.

Descripción

Para visualizar la motivación detrás del muestreo de rechazo, imagine graficar la función de densidad de una variable aleatoria en un tablero rectangular grande y lanzarle dardos. Suponga que los dardos están distribuidos uniformemente por el tablero. Ahora retira todos los dardos que están fuera del área debajo de la curva. Los dardos restantes se distribuirán uniformemente dentro del área bajo la curva y las posiciones x de estos dardos se distribuirán según la densidad de la variable aleatoria. Esto se debe a que hay más espacio para que los dardos caigan donde la curva es más alta y, por lo tanto, la densidad de probabilidad es mayor.

La visualización como se acaba de describir es equivalente a una forma particular de muestreo de rechazo donde la "distribución de la propuesta" es uniforme (por lo tanto su gráfica es un rectángulo). La forma general de muestreo de rechazo supone que el tablero no es necesariamente rectangular sino que tiene forma de acuerdo con la densidad de alguna distribución de propuesta de la que sabemos cómo tomar muestras (por ejemplo, usando muestreo inverso), y que es al menos tan alta en cada punto. punto como la distribución de la que queremos tomar la muestra, de modo que el primero encierre completamente al segundo. (De lo contrario, habría partes del área curva de la que queremos tomar muestras que nunca podrían alcanzarse).

El muestreo de rechazo funciona de la siguiente manera:

1. Muestra un punto sobre el eje x de la distribución de la propuesta.
2. Dibujar una línea vertical en esta x-posición, hasta el máximo valor y de la función de densidad de probabilidad de la distribución de la propuesta.
3. Muestra uniformemente a lo largo de esta línea de 0 al máximo de la función de densidad de probabilidad. Si el valor muestrado es mayor que el valor de la distribución deseada en esta línea vertical, rechaza el valor x y vuelve al paso 1; de lo contrario el valor x es una muestra de la distribución deseada.

Este algoritmo se puede utilizar para tomar muestras del área bajo cualquier curva, independientemente de si la función se integra a 1. De hecho, escalar una función mediante una constante no tiene ningún efecto en las posiciones x muestreadas. Por tanto, el algoritmo se puede utilizar para tomar muestras de una distribución cuya constante de normalización se desconoce, lo cual es común en estadística computacional.

Teoría

El método de muestreo de rechazo genera valores de muestreo de una distribución objetivo ${displaystyle X}$ con función de densidad de probabilidad arbitraria ${displaystyle f(x)}$ utilizando una distribución de propuestas ${displaystyle Sí.$ con densidad de probabilidad ${displaystyle g(x)}$. La idea es que se puede generar un valor de muestra ${displaystyle X}$ por en lugar del muestreo ${displaystyle Sí.$ y aceptar la muestra de ${displaystyle Sí.$ con probabilidad ${displaystyle f(x)/(Mg(x)}$, repitiendo los sorteos de ${displaystyle Sí.$ hasta que se acepte un valor. ${displaystyle M}$ aquí es una constante, finita ligada a la relación de probabilidad ${displaystyle f(x)/g(x)}$, satisfacción ${displaystyle 1 = 0 = 0$ sobre el apoyo de ${displaystyle X}$; en otras palabras, M debe satisfacer ${displaystyle f(x)leq Mg(x)}$ para todos los valores ${displaystyle x}$. Tenga en cuenta que esto requiere que el apoyo de ${displaystyle Sí.$ debe incluir el apoyo ${displaystyle X}$—en otras palabras, ${displaystyle g(x)}$ siempre ${displaystyle f(x)}$.

La validación de este método es el principio del sobre: al simular el par ${textstyle (x,v=ucdot Mg(x)}$, uno produce una simulación uniforme sobre el subgrafo de ${textstyle Mg(x)}$. Aceptar sólo pares tales que ${textstyle u efectuadaf(x)/(Mg(x)}$ entonces produce pares ${displaystyle (x,v)}$ distribuido uniformemente sobre el subgrafo ${displaystyle f(x)}$ y por lo tanto, marginalmente, una simulación de ${displaystyle f(x).}$

Esto significa que, con suficientes réplicas, el algoritmo genera una muestra de la distribución deseada ${displaystyle f(x)}$. Hay una serie de extensiones a este algoritmo, como el algoritmo de Metropolis.

Este método se relaciona con el campo general de las técnicas de Monte Carlo, incluyendo algoritmos de la cadena Markov Monte Carlo que también utilizan una distribución proxy para lograr la simulación de la distribución objetivo ${displaystyle f(x)}$. Forma la base para algoritmos como el algoritmo de Metropolis.

La probabilidad de aceptación incondicional es la proporción de muestras propuestas que se aceptan, que es

$${displaystyle {begin{aligned}mathbb {P} left(Uleq {f(Y)}{Mg(Y)}right) Pulse=operatorname {E} mathbf {1} _{left[Uleq {frac {f(Y)}{Mg(Y)}}right]}\[6pt] {E} [mathbf {1} _{left [Uleq {frac {f(Y)}{Mg(Y)}}}right]right]right] {fnMicrosoft} {fnMicrosoft} {fnMicrosoft} {f} {f}} {f}}} {fnMicrosoft}}} {f}} {fnMicrosoft} {fnMicrosoft}}} {f}}} {f}f}f}f}f}f}f}f} {f}f}f}f}f}f}f}f}f}f}f}f}f} {f}f}f}f}f}f}f}f}f}f}f} {f}f}f}f}f}f} {f}f} {f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}$$

${displaystyle Usim mathrm {Unif} (0,1)}$

${displaystyle y}$

${displaystyle g(cdot)}$

${displaystyle Sí.$

Número de muestras requeridas ${displaystyle Sí.$ para obtener un valor aceptado así sigue una distribución geométrica con probabilidad ${displaystyle 1/M}$, que tiene significado ${displaystyle M}$. Intuitivamente, ${displaystyle M}$ es el número esperado de las iteraciones que se necesitan, como medida de la complejidad computacional del algoritmo.

Reescribe la ecuación anterior,

$${displaystyle M={frac {1}{mthbb {}left(Uleq {frac {f(Y)}{Mg(Y)}}}}}}}}}}$$

${textstyle 1leq - No.$

${textstyle mathbb {P} left(Uleq {frac {f(Y)}{Mg(Y)}}right)}}$

${displaystyle [0,1]}$

${displaystyle M}$

${displaystyle M}$

${textstyle f(x)/g(x)}$

${displaystyle M}$

${displaystyle M}$

${displaystyle f(x)leq Mg(x)}$

${displaystyle g(x)}$

${displaystyle f(x)}$

${displaystyle M}$

${displaystyle f(x)=g(x)}$

El muestreo de rechazo se utiliza más a menudo en casos en que la forma de ${displaystyle f(x)}$ hace que el muestreo sea difícil. Una única iteración del algoritmo de rechazo requiere muestreo de la distribución de la propuesta, a partir de una distribución uniforme, y evaluación de la ${displaystyle f(x)/(Mg(x)}$ expresión. Por lo tanto, el muestreo de rechazo es más eficiente que cualquier otro método cada vez que el coste de estas operaciones —que es el costo esperado de obtener una muestra con muestreo de rechazo— es menor que el costo de obtener una muestra usando el otro método.

Algoritmo

El algoritmo, que fue utilizado por John von Neumann y se remonta a Buffon y su aguja, obtiene una muestra de la distribución ${displaystyle X}$ con densidad ${displaystyle f}$ utilizando muestras de distribución ${displaystyle Sí.$ con densidad ${displaystyle g}$ como sigue:

• Obtenga una muestra ${displaystyle y}$ de distribución ${displaystyle Sí.$ y una muestra ${displaystyle u}$ desde ${displaystyle mathrm {Unif} (0,1)}$ (la distribución uniforme en el intervalo de unidad).
• Compruebe si ${textstyle u wonf(y)/Mg(y)}$.
  • Si esto sostiene, acepta ${displaystyle y}$ como muestra extraída ${displaystyle f}$;
  • si no, rechazar el valor ${displaystyle y}$ y volver al escalón.

El algoritmo tomará un promedio de ${displaystyle M}$ iteraciones para obtener una muestra.

Ventajas sobre el muestreo utilizando métodos ingenuos

El muestreo de rechazo puede ser mucho más eficiente en comparación con los métodos ingenuos en algunas situaciones. Por ejemplo, dado un problema como muestreo ${textstyle Xsim F(cdot)}$ condicionalmente ${displaystyle X}$ dado el conjunto ${displaystyle A}$, es decir, ${textstyle X TorturaXin A}$, a veces ${textstyle X}$ se puede simular fácilmente, utilizando los métodos ingenuos (por ejemplo, mediante muestreo de transformación inversa):

• Muestra ${textstyle Xsim F(cdot)}$ independientemente, y aceptar a los que satisfacen ${displaystyle {ngeq 1:X_{n}in A}$
• Producto: ${displaystyle ¿Qué? A,i=1,...$ (véase también truncación (estadística))

El problema es que este muestreo puede ser difícil e ineficiente, si ${textstyle mathbb {P} (Xin A)approx 0}$. El número previsto de iteraciones sería ${fnMicroc} {fnMithbb {P} (Xin A)}}$, que podría estar cerca del infinito. Además, incluso cuando se aplica el método de muestreo de rechazo, siempre es difícil optimizar el límite ${displaystyle M}$ por la relación de probabilidad. Más a menudo que no, ${displaystyle M}$ es grande y la tasa de rechazo es alta, el algoritmo puede ser muy ineficiente. La Familia Exponencial Natural (si existe), también conocida como inclinación exponencial, proporciona una clase de distribuciones de propuestas que pueden reducir la complejidad de la computación, el valor de ${displaystyle M}$ y acelerar las computaciones (ver ejemplos: trabajar con Familias Exponenciales Naturales).

Muestreo de rechazo mediante inclinación exponencial

Dada una variable aleatoria ${displaystyle Xsim F(cdot)}$, ${displaystyle F(x)=mathbb {P} (Xleq x)}$ es la distribución de objetivos. Asumo para la simplicidad, la función de densidad se puede escribir explícitamente como ${displaystyle f(x)}$. Elija la propuesta como

$${fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {f} {fnMicrosoft Sans Serif} {f} {fnMicrosoft Sans Serif} {f} {f} {f}f}fnMicrosoft} {fnMicrosoft}f}f}f} {f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}fnMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMinMin$$

Donde ${displaystyle psi (theta)=log left(mathbb {E} exp(theta X)right)}$ y ${displaystyle Theta ={theta:psi (theta)traducido$. Claramente, ${displaystyle {F_{theta }(cdot)}_{theta in Theta }$, es de una familia exponencial natural. Además, la relación de probabilidad es

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

Note que ${displaystyle psi (theta) se llevó a cabo$ implica que de hecho es una función de generación acumulada, es decir,

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

Es fácil derivar la función de generación de acumuladores de la propuesta y, por lo tanto, los cumulantes de la propuesta.

${displaystyle {begin{aligned}psi _{theta }(eta)}=log left(mathbb {E} _{theta }exp(eta X)right)=psi (theta +eta)-psi (theta)se hizo \infty\\Mathbb {E} _{theta }(X) Alguien=left.{frac {partial psi _{theta }(eta)}{partial eta. =0}\\mathrm [Var] _{theta }(X) Sentido=left.{frac {partial ^{2}psi _{theta } {eta)}{partial ^{2}eta]. {fnMicrosoft Sans}$

Como ejemplo simple, supongamos bajo ${displaystyle F(cdot)}$, ${displaystyle Xsim mathrm {N}$, con ${textstyle psi (theta)=mu theta {fnK}$. El objetivo es probar ${displaystyle X TorturaXin left[b,infty right]}$, donde ${displaystyle b$. El análisis es el siguiente:

• Elija la forma de la distribución de la propuesta ${displaystyle F_{theta}(cdot)}$, con función acumulativa como

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

que implica además que es una distribución normal ${displaystyle mathrm {N} (mu +theta sigma ^{2},sigma ^{2})}$.

• Decidir el pozo elegido ${displaystyle theta ^{*}$ para la distribución de la propuesta. En esta configuración, la forma intuitiva de elegir ${displaystyle theta ^{*}$ es para establecer

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

eso es ${displaystyle theta ^{*}={frac {b-mu - Sí.$ Así pues, la distribución de la propuesta es ${displaystyle g_{theta ^{*}(x)=mathrm {N} (b,sigma ^{2})}$.

• Escribe explícitamente el objetivo, la propuesta y la relación de probabilidad

$${fnMicrosoft Sans Serif}*$$

• Derive el límite ${displaystyle M}$ para la relación de probabilidad ${displaystyle Z(x)}$, que es una función decreciente para ${displaystyle xin [b,infty]$, por lo tanto

$${fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif}$$

• criterio de muestreo de rechazo: para ${displaystyle Usim mathrm {Unif} (0,1)}$, si

$${displaystyle ¿Qué?$$

posee, acepta el valor ${displaystyle X}$; si no, continuar muestreando nuevo ${textstyle Xsim _{i.i.d.}mathrm {N} (mu +theta ^{*}sigma ^{2},sigma ^{2})}$ y nuevos ${textstyle Usim mathrm {Unif} (0,1)}$ hasta la aceptación.

Para el ejemplo anterior, como la medición de la eficiencia, el número esperado de las iteraciones el método de muestreo de rechazo exponencial natural de la familia es de orden ${displaystyle b}$, eso es ${displaystyle M(b)=O(b)}$, mientras que bajo el método ingenuo, el número esperado de las iteraciones es ${textstyle {frac {1}{mthbb {fnMicrosoft Sans Serif}}{2}}}}$, que es mucho más ineficiente.

En general, la inclinación exponencial de una clase paramétrica de distribución de propuestas, resuelve convenientemente los problemas de optimización, con sus propiedades útiles que caracterizan directamente la distribución de la propuesta. Para este tipo de problema, simular ${displaystyle X}$ condicionalmente ${displaystyle Xin A}$, entre la clase de distribuciones simples, el truco es utilizar la familia exponencial natural, que ayuda a ganar cierto control sobre la complejidad y a acelerar considerablemente el cálculo. De hecho, hay razones matemáticas profundas para usar la familia exponencial natural.

Inconvenientes

El muestreo de rechazo requiere conocer la distribución objetivo (específicamente, la capacidad de evaluar el PDF objetivo en cualquier punto).

El muestreo de rechazo puede dar lugar a que se tomen muchas muestras no deseadas si la función que se está muestreando está muy concentrada en una región determinada, por ejemplo, una función que tiene un pico en algún lugar. Para muchas distribuciones, este problema se puede resolver utilizando una extensión adaptativa (ver muestreo de rechazo adaptativo), o con un cambio apropiado de variables con el método de la proporción de uniformes. Además, a medida que las dimensiones del problema aumentan, la relación entre el volumen incrustado y las "esquinas" aumenta. del volumen de inclusión tiende a cero, por lo que pueden producirse muchos rechazos antes de que se genere una muestra útil, lo que hace que el algoritmo sea ineficiente y poco práctico. Ver maldición de la dimensionalidad. En dimensiones altas, es necesario utilizar un enfoque diferente, típicamente un método Monte Carlo de cadena de Markov, como el muestreo de Metropolis o el muestreo de Gibbs. (Sin embargo, el muestreo de Gibbs, que descompone un problema de muestreo multidimensional en una serie de muestras de baja dimensión, puede utilizar el muestreo de rechazo como uno de sus pasos).

Muestreo de rechazo adaptativo

Para muchas distribuciones, es difícil encontrar una distribución propuesta que incluya la distribución dada sin desperdiciar mucho espacio. Una extensión del muestreo de rechazo que se puede utilizar para superar esta dificultad y muestrear eficientemente a partir de una amplia variedad de distribuciones (siempre que tengan funciones de densidad log-cóncavas, que de hecho es el caso de la mayoría de las distribuciones comunes, incluso aquellas cuyo < Las funciones i>densidad no son cóncavas en sí mismas) se conoce como muestreo de rechazo adaptativo (ARS).

Hay tres ideas básicas para esta técnica, tal como la introdujo finalmente Gilks en 1992:

1. Si ayuda, defina la distribución de su sobre en el espacio de registro (por ejemplo, la probabilidad de registro o la densidad de registro). Es decir, trabajar con ${displaystyle hleft(xright)=log gleft(xright)}$ en lugar de ${displaystyle gleft(xright)}$ directamente.
   • A menudo, las distribuciones que tienen funciones de densidad algebraicamente desordenada tienen funciones de densidad de troncos razonablemente más simples (es decir, cuando ${displaystyle fleft(xright)}$ es desordenado, ${displaystyle log fleft(xright)}$ puede ser más fácil trabajar con o, por lo menos, más cerca de linear de la pieza.
2. En lugar de una única función uniforme de densidad de sobre, utilice una función de densidad lineal en sentido estricto como su sobre.
   • Cada vez que usted tiene que rechazar una muestra, puede utilizar el valor de ${displaystyle fleft(xright)}$ que usted evaluó, para mejorar la aproximación de la pieza ${displaystyle hleft(xright)}$. Esto reduce la posibilidad de que su próximo intento sea rechazado. Asintotically, the probability of needing to reject your sample should converge to cero, and in practice, often very quickly.
   • Como se propone, cada vez que elegimos un punto que es rechazado, apretamos el sobre con otro segmento de línea que es tangente a la curva en el punto con el mismo x-coordinado como el punto elegido.
   • Un modelo lineal de la distribución de troncos de la propuesta resulta en un conjunto de distribuciones exponenciales en sentido parcial (es decir, segmentos de una o más distribuciones exponenciales, final a final). Las distribuciones exponenciales son bien comportadas y bien comprendidas. El logaritmo de una distribución exponencial es una línea recta, y por lo tanto este método implica esencialmente encerrar el logaritmo de la densidad en una serie de segmentos de línea. Esta es la fuente de la restricción log-concave: si una distribución es log-concave, entonces su logarithm es concave (en forma como un U al revés), lo que significa que un segmento de línea tangente a la curva siempre pasará sobre la curva.
   • Si no funciona en el espacio de registro, una función de densidad lineal en sentido parcial también puede ser muestreada a través de distribuciones triangulares
3. Podemos aprovechar aún más el requisito de concavidad (log) para evitar potencialmente el costo de evaluar ${displaystyle fleft(xright)}$ cuando su muestra es aceptado.
   • Al igual que podemos construir un borde superior lineal (la función "envelope") usando los valores de ${displaystyle hleft(xright)}$ que hemos tenido que evaluar en la cadena actual de rechazos, también podemos construir un límite inferior lineal (la función de "squeezing") a partir de estos valores también.
   • Antes de evaluar (el potencialmente caro) ${displaystyle fleft(xright)}$ para ver si su muestra será aceptada, podemos Ya lo sé. si será aceptado comparando con el (idealmente más barato) ${displaystyle g_{l}left(xright)}$ (o ${displaystyle h_{l}left(xright)}$ en este caso) función de squeezing que tienen disponible.
   • Este escalón es opcional, incluso cuando lo sugiere Gilks. En el mejor de los casos, le ahorra de sólo una evaluación adicional de su densidad de destino (messy/o costosa). Sin embargo, presumiblemente para funciones de densidad particularmente costosas (y asumiendo la rápida convergencia de la tasa de rechazo hacia cero) esto puede hacer una diferencia considerable en el tiempo de ejecución final.

El método consiste esencialmente en determinar sucesivamente una envolvente de segmentos de línea recta que se aproxima cada vez mejor al logaritmo sin dejar de permanecer por encima de la curva, comenzando con un número fijo de segmentos (posiblemente solo una línea tangente). El muestreo a partir de una variable aleatoria exponencial truncada es sencillo. Simplemente tome el registro de una variable aleatoria uniforme (con el intervalo apropiado y el truncamiento correspondiente).

Desafortunadamente, ARS solo se puede aplicar para el muestreo de densidades objetivo logarítmicas cóncavas. Por esta razón, en la literatura se han propuesto varias extensiones de ARS para abordar distribuciones objetivo no logarítmicas cóncavas. Además, se han diseñado diferentes combinaciones de ARS y el método Metropolis-Hastings para obtener un sampler universal que construya una propuesta de densidades autoajustable (es decir, una propuesta construida y adaptada automáticamente al objetivo). Esta clase de métodos a menudo se denomina algoritmos de muestreo de metrópolis de rechazo adaptativo (ARMS). Las técnicas adaptativas resultantes siempre se pueden aplicar, pero en este caso las muestras generadas están correlacionadas (aunque la correlación desaparece rápidamente a cero a medida que crece el número de iteraciones).