Odwrotne próbkowanie CDF dla mieszanego rozkładu

9

Krótka wersja poza kontekstem

Pozwolić y być zmienną losową z CDF

F(){θ y = 0 θ+(1θ)×CDFlog-normal(;μ,σ) y > 0

Powiedzmy, że chciałem symulować losowania yprzy 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, yja=(y1,,yK.)ja, używając MCMC (konkretnie Stan):

θkjalogit-1(αkxja),μkjaβkxja-σk2)2)fa(){θ y = 0 θ+(1-θ)×CDFlog-normal(;μ,σ) y> 0ukfa(yk),zkΦ-1(uk)zN.(0,R)×kfa(yk)(α,β,σ,R)priors

gdzie ja 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(y;μ,φ,Γ)=c{G1(y1),,Gm(ym)|Γ}i=1mg(yi;μi,φi)c(u|Γ)=|Γ|1/2exp(12qT(ImΓ1)q)q=(q1,,qm)T.,qja=Φ-1(uja)
Mój faK. odpowiada ich solmmój z odpowiada ich q, i mój R odpowiada ich Γ; 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 Rw 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 mojefa-1tutaj? 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ć.

Shadowtalker
źródło
@ Xi'an to kopuła Gaussa, służąca do oszacowania zależności między yskładniki.
shadowtalker
1
Wątek, do którego odwołujesz się do próbkowania z mieszanin normalnych, odnosi się bezpośrednio do twojego problemu bez istotnych modyfikacji: zamiast używać odwrotnych CDF normalnych, użyj odwrotnych CDF twoich dwóch składników. Odwrotny CDF atomu przyy=0 jest funkcją stałą, zawsze równą 0.
whuber
@ whuber Jestem po prostu zdezorientowany, jak korzystać z odwrotnych CDF dwóch komponentów: co rysuję, z czego czerpię, a następnie do czego podłączam każdą rzecz?
shadowtalker,
1
@ Xi'an ładnie wyjaśnia, że ​​w swojej odpowiedzi na pytanie o mieszankę normalną: używasz zmiennej uniform, aby wybrać składnik mieszanki, a następnie rysujesz wartość z tego składnika (w dowolny sposób). W twoim przypadku wyjątkowo łatwo jest wyciągnąć wartość z pierwszego komponentu: zawsze0! Aby narysować wartość z drugiego komponentu, użyj dowolnego lognormalnego generatora liczb losowych, który ci się podoba. W każdym przypadku kończy się liczba: nie ma „podłączania” do osiągnięcia; głównym celem generowania liczb losowych jest uzyskanie tej liczby.
whuber
@ Whuber nowa odpowiedź wyjaśniła mi to. Dziękuję wam obu.
shadowtalker,

Odpowiedzi:

5

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 1jaja, celem jest symulacja wektorów yja=(y1ja,,yK.ja) takie, które zależą od współzmiennej xja,

ykja={0 z prawdopodobieństwem logit-1(αkxja)log(σkzkja+βkxja) z prawdopodobieństwem 1-logit-1(αkxja)
z zja=(z1ja,,zK.ja)N.K.(0,R). Dlatego jeśli chcemy symulować dane z tego modelu, można postępować w następujący sposób:

Dla 1jaja,

  1. Generować zja=(z1ja,,zK.ja)N.K.(0,R)
  2. Generować u1ja,,uK.jaiidU(0,1)
  3. Czerpać ykja=ja{ukja>logit-1(αkxja)}log{σkzkja+βkxja} dla 1kK.

Jeśli ktoś jest zainteresowany pokoleniem od tyłu (α,β,μ,σ,R) biorąc pod uwagę ykja, jest to trudniejszy problem, aczkolwiek wykonalny przez próbkowanie Gibbs lub ABC.

Xi'an
źródło
1
Wiedziałem, że czegoś mi brakuje. „Z perspektywy czasu wszystko jest oczywiste”. Moja intencja: interesuje mnie wartośćfa(yja|xja), więc tak, jestem zainteresowany rysowaniem ze wspólnego tylnego parametru. Chcę symulowaćySprawdź, czy model pasuje.
shadowtalker
1
W jaki sposób drugi problem jest znacznie trudniejszy? Oszacowałem już model i mam tylne rysunki. Jeśli chcesz, możemy kontynuować czat, aby nie zaśmiecać tutaj komentarzy.
shadowtalker,
1
Och, ogólnie tak. Na szczęście mam tam Stana i No-U-Turn Samplera, którzy wykonali dla mnie ciężką pracę.
shadowtalker,
7

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ść

fa-(u)=inf{xR; fa(x)u}
następnie
Xfa jest równa X=fa-(U) kiedy UU(0,1).
Oznacza to, że kiedy fa(y) ma skok θ w y=0, fa-(u)=0 dla uθ. Innymi słowy, jeśli narysujesz mundurU(0,1) i kończy się mniejszy niż θ, wasze pokolenie X jest x=0. W przeciwnym razieu>θ, 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(μ+σΦ-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θkja a jeśli to jest równe 1, zamieniasz go w log-normalny. Ponieważ jest to równe 1 z prawdopodobieństwemθkjapowinieneś 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])
Xi'an
źródło
Tak więc razem moją procedurą symulacyjną byłoby: 1) losowanie z, 2) oblicz uk=Φ(zk), a następnie 3) oblicz yk=0 gdyby ukθ i yk=falog-normal-1(uk)Inaczej. Poprawny?
shadowtalker
Nie, niepoprawnie. Najpierw rysujesz pierwszy mundur0i log-normal, a następnie drugi mundur na wypadek, gdybyś zdecydował się na log-normal. Zobacz zmodyfikowaną wersję mojej odpowiedzi.
Xi'an,
Ale to ignoruje zskładnik; stąd moje pytanie. Dokonałem edycji wyjaśniającej, a także rozwiązałem błąd w moim pseudokodzie.
shadowtalker
Moja odpowiedź dotyczy krótkiej wersji i podanego kodu R. Mam nadzieję, że to pomaga w długiej wersji, ale twoja formuła dla modelu połączenia jest nadal niepoprawna. Powinieneś zdefiniować model wybez użycia mundurów ...
Xi'an
Jak ten model jest nieprawidłowy? Właśnie podłączyłem mójfa1,,faK. do wzoru podanego w cytowanym przeze mnie artykule (odpowiadającym sol1,,solmw ich notacji). Czy to nieważne?
shadowtalker