Odwrotne próbkowanie CDF dla mieszanego rozkładu

Krótka wersja poza kontekstem

Pozwolić $y$ być zmienną losową z CDF

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

Powiedzmy, że chciałem symulować losowania $y$przy użyciu odwrotnej metody CDF. Czy to jest możliwe? Ta funkcja nie ma dokładnie odwrotności. Z drugiej strony istnieje próbkowanie z transformacją odwrotną dla rozkładu mieszanki dwóch rozkładów normalnych, co sugeruje, że istnieje znany sposób zastosowania próbkowania z transformacją odwrotną.

Znam metodę dwuetapową, ale nie wiem, jak zastosować ją w mojej sytuacji (patrz poniżej).

---

Długa wersja z tłem

Dopasowałem następujący model do odpowiedzi o wartości wektorowej, $y^i = \left( y_1 , \dots , y_K \right)^i$, używając MCMC (konkretnie Stan):

$$\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}$$

gdzie $i$ indeksy $N$ obserwacje, $R$ jest macierzą korelacji, oraz $x$ jest wektorem predyktorów / regresorów / funkcji.

Oznacza to, że mój model jest modelem regresyjnym, w którym zakłada się, że warunkowy rozkład odpowiedzi jest kopulą Gaussa z zerowymi napełnieniami logarytmiczno-normalnymi. Wcześniej pisałem o tym modelu; okazuje się, że Song, Li i Yuan (2009, gated ) opracowali go i nazywają go wektorem GLM lub VGLM. Poniżej podano ich specyfikację tak bliską dosłownym, jak to tylko możliwe:

$$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)$$ Mój $F_K$ odpowiada ich $G_m$mój $z$ odpowiada ich $\mathbf{q}$, i mój $R$ odpowiada ich $\Gamma$; szczegóły znajdują się na stronie 62 (strona 3 pliku PDF), ale poza tym są identyczne z tym, co tu napisałem.

Część z napompowaniem zerowym z grubsza jest zgodna ze specyfikacją Liu i Chana (2010 r., Niestrzeżona ).

Teraz chciałbym zasymulować dane z szacowanych parametrów, ale jestem trochę zdezorientowany, jak sobie z tym poradzić. Najpierw pomyślałem, że mogę po prostu symulować$y$ bezpośrednio (w kodzie R):

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])
    }
}

który nie używa $R$w ogóle. Chciałbym spróbować użyć oszacowanej macierzy korelacji.

Mój następny pomysł polegał na losowaniu $z$ a następnie przekonwertować je z powrotem na $y$. Wydaje się, że jest to również zbieżne z odpowiedziami w Generowanie próbek z Copuli w R i Bivariate próbkowania do dystrybucji wyrażonych w twierdzeniu o kopule Sklara? . Ale co do cholery jest moje$F^{-1}$tutaj? Odwrotne próbkowanie transformacji dla rozkładu mieszanki dwóch rozkładów normalnych sprawia, że ​​brzmi to tak, jak to jest możliwe, ale nie mam pojęcia, jak to zrobić.

Odpowiedź na długą wersję z tłem:

Ta odpowiedź na długą wersję nieco rozwiązuje inny problem, a ponieważ wydaje się, że mamy trudności z sformułowaniem modelu i problemu, postanowiłem sformułować go tutaj, mam nadzieję, że poprawnie.

Dla $1\le i\le I$, celem jest symulacja wektorów $y^i=(y^i_1,\ldots,y^i_K)$ takie, które zależą od współzmiennej $x^i$,

$$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 $z^i=(z^i_1,\ldots,z^i_K)\sim\mathcal{N}_K(0,R)$. Dlatego jeśli chcemy symulować dane z tego modelu, można postępować w następujący sposób:

Dla $1\le i\le I$,

1. Generować $z^i=(z^i_1,\ldots,z^i_K)\sim\mathcal{N}_K(0,R)$
2. Generować $u^i_1,\ldots,u^i_K \stackrel{\text{iid}}{\sim} \mathcal{U}(0,1)$
3. Czerpać $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\}$ dla $1\le k\le K$

Jeśli ktoś jest zainteresowany pokoleniem od tyłu $(\alpha,\beta,\mu,\sigma,R)$ biorąc pod uwagę $y^i_ k$, jest to trudniejszy problem, aczkolwiek wykonalny przez próbkowanie Gibbs lub ABC.

Odpowiedź na krótką wersję poza kontekstem:

„Odwracanie” pliku cdf, który nie jest odwracalny w sensie matematycznym (podobnie jak dystrybucja mieszana), jest wykonalne, jak opisano w większości podręczników Monte Carlo. (Podobnie jak nasz , patrz Lemat 2.4.) Jeśli zdefiniujesz uogólnioną odwrotność

$$F^{-}(u) = \inf\left\{ x\in\mathbb{R};\ F(x)\ge u \right\}$$ następnie $$X \sim F \text{ is equivalent to } X=F^{-}(U)\text{ when } U\sim\mathcal{U}(0,1)\,.$$ Oznacza to, że kiedy $F(y)$ ma skok $\theta$ w $y=0$, $F^{-}(u)=0$ dla $u\le\theta$. Innymi słowy, jeśli narysujesz mundur$\mathcal{U}(0,1)$ i kończy się mniejszy niż $\theta$, wasze pokolenie $X$ jest $x=0$. W przeciwnym razie$u>\theta$, w końcu generujesz z części ciągłej, a mianowicie log-normal w twoim przypadku. Oznacza to użycie drugiej jednolitej generacji losowej,$v$, niezależnie od pierwszego jednolitego losowania i ustawienia $y=\exp(\mu+\sigma\Phi^{-1}(v))$ w celu uzyskania log-normalnej generacji.

To prawie to, co twój kod R.

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

to robi. Generujesz Bernoulliego z prawdopodobieństwem$\theta_k^i$ a jeśli to jest równe $1$, zamieniasz go w log-normalny. Ponieważ jest to równe 1 z prawdopodobieństwem$\theta_k^i$powinieneś zamiast tego przekształcić go w logarytmiczną symulację, gdy jest ona równa zero , kończąc na zmodyfikowanym kodzie R:

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