¿Cómo puedo muestrear de una distribución con CDF no computable?

Problema relacionado con la simulación de semi-informática aquí.

Tengo una distribución donde

P (x) =$\frac{(e^b-1) e^{b (n-x)}}{e^{b n+b}-1}$

para algunas constantes byn, yx es un número entero tal que .$0\leq x \leq n$

Ahora, necesito probar de esta distribución. Tiene un CDF invertible, por lo que es posible hacer esto directamente en teoría. El problema es que los números involucrados son GRANDES. Tan grande, de hecho, que ambos desbordan las variables formateadas convencionalmente, y toman al menos minutos (en algún momento me di por vencido ...) para calcular usando formatos de precisión arbitrarios. Básicamente, el CDF inverso todavía implica un término de , para . A pesar de esto, los números de salida seguirán en el rango , por lo que parece que debería haber una forma de hacerlo.$e^{b(n+1)}$$350 < n < 3500$$0-n$

Lo que estoy buscando es una forma de muestreo aproximadamente de esta distribución que sea computable. ¿Existen métodos alternativos de muestreo? ¿Qué son?

El CDF es fácilmente invertible. Una fórmula para la inversión conduce a lo que tiene que ser una de las soluciones más simples y convenientes posibles.

Comience observando que la probabilidad del resultado , , es proporcional a . Por lo tanto, si generamos un valor uniforme entre y = , solo necesitamos encontrar la más grande para la cual$k$$0 \le k \le n$$e^{-b k}$$q$$0$$q_{\max}=\sum_{k=0}^{n} e^{-b k}$$(1 - e^{-b(n+1)})/(1 - e^{-b})$$k$

Álgebra simple da la solución

Aquí hay una Rimplementación construida como todos los otros generadores de números aleatorios: su primer argumento especifica cuántos valores de iid generar y el resto de los argumentos nombran los parámetros ( as y as ):$b$b$n$n.max

rgeom.truncated <- function(n=1, b, n.max) {
  a <- 1 - exp(-b)
  q.max <- (1 - exp(-b*(n.max+1))) / a
  q <- runif(n, 0, q.max)
  return(-ceiling(log(1 - q*a) / b))
}

Como ejemplo de su uso, generemos un millón de variantes aleatorias de acuerdo con esta distribución:

b <- 0.001
n.max <- 3500
n.sim <- 10^6
set.seed(17)
system.time(sim <- rgeom.truncated(n.sim, b,n.max))

( Se necesitaron segundos).$0.10$

h <- hist(sim+1, probability=TRUE, breaks=50, xlab="Outcome+1")
pmf <- exp(-b * (0: n.max)); pmf <- pmf / sum(pmf)
lines(0:n.max, pmf, col="Red", lwd=2)

( Se agregó a cada valor para crear un mejor histograma: el procedimiento tiene una idiosincrasia (= error) en la que la primera barra es demasiado alta cuando el punto final izquierdo se establece en cero). La curva roja es la distribución de referencia que esta simulación intenta reproducir. Vamos a evaluar la bondad del ajuste con una prueba de chi-cuadrado:$1$Rhist

observed <- table(sim)
expected <- n.sim * pmf
chi.square <- (observed-expected)^2 / expected
pchisq(sum(chi.square), n.max, lower.tail=FALSE)

El valor p es : un ajuste hermoso.$0.84$

Se trata de una distribución geométrica truncada con . Hay una variedad de formas de abordar esto.$p = 1-e^{-b}$

Aconsejaría diferentes opciones en diferentes situaciones; Algunas opciones implicarían simular a partir de una geometría y regeneración cuando está fuera del rango, tomar la parte entera de un exponencial truncado apropiado ( como aquí ), o usar cualquiera de varias técnicas rápidas adaptadas a distribuciones discretas en un rango finito. Dado que es grande, tomar el piso de un exponencial truncado probablemente sea relativamente rápido, pero si es la mejor opción también depende de .$n$$b$

Aquí hay una pregunta relacionada con las matemáticas.

Antes de intentar sugerencias específicas, ¿cuál es un rango típico de valores para ?$b$

Primero, tenga en cuenta que que, si fuera continuo, estaría relacionado con una distribución exponencial. Entonces, lo que puede hacer es simular a partir de una distribución exponencial truncada y tomar la (parte entera) de las observaciones.$P(x)\propto e^{-bx}$$x$floor()

El cdf de un exponencial truncado es

Entonces, si hacemos , obtenemos que . Si es grande, entonces que sugiere aproximar .$F(x;n,b)=u$$x=-\dfrac{1}{b}\log[1-u(1-e^{-bn})]$$bn$$e^{-bn}\approx 0$$x\approx -\dfrac{1}{b}\log[1-u]$

rweirdp <- function(ns,n,b){
u <- runif(ns)
samp <- - log(1-u*(1-exp(-n*b)))/b
return(floor(samp))
}

rweirdp(1000,10,1)

Una forma de tomar muestras de la distribución objetivo $p(k)\propto \exp\{-bk\}$ Es para

1. ejecutar un experimento Metropolis-Hastings para determinar el soporte (interesante) de la distribución, es decir, en qué subconjunto de $\{0,1,\ldots,n\}$ se concentra;

   metro=function(N,b,n){
   x=sample(0:n,N,rep=TRUE)
   for (t in 2:N){
     x[t]=prop=x[t-1]+sample(c(-1,1),1)
   
     if ((prop<0)||(prop>n)||(log(runif(1))>b*(x[t]-prop)))
         x[t]=x[t-1]
     }
   return(x)
   }
   

2. Utilice el soporte así determinado, $\{k_0,\ldots,k_1\}$ decir, para calcular las probabilidades exactas como $p(k)\propto \exp\{-bk+bk_0\}$ para evitar desbordamientos.

Actualización: al pensar más en ello, ya que$p(\cdot)$ está disminuyendo en k, el soporte efectivo de la distribución siempre comenzará en $k_0=0$. Si$b$ es bastante grande, este soporte terminará muy rápidamente, en cuyo caso $n$ no importa tanto como valores grandes de $k$Nunca será visitado. Si$b$ es muy pequeño, el pdf es casi plano, lo que significa que se puede usar una distribución uniforme en $\{0,1,\ldots,n\}$como una propuesta de aceptar-rechazar. Y use registros en el paso de aceptación para evitar desbordamientos.