Muestreo eficiente de una distribución Beta umbralizada

10

¿Cómo debo muestrear eficientemente de la siguiente distribución?

xB(α,β), x>k

Si no es demasiado grande, entonces el muestreo de rechazo puede ser el mejor enfoque, pero no estoy seguro de cómo proceder cuando es grande. ¿Quizás haya alguna aproximación asintótica que se pueda aplicar?kkk

user1502040
fuente
1
No está inequívocamente claro lo que pretende allí con " ". ¿Te refieres a una distribución beta truncada (truncada a la izquierda en )? xB(α,β), x>kk
Glen_b -Reinstate Monica el
@Glen_b exactamente.
user1502040
55
Para ambos parámetros de forma mayor que 1, la distribución beta es cóncava logarítmica, por lo que se pueden usar envolventes exponenciales para el muestreo de rechazo. Para generar variantes beta no truncadas que ya está tomando muestras de distribuciones exponenciales truncadas (que es fácil de hacer), debería ser sencillo adaptar este método.
Scortchi - Restablece a Monica

Respuestas:

14

La forma más simple y más general, que se aplica a cualquier distribución truncada (también se puede generalizar al truncamiento en ambos lados), es usar muestreo de transformación inversa . Si es la distribución acumulativa de interés, establezca y tomep 0 = F ( k )Fp0=F(k)

UU(p0,1)X=F1(U)

donde es una muestra de truncada a la izquierda en . La función cuantil trazará un mapa de probabilidades a partir de muestras . Dado que tomamos valores de solo del "área" que coincide con los valores de distribución beta de la región no truncada, solo tomará muestras de esos valores.F k F - 1 F UXFkF1FU

Este método se ilustra en la imagen a continuación, donde el área truncada se marca con un rectángulo gris, los puntos en rojo se dibujan desde la distribución y luego se transforman en muestrasB ( 2 , 8 )U(p0,1)B(2,8)

Transformación inversa de muestreo de distribución truncada

Tim
fuente
55
(+1) Vale la pena señalar que la función cuantil no se evalúa tan fácilmente.
Scortchi - Restablece a Monica
1
@Scortchi Si aob son 1 o al menos un número entero, hay una forma no tan mala (ver wikipedia ). Y en Python hay scipy.special.betaincpara el inverso y en R hay pbeta.
Graipher
3
@Graipher: Debería haber dicho "a bajo precio, en general", sería mejor evitar Newton-Raphson u otras soluciones iterativas si es posible. (Por cierto, es qbetapara la función cuantil en R.)
Scortchi - Restablecer Monica
1
@Scortchi tiene razón, pero en la mayoría de los casos, para las computadoras modernas esto no debería ser un problema importante. También recomiendo este enfoque ya que está disponible directamente en la mayoría de los programas y puede generalizarse a cualquier distribución truncada, solo si uno tiene acceso a la función cuantil.
Tim
1
Indudablemente, es bueno tener a mano un método generalmente aplicable y fácil de implementar cuyo tiempo de ejecución no crezca con ; & para distribuciones con funciones cuantiles de forma cerrada, por ejemplo, Weibull, debe ser tan bueno como sea posible. Sin embargo, sospecho que tendrá que configurarse para cortar una parte bastante grande de la distribución beta antes de que supere los algoritmos eficientes de muestreo de rechazo que también están disponibles en la mayoría del software y que se basan solo en el cálculo de la densidad de probabilidad de la beta. kk
Scortchi - Restablece a Monica
8

La respuesta de @ Tim muestra cómo el muestreo de transformación inversa puede adaptarse para distribuciones truncadas, liberando el tiempo de ejecución de la dependencia del umbral . Se pueden obtener mayores eficiencias evitando la costosa evaluación numérica de la función de cuantiles beta y utilizando el muestreo de transformación inversa como parte del muestreo de rechazo.k

La función de densidad de una distribución beta con parámetros de forma y doblemente truncada en (para un poco más de generalidad) esαβk1<k2

f(x)=x(α1)(1x)(β1)B(k2,α,β)B(k1,α,β)

Tome cualquier parte monótonamente creciente de la densidad entre y : para es log-cóncavo, por lo que puede envolverlo con una función exponencial dibujada en una tangente a cualquier punto a lo largo de ella:xLxUα,β>1

g(x)=cλeλ(xxL)

Encuentre estableciendo los gradientes de las densidades de registro igualesλ

λ=a1xb11x
y encuentre calculando cuánto necesita aumentar la densidad exponencial para cumplir con la densidad en ese punto c
c=f(x)λeλ(xxL)

ingrese la descripción de la imagen aquí

El mejor sobre para el muestreo de rechazo es el que tiene el área más pequeña debajo. El área es Sustituyendo expresiones en por & , & dropping factores constantes da

A=c(1eλ(xUxL))
xλc

Q(x)=xa(1x)b(a+b2)xa+1[exp((b1)(xxL)1x+xL(a1)x(a1))exp((b1)(xxU)1x+xU(a1)x(a1))]

Escribir la derivada se deja como un ejercicio para los lectores o sus computadoras. Cualquier algoritmo de búsqueda de raíz se puede utilizar para encontrar la para la cual . xdQdQdxxdQdx=0

Lo mismo ocurre, mutatis mutandis, para partes de la densidad que disminuyen monotónicamente. Por lo tanto, la distribución beta truncada puede estar perfectamente envuelta por dos funciones exponenciales, por ejemplo, una de al modo y otra del modo a , listas para el muestreo de rechazo. (Para una variable aleatoria uniforme truncada , tiene una distribución exponencial truncada con el parámetro de velocidad .)k 2 U - log ( 1 - U )k1k2U λlog(1U)λλ

ingrese la descripción de la imagen aquí

La belleza de este enfoque es que todo el trabajo duro está en la configuración. Una vez que se define la función envolvente, se calcula la constante de normalización para la densidad beta truncada, todo lo que queda es generar variables aleatorias uniformes y realizar en ellas algunas operaciones aritméticas simples, registros y potencias, y comparaciones. Ajustar la función de envolvente —con líneas horizontales o curvas más exponenciales— puede, por supuesto, reducir el número de rechazos.

Scortchi - Restablece a Monica
fuente
1
+1 Buena idea. Debido a que Beta es aproximadamente normal para valores modestos a grandes de sus parámetros, dependiendo de qué tan cerca estén unos de otros, el uso de una envoltura gaussiana podría ser bastante más eficiente todavía.
whuber
@whuber: quería una distribución con una distribución cuantil de forma cerrada para la envolvente, pero supongo que podría generar las variables gaussianas truncadas de manera eficiente utilizando una de esas buenas aproximaciones a la función cuantil gaussiana. Todavía estoy interesado en lo que harías cuando o . β < 1α<1β<1
Scortchi - Restablece a Monica
1
Para estos pequeños y , querrás cambiar a una cola exponencial. Sin embargo, no estoy seguro de lo que quieres decir con "forma cerrada": cuando miras detenidamente las implementaciones informáticas de exponenciales, gaussianas, funciones hipergeométricas, etc., es decir, todas las funciones no algebraicas, descubres que ninguna de las tienen una "forma cerrada" general: se calculan a través de aproximaciones sucesivas como series de Taylor, fracciones parciales o expansiones asintóticas. No hay mucha diferencia entre calcular un exponencial gaussiano y calcular un exponencial gaussiano. βαβ
whuber
@whuber: (1) El enfoque que he tomado aquí para construir sobres no funcionaría porque las densidades no son cóncavas logarítmicas. (2) (a) Me refería ciertamente a funciones algebraicas + registros y potencias, trig. funciones si me hubieran preguntado, y tal vez incluso funciones gamma. Confieso que no tenía una noción precisa. (b) Punto tomado: las evaluaciones rápidas de funciones no se limitan a aquellas con formas cerradas.
Scortchi - Restablece a Monica
1
Buen punto sobre el fracaso de la concavidad logarítmica. Sospecho que una distribución de la ley de potencia debería ser un buen sobre para o . β < 1α<1β<1
whuber