Muestreo CDF inverso para una distribución mixta

La versión corta fuera de contexto

Sea una variable aleatoria con CDF $y$

F (\cdot) \equiv {\begin{cases} θ & y = 0 \\ θ + (1 - θ) \times {CDF}_{log-normal} (\cdot; μ, σ) & y > 0 \end{cases}

$F(\cdot) \equiv \cases{\theta & y = 0 \\ \theta + (1-\theta) \times \text{CDF}_{\text{log-normal}}(\cdot; \mu, \sigma) & y > 0}$

Digamos que quería simular dibujos de usando el método inverso de CDF. ¿Es eso posible? Esta función no tiene exactamente una inversa. Por otra parte, hay un muestreo de transformación inversa para la distribución de mezclas de dos distribuciones normales que sugiere que hay una forma conocida de aplicar el muestreo de transformación inversa aquí. $y$

Soy consciente del método de dos pasos, pero no sé cómo aplicarlo a mi situación (ver más abajo).

La versión larga con fondo

Ajusté el siguiente modelo para una respuesta con valor vectorial, , usando MCMC (específicamente, Stan): $y^i = \left( y_1 , \dots , y_K \right)^i$

θ_{k}^{i} \equiv {logit}^{- 1} (α_{k} x^{i}), μ_{k}^{i} \equiv β_{k} x^{i} - \frac{σ_{k}^{2}}{2} F (\cdot) \equiv {\begin{cases} θ & y = 0 \\ θ + (1 - θ) \times {CDF}_{log-normal} (\cdot; μ, σ) & y > 0 \end{cases} u_{k} \equiv F (y_{k}), z_{k} \equiv Φ^{- 1} (u_{k}) z \sim N (0, R) \times \prod_{k} f (y_{k}) (α, β, σ, R) \sim priors

$\theta_k^i \equiv \operatorname{logit}^{-1}\left( \alpha_k x^i \right), \quad \mu_k^i \equiv \beta_k x^i - \frac{ \sigma^2_k }{ 2 } \\ F(\cdot) \equiv \cases{\theta & y = 0 \\ \theta + (1-\theta) \times \text{CDF}_{\text{log-normal}}(\cdot; \mu, \sigma) & y > 0} \\ u_k \equiv F(y_k), \quad z_k \equiv\Phi^{-1}{\left( u_k \right)} \\ z \sim \mathcal{N}(\mathbf{0}, R) \times \prod_k f(y_k) \\ \left( \alpha, \beta, \sigma, R \right) \sim \text{priors}$

donde indexa observaciones, es una matriz de correlación es un vector de predictores / regresores / características. $i$ $N$ $R$ $x$

Es decir, mi modelo es un modelo de regresión en el que se supone que la distribución condicional de la respuesta es una cópula gaussiana con márgenes log-normales inflados a cero. He publicado sobre este modelo antes; Resulta que Song, Li y Yuan (2009, cerrada ) lo han desarrollado y lo llaman un vector GLM o VGLM. La siguiente es su especificación tan cercana a textual como pude obtenerla: Mi

f (y; μ, φ, Γ) = c {G_{1} (y_{1}), \dots, G_{m} (y_{m}) | Γ} \prod_{i = 1}^{m} g (y_{i}; μ_{i}, φ_{i}) c (u | Γ) = {| Γ |}^{- 1 / 2} \exp (\frac{1}{2} q^{T} (I_{m} - Γ^{- 1}) q) q = {(q_{1}, \dots, q_{m})}^{T}, q_{i} = Φ^{- 1} (u_{i})

$f(\mathbf{y}; \mathbf{\mu}, \mathbf{\varphi}, \Gamma) = c\{ G_1(y_1), \dots, G_m(y_m) | \Gamma \} \prod_{i=1}^m g(y_i; \mu_i, \varphi_i) \\ c(\mathbf{u} | \Gamma) = \left| \Gamma \right|^{-1/2}\exp\left( \frac{1}{2} \mathbf{q}^T \left( I_m - \Gamma^{-1} \right) \mathbf{q} \right) \\ \mathbf{q} = \left( q_1, \dots, q_m \right)^T, \quad q_i = \Phi^{-1}(u_i)$

F_{K}

$F_K$ corresponde a su , mi corresponde a su , y mi corresponde a su ; los detalles se encuentran en la página 62 (página 3 del archivo PDF) pero son idénticos a lo que escribí aquí.

G_{m}

$G_m$

z

$z$

q

$\mathbf{q}$

R

$R$

Γ

$\Gamma$

La parte inflada a cero sigue aproximadamente las especificaciones de Liu y Chan (2010, sin delegar ).

Ahora me gustaría simular datos de los parámetros estimados, pero estoy un poco confundido sobre cómo hacerlo. Primero pensé que podía simular directamente (en código R): $y$

for (i in 1:N) {
    for (k in 1:K) {
        Y_hat <- rbinom(1, 1, 1 - theta[i, k])
        if (Y_hat == 1)
            Y_hat <- rlnorm(1, mu[i, k], sigma[k])
    }
}

que no usa en absoluto. Me gustaría intentar usar la matriz de correlación que estimé. $R$

Mi siguiente idea fue tomar sorteos de y luego convertirlos nuevamente a . Esto también parece coincidir con las respuestas en Generar muestras de Cópula en R y Muestreo bivariado para distribución expresadas en el teorema de cópula de Sklar? . Pero, ¿qué diablos es mi aquí? El muestreo de transformación inversa para la distribución de mezclas de dos distribuciones normales hace que parezca posible, pero no tengo idea de cómo hacerlo. $z$ $y$ $F^{-1}$

r distributions sampling simulation copula Shadowtalker
fuente

@ Xi'an es una cópula gaussiana, para estimar la dependencia entre los componentes .

y

$y$

shadowtalker

El hilo al que hace referencia sobre el muestreo de mezclas de normales se aplica directamente a su problema sin ninguna modificación esencial: en lugar de usar los CDF inversos de las normales, use los CDF inversos de sus dos componentes. El CDF inverso del átomo en es la función constante, siempre igual a .

y = 0

$y=0$

0

$0$

whuber

@whuber Estoy confundido acerca de cómo usar los CDF inversos de los dos componentes: ¿qué dibujo, de qué lo saco y luego en qué enchufo cada cosa?

shadowtalker

@ Xi'an explica muy bien eso en su respuesta a la pregunta de la mezcla normal: utiliza una variante uniforme para seleccionar el componente de la mezcla y luego extrae un valor de ese componente (de la forma que desee). En su caso, es excepcionalmente fácil extraer un valor del primer componente: ¡siempre es ! Para dibujar un valor del segundo componente, use cualquier generador de números aleatorios lognormal que desee. En cada caso, terminas con un número: no hay que "enchufarlo" para lograrlo; El objetivo de la generación de números aleatorios es obtener ese número.

0

$0$

whuber

@whuber la nueva respuesta me lo aclaró. Gracias a los dos.

shadowtalker

Respuestas:

La respuesta a la versión larga con antecedentes:

Esta respuesta a la versión larga de alguna manera aborda otro problema y, dado que parece que tenemos dificultades para formular el modelo y el problema, elijo reformularlo aquí, con suerte correctamente.

Para , el objetivo es simular los vectores modo que, condicional a una covariable , con . Por lo tanto, si uno quiere simular datos de este modelo, podría proceder de la siguiente manera: $1\le i\le I$ $y^i=(y^i_1,\ldots,y^i_K)$ $x^i$

y_{k}^{i} = {\begin{cases} 0 & with probability {logit}^{- 1} (α_{k} x^{i}) \\ \log (σ_{k} z_{k}^{i} + β_{k} x^{i}) & with probability 1 - {logit}^{- 1} (α_{k} x^{i}) \end{cases}

$y_k^i = \begin{cases} 0 &\text{ with probability }\operatorname{logit}^{-1}\left( \alpha_k x^i \right)\\ \log(\sigma_k z_k^i + \beta_k x^i) &\text{ with probability }1-\operatorname{logit}^{-1}\left( \alpha_k x^i \right)\end{cases}$

z^{i} = (z_{1}^{i}, \dots, z_{K}^{i}) \sim N_{K} (0, R)

$z^i=(z^i_1,\ldots,z^i_K)\sim\mathcal{N}_K(0,R)$

Para , $1\le i\le I$

Generar $z^i=(z^i_1,\ldots,z^i_K)\sim\mathcal{N}_K(0,R)$
Genere $u^i_1,\ldots,u^i_K \stackrel{\text{iid}}{\sim} \mathcal{U}(0,1)$
Derive para $y^i_k=\mathbb{I}\left\{u^i_k>\operatorname{logit}^{-1}\left( \alpha_k x^i \right)\right\}\, \log\{ \sigma_k z_k^i + \beta_k x^i\}$ $1\le k\le K$

Si uno está interesado en la generación desde la parte posterior de dada la , este es un problema más difícil, aunque factible por muestreo de Gibbs o ABC. $(\alpha,\beta,\mu,\sigma,R)$ $y^i_ k$

Xi'an
fuente

Sabía que me faltaba algo. "Todo es obvio en retrospectiva". Mi intención: estoy interesado en el valor de , así que sí, estoy interesado en dibujar desde la articulación posterior de los parámetros. Quiero las simuladas para ver si el modelo se ajusta.

F (y^{i} | x^{i})

$F(y^i | x^i)$

y

$y$

shadowtalker

¿Cómo es el segundo problema mucho más difícil? Ya he estimado el modelo y tengo dibujos posteriores. Podemos continuar en el chat si lo desea, para evitar abarrotar los comentarios aquí.

shadowtalker

Oh, en general, si. Afortunadamente tengo a Stan y al Sampler No-U-Turn haciendo el trabajo duro por mí allí.

shadowtalker

La respuesta a la versión corta fuera de contexto:

"Invertir" un cdf que no es invertible en el sentido matemático (como su distribución mixta) es factible, como se describe en la mayoría de los libros de texto de Monte Carlo. (Al igual que el nuestro , vea el Lema 2.4.) Si define el inverso generalizado entonces Esto significa que, cuando tiene un salto de en , para . En otras palabras, si dibujas un uniforme y termina más pequeño que , tu generación de

F^{-} (u) = inf {x \in R; F (x) \geq u}

$F^{-}(u) = \inf\left\{ x\in\mathbb{R};\ F(x)\ge u \right\}$

X \sim F is equivalent to X = F^{-} (U) when U \sim U (0, 1) .

$X \sim F \text{ is equivalent to } X=F^{-}(U)\text{ when } U\sim\mathcal{U}(0,1)\,.$

F (y)

$F(y)$

θ

$\theta$

y = 0

$y=0$

F^{-} (u) = 0

$F^{-}(u)=0$

u \leq θ

$u\le\theta$

U (0, 1)

$\mathcal{U}(0,1)$

θ

$\theta$

X

$X$ es . De lo contrario, cuando , terminas generando a partir de la parte continua, es decir, el log-normal en tu caso. Esto significa usar una segunda generación aleatoria uniforme, , independiente del primer sorteo uniforme y establecer para obtener una generación logarítmica normal.

x = 0

$x=0$

u > θ

$u>\theta$

v

$v$

y = \exp (μ + σ Φ^{- 1} (v))

$y=\exp(\mu+\sigma\Phi^{-1}(v))$

Esto es casi lo que su código R

Y_hat <- rbinom(1, 1, theta[i, k]) if (Y_hat == 1) Y_hat <- rlnorm(1, mu[i, k], sigma[k])

está haciendo. un Bernoulli con probabilidad y si es igual a , lo conviertes en un log-normal. Dado que es igual a 1 con probabilidad , en su lugar, debe convertirlo en una simulación logarítmica normal cuando sea igual a cero , terminando con el código R modificado: $\theta_k^i$ $1$ $\theta_k^i$

Y_hat <- rbinom(1, 1, theta[i, k])
    if (Y_hat == 0)
        Y_hat <- rlnorm(1, mu[i, k], sigma[k])

Xi'an
fuente

Entonces, en conjunto, mi procedimiento de simulación sería: 1) dibujar , 2) calcular , luego 3) calcular si y contrario. ¿Correcto?

z

$z$

u_{k} = Φ (z_{k})

$u_k = \Phi(z_k)$

y_{k} = 0

$y_k = 0$

u_{k} \leq θ

$u_k \leq \theta$

y_{k} = F_{log-normal}^{- 1} (u_{k})

$y_k = F^{-1}_\text{log-normal}(u_k)$

shadowtalker

No, incorrecto Dibuja un primer uniforme para decidir entre y log-normal, luego un segundo uniforme en caso de que haya decidido un log-normal. Vea la versión editada de mi respuesta.

0

$0$

Xi'an

Pero eso ignora el componente ; De ahí mi pregunta. Hice una edición aclaratoria y también abordé el error en mi pseudocódigo.

z

$z$

shadowtalker

Mi respuesta es para la versión corta y para el código R que proporcionó. Espero que ayude para la versión larga, pero su fórmula para el modelo conjunto sigue siendo incorrecta. Debe definir el modelo en las sin usar los uniformes ...

y

$y$

Xi'an

¿Cómo es incorrecto ese modelo? Acabo de conectar mi en la fórmula proporcionada por el documento que cité (correspondiente a en su notación). ¿Es eso inválido?

F_{1}, \dots, F_{K}

$F_1, \dots, F_K$

G_{1}, \dots, G_{m}

$G_1, \dots, G_m$

shadowtalker