Generación de variantes aleatorias al azar doble Triangular distribuido

Question

Wikipedia muestra cómo generar Triangular-variables aleatorias distribuidas mediante una varia $U$ extraídas de la distribución uniforme.

Un "Doble Triangular" la distribución es un caso especial de una mezcla de dos distribuciones triangulares. Específicamente, se determina por tres números de $a \lt c \lt b$ y una proporción$p$$0 \lt p \lt 1$. Se apoya en el intervalo de $[a, b]$.

En el intervalo de $[a,c]$ su función de densidad está dada por

$$f(x) = \lambda(x-a)$$

donde $\lambda(c-a)^2 = 2p$.

En el intervalo de $[c, b]$ su función de densidad está dada por

$$f(x) = \mu(b-x)$$

donde $\mu(b-c)^2 = 2(1-p)$.

¿Cuál sería un buen enfoque para la generación de Doble Triangular-variables aleatorias distribuidas?

Answer 1

Siguiente whuber comentarios, aquí es una manera de simular a partir de la mezcla utilizando un único y uniforme:

Desde el pdf de esta mezcla es $$f(x)=p\dfrac{2(x-a)}{(c-a)^2}\mathbb{I}_{(a,c)}(x)+(1-p)\dfrac{2(b-x)}{(b-c)^2}\mathbb{I}_{(c,b)}(x)$$ el cdf es$$F(x)=p\dfrac{(x-a)^2}{(c-a)^2}\mathbb{I}_{(a,c)}(x)+\left\{p+(1-p)\left[1-\dfrac{(b-x)^2}{(b-c)^2}\right]\right\}\mathbb{I}_{(c,b)}(x)+\mathbb{I}_{(b,\infty)}(x)$$ Dibujo de un único y uniforme $U\sim\mathcal{U}(0,1)$ y la inversión de la ecuación$$u=p\dfrac{(x-a)^2}{(c-a)^2}\mathbb{I}_{(a,c)}(x)+\left\{p+(1-p)\left[1-\dfrac{(b-x)^2}{(b-c)^2}\right]\right\}\mathbb{I}_{(c,b)}(x)$$produces a generation of $X\sim f(x)$. Esta inversión puede ser descompuesto en dos pasos:

1. Determinar si $u\le p$ o $u>p$ [esta es la descomposición de los uniformes sugerido por whuber]
2. Si $u\le p$, solucionar $(x-a)^2=u(c-a)^2$, es decir, tomar $$x=a+\sqrt{u}(c-a)$$
3. Otra cosa, si $u>p$, $v=(u-p)/(1-p)$ y solucionar $(b-x)^2=(1-v)(b-c)^2$, es decir, tomar $$x=b-\sqrt{1-v}(b-c)$$or equivalently$$x=b-\sqrt{v}(b-c)$$

Curiosamente, la última ecuación es no equivalente a $$x=c+\sqrt{1-v}(b-c)$$porque de el cuadrado de la raíz.

Answer 2

El rechazo de muestreo es muy simple aquí (ver la aceptación o de respuesta al final de este post para alternativas más eficientes, a pesar de que):

1. Construir un rectángulo alrededor de la densidad (que es fácil aquí, debido a la limitada apoyo) de la caja verde en la figura.
2. Dibujar uniforme de variables en el cuadro (por lo que sólo dos uniformes con los vectores de la derecha apoye en ambas direcciones) - todos los puntos de la figura, tanto en blanco y rojo.
3. Mantener los puntos que están por debajo de la densidad de la doble triangular (los puntos rojos en la figura).
4. Tomar el "$x$-eje" de las coordenadas de estos puntos como realizaciones de la doble triangular.

En una imagen:

Creado con

a <- 0
b <- 2
c <- 1
lambda <- .5
x <- seq(a,b,by=.01)
ddoubletriangular <- function(x,a,b,c,lambda){
  p <- lambda/2*(c-a)^2
  mu <- 2*(1-p)/(b-c)^2
  return(lambda*(x-a)*(x<=c)+mu*(b-x)*(x>c))
}
plot(x,ddoubletriangular(x,a,b,c,lambda),type="l")
upperbound <- max(ddoubletriangular(x,a,b,c,lambda))
lines(c(b,b,a,a,b),c(0,upperbound,upperbound,0,0),col="green",lwd=4)

px <- runif(100,min=a,max=b)
py <- runif(100,min=0,max=upperbound)
points(px,py)

ikeep <- (py < ddoubletriangular(px,a,b,c,lambda))
px <- px[ikeep]
py <- py[ikeep]
points(px,py,col="red",pch=19)  

Vamos a tratar más que $N=100$ sorteos para el estudio de la calidad del algoritmo:

N <- 1500000
px <- runif(N,min=a,max=b)
py <- runif(N,min=0,max=upperbound)
ikeep <- (py < ddoubletriangular(px,a,b,c,lambda))
px <- px[ikeep]
truehist(px,main="Rejection method",sub=paste("n=",length(px),"of",N))
lines(x,ddoubletriangular(x,a,b,c,lambda),type="l")

Esto se traduce en

ACTUALIZACIÓN:

Fuerza documentado en los comentarios, el algoritmo anterior es, "fuera de la plataforma", innecesariamente ineficiente. Aquí es una ilustración de lo que @whuber propone (o eso espero):

Primera mejora:

La construcción de dos cajas rectangulares y realizar rechazo de muestreo en cada región por separado (con una probabilidad de $p$ para el primer segmento, y $1-p$ para el segundo), ver la siguiente figura. Esto mediante la construcción mejora la probabilidad de aceptación a 1/2.

Segunda mejora:

Reconocer que las regiones en las casillas de cada uno son lo @whuber llamadas "rectangular simétrica": el triángulo situado encima de un segmento de la densidad lineal es una versión rotada de la aceptación triángulo debajo de la densidad. Por lo tanto, podemos usar el "rechazado" muestras del algoritmo básico y el uso de estas como legítimo realizaciones para el otro extremo del segmento.

Por ejemplo, el rechazo de un sorteo con $x$-coordinar $c+\epsilon$ puede ser utilizado como un aceptados sorteo de $b-\epsilon$, debido a la "altura" de la región de rechazo en $c+\epsilon$ es por la construcción idéntica a la "altura" de la aceptación de la región en $b-\epsilon$. Esto tiene la interesante característica de aumento de la tasa de aceptación del 100%.

Aquí es una forma de código:

a <- 0
b <- 2
c <- 1
lambda <- .5
x <- seq(a,b,by=.002)
p <- lambda/2*(c-a)^2
mu <- 2*(1-p)/(b-c)^2
ddoubletriangular <- lambda*(x-a)*(x<=c)+mu*(b-x)*(x>c)

upperbound1 <- lambda*c
upperbound2 <- mu*x[min(which(x>c))+1]

N <- 1500000
px1 <- runif(N*p,min=a,max=c)
py1 <- runif(N*p,min=0,max=upperbound1)

px2 <- runif(N*(1-p),min=c,max=b)
py2 <- runif(N*(1-p),min=0,max=upperbound2)

ikeep1 <- (py1 < lambda*(px1-a))
px1a <- px1[ikeep1]
px1b <- a+c-px1[!ikeep1]

ikeep2 <- (py2 < mu*(b-px2))
px2a <- px2[ikeep2]
px2b <- c+b-px2[!ikeep2]

px <- c(px1a,px1b,px2a,px2b)

truehist(px,main="Modified rejection method",sub=paste("n=",length(px),"of",N))
lines(x,lambda*(x-a)*(x<=c)+mu*(b-x)*(x>c),type="l")