Ile stron ma kość? Wnioskowanie bayesowskie w JAGS

Problem

Chciałbym wyciągnąć wnioski na temat systemu analogicznego do śmierci z nieznaną liczbą stron. Kostka jest walcowana kilka razy, po czym chciałbym wywnioskować rozkład prawdopodobieństwa dla parametru odpowiadającego liczbie boków matrycy, θ.

Intuicja

Jeśli po 40 rzutach zaobserwowałeś 10 czerwonych, 10 niebieskich, 10 zielonych i 10 żółtych, wydaje się, że θ powinno osiągać wartość szczytową przy 4, a stronniczość przewracania z każdej strony to rozkłady wyśrodkowane na 1/4.

θ ma trywialną dolną granicę, będącą liczbą różnych stron zaobserwowanych w danych.

Górna granica jest nadal nieznana. Może istnieć piąta strona, która prawdopodobnie miałaby niski poziom uprzedzeń. Im więcej danych obserwujesz bez piątej kategorii, tym wyższe prawdopodobieństwo z tyłu θ = 4.

Podejście

Użyłem JAGS do podobnych problemów (przez R i rjags), co wydaje się odpowiednie tutaj.

W odniesieniu do danych, powiedzmy, obs <- c(10, 10, 10, 10)że odpowiadają one obserwacjom w powyższym przykładzie.

Myślę, że obserwacje należy modelować z rozkładem wielomianowym obs ~ dmulti(p, n), gdzie p ~ ddirch(alpha)i n <- length(obs).

θ jest powiązany z liczbą kategorii implikowanych przez alpha, więc w jaki sposób mogę modelować w alphacelu uwzględnienia różnej możliwej liczby kategorii?

Alternatywy?

Jestem całkiem nowy w analizach bayesowskich, więc może całkowicie szczekam na niewłaściwym drzewie, czy istnieją alternatywne modele, które mogą dostarczyć różnych informacji na temat tego problemu?

Wielkie dzięki! David

Jest to interesujący problem zwany „próbkowaniem gatunków”, który przez lata cieszył się dużym zainteresowaniem i obejmuje wiele innych problemów z szacowaniem (takich jak przechwytywanie znaków). Wystarczy powiedzieć, że JAGS nie pomoże ci w tym przypadku - JAGS nie obsługuje łańcuchów Markowa o zmiennym wymiarze w iteracjach. Należy skorzystać ze schematu MCMC zaprojektowanego dla takich problemów, takich jak MCMC z odwracalnym skokiem.

Oto jedno podejście, które jest odpowiednie dla opisywanego modelu, które po raz pierwszy spotkałem w pracy Jeffa Millera ( opracowane ).

Część I (oryginalne pytanie)

Przyjmę jedno założenie, że obserwacja danej kategorii implikuje istnienie kategorii o niższej randze. Oznacza to, że obserwowanie rzutu kostką na stronie 9 implikuje istnienie stron 1-8. To nie muszą być w ten sposób - kategorie może być dowolna - ale będę zakładać, więc w moim przykładzie. Oznacza to, że można zaobserwować wartości 0, w przeciwieństwie do innych problemów z szacowaniem gatunków.

Powiedzmy, że mamy próbkę wielomianową,

$$Y = \{y_1, y_2, \dots, y_m, y_{m+1}, \dots, y_{n} \} \sim \mathcal{M}(\{p_1, p_2, \dots, p_m, p_{m+1}, \dots, p_n\})$$

gdzie to maksymalna zaobserwowana kategoria, to (nieznana) liczba kategorii, a wszystkie równe 0. Parametr jest skończony i potrzebujemy przeor za to. Każdy dyskretny, odpowiedni wcześniej ze wsparciem dla będzie działał; weźmy na przykład skróconą przez zero Poissona:$m$$n$$\{y_{m+1},\dots,y_{n}\}$$n$$[1, \infty)$

$$n \sim \mathcal{P}(\lambda), n > 0$$

Dirichlet jest wygodnym wstępem dla prawdopodobieństw wielomianowych,

$$P = \{ p_1, \dots, p_n \} \sim \mathcal{D}(\{ \alpha_1, \dots, \alpha_n \})$$

I dla uproszczenia załóżmy, że .$\alpha_1 = \alpha_2 = \dots = \alpha_n = \tilde{\alpha}$

Aby uprościć problem, marginalizujemy wagi:

$$p(Y|\tilde{\alpha}, n) = \int_P p(Y|P, n)p(P|\tilde{\alpha}, n) dP$$

Co w tym przypadku prowadzi do dobrze zbadanego rozkładu wielomianowego Dirichleta . Celem jest zatem oszacowanie warunkowego tylnego,

$$p(n|Y, \tilde{\alpha}, \lambda) = \frac{ p(Y|n, \tilde{\alpha}) p(n|\lambda) }{ p(Y|\tilde{\alpha}, \lambda) }$$

Tam, gdzie wyraźnie zakładam, że i są stałymi hiperparametrami. Łatwo zauważyć, że:$\tilde{\alpha}$$\lambda$

$$p(Y|\tilde{\alpha}, \lambda) = \sum_{n=1}^\infty p(Y|n, \tilde{\alpha}) p(n|\lambda)$$

Gdzie gdzie . Ta nieskończona seria powinna zbiegać się dość szybko (o ile ogon poprzedzający nie jest zbyt ciężki), a zatem jest łatwy do przybliżenia. Dla obciętego Poissona ma on postać:$p(Y|n, \tilde{\alpha}) = 0$$n < m$

$$p(Y|\tilde{\alpha}, \lambda) = \frac{1}{(e^\lambda - 1)} \sum_{n=m}^\infty \frac{\Gamma(n\tilde{\alpha})\prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})}{\Gamma(n\tilde{\alpha} + \sum_{i=1}^n y_i) \Gamma(\tilde{\alpha})^n} \cdot \frac{\lambda^n}{n!}$$

Prowadzący do:

$$p(n|Y,\tilde{\alpha}, \lambda) = \frac{\Gamma(n\tilde{\alpha})\prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})}{\Gamma(n\tilde{\alpha} + \sum_{i=1}^n y_i) \Gamma(\tilde{\alpha})^n} \cdot \frac{\lambda^n}{n!} \cdot \left(\sum_{j=m}^\infty \frac{\Gamma(j\tilde{\alpha})\prod_{i=1}^j \Gamma(y_i + \tilde{\alpha})}{\Gamma(j\tilde{\alpha} + \sum_{i=1}^j y_i) \Gamma(\tilde{\alpha})^j} \cdot \frac{\lambda^j}{j!}\right)^{-1}$$

Który ma wsparcie dla . W tym przypadku nie ma potrzeby MCMC, ponieważ nieskończoną serię w mianowniku reguły Bayesa można przybliżać bez większego wysiłku.$[m, \infty)$

Oto niechlujny przykład w języku R:

logPosteriorN <- function(max, Y, lambda, alpha){
    m <- length(Y)
    sumy <- sum(Y)
    pp <- sapply(1:max, function(j){
        prior <- log(lambda)*j - log(exp(lambda)-1) - lgamma(j+1)
        posterior <- lgamma(alpha*j) + sum(lgamma(Y + alpha)) - j*lgamma(alpha) - lgamma(sumy + j*alpha)
        if( j > m ) { posterior <- posterior + (j-m)*lgamma(alpha) } 
        else if( j < m ) { posterior = -Inf }
        prior + posterior
        })
    evidence <- log(sum(exp(pp))) # there's no check that this converges
    pp - evidence
}

## with even representation of sides
Y <- c(10, 10, 10, 10)
post <- logPosteriorN(30, Y, 10, 1.2)
plot(1:30, exp(post), pch=19, type="b")

## with uneven representation of sides
Y <- c(1, 2, 1, 0, 0, 2, 1, 0, 1)
post <- logPosteriorN(30, Y, 10, 1.2)
plot(1:30, exp(post), pch=19, type="b")

Twoja intuicja jest prawidłowa: rzadkie pobieranie próbek między kategoriami prowadzi do większej niepewności co do całkowitej liczby kategorii. Jeśli chcesz traktować jako nieznany parametr, musisz użyć MCMC i alternatywnych aktualizacji i .$\tilde{\alpha}$$n$$\tilde{\alpha}$

Oczywiście jest to jedno podejście do szacowania. Z łatwością znajdziesz inne (smaki bayesowskie i nie bayesowskie) z odrobiną poszukiwań.

Część II (odpowiedź na komentarz)

$Y = \{y_1, \dots, y_m, y_{m+1}, \dots, y_n \}$ to częściowo obserwowany wektor wielomianowy z odpowiednimi prawdopodobieństwami : $\Omega = \{\omega_1, \dots, \omega_m, \omega_{m+1}, \dots, \omega_n\}$

$$\mathrm{Pr}(Y|\Omega, n) = \frac{\Gamma(\sum_{i=1}^n y_i + 1)}{\prod_{i=1}^n \Gamma(y_i + 1) } \prod_{i=1}^n \omega_i^{y_i}$$

Gdzie , i ale poza tym indeksy są arytmetyczne. Tak jak poprzednio, problemem jest wnioskowanie o prawdziwej liczbie kategorii , a zaczynamy od znaku poprzedzającego na takiego jak Poisson o zerowym skróceniu: $y \in \mathbb{N}$$y_1 \dots y_m > 0$$y_{m+1} \dots y_n = 0$$n$$n$

$$\mathrm{Pr}(n|\lambda) = \frac{\lambda^{n}}{(\exp\{\lambda\} - 1)n!},~n \in \mathbb{Z}^+$$

Podobnie jak poprzednio, traktujemy prawdopodobieństwa wielomianowe jak Dirichleta dystrybuowanego z hiperparametrem symetrycznym , tj. Dla danego , $\Omega$$\tilde{\alpha}$$n$

$$\mathrm{Pr}(\Omega|\tilde{\alpha}, n) = \frac{\Gamma(n\tilde{\alpha})}{\Gamma(\tilde{\alpha})^n} \prod_{i=1}^n \omega_i^{\tilde{\alpha}-1}$$

Całkowanie (marginalizacja) ponad wektorem prawdopodobieństw daje wielomianowy Dirichlet:

$$\mathrm{Pr}(Y|\tilde{\alpha}, n) = \int \mathrm{Pr}(Y|\Omega, n) \mathrm{Pr}(\Omega|\tilde{\alpha}, n) = \frac{\Gamma(n \tilde{\alpha})} {\Gamma(\sum_{i=1}^n y_i + n \tilde{\alpha}) \Gamma(\tilde{\alpha})^n} \prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})$$

Tutaj odbiegamy od modelu w części I powyżej. W części I zawarto niejawne uporządkowanie kategorii: na przykład w matrycy stronnej kategorie (boki) mają domyślny porządek, a obserwacja dowolnej kategorii implikuje istnienie mniejszych kategorii . W części II mamy częściowo zaobserwowany wielomianowy losowy wektor, który nie ma niejawnego uporządkowania. Innymi słowy, dane reprezentują nieuporządkowany podział punktów danych na obserwowanych kategorii. Oznaczę nieuporządkowaną partycję wynikającą z powiększonej o nieobserwowane kategorie, jako .$n$$i \in \{1 \dots n\}$$j < i$$m \leq n$$Y$$n-m$$\mathcal{P}[Y]$

Prawdopodobieństwo nieuporządkowanej partycji uwarunkowanej prawdziwą liczbą kategorii można znaleźć, biorąc pod uwagę liczbę permutacji kategorii, które prowadzą do tej samej partycji: $n$

$$\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, n) = \frac{n!}{(n-m)!} \mathrm{Pr}(Y|\tilde{\alpha}, n)$$

I można to zintegrować za pomocą aby uzyskać: $n$

$$\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, \lambda) = \sum_{j=m}^{\infty} \mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, n) \mathrm{Pr}(n|\lambda)$$

Użycie reguły Bayesa do pobrania tylnej:

$$\mathrm{Pr}(n|\mathcal{P}[Y], \tilde{\alpha}, \lambda) = \frac{\mathrm{Pr}(\mathcal{P}[Y]|n, \tilde{\alpha}) \mathrm{Pr}(n|\lambda)}{\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, \lambda)}$$

Po prostu podłącz powyższe definicje. Ponownie mianownik jest nieskończoną serią, która szybko się zbiegnie: w tym prostym modelu MCMC nie wymaga odpowiedniego przybliżenia.

Zmieniając kod R z części I:

logPosteriorN_2 <- function(max, Y, lambda, alpha){
    m <- length(Y)
    sumy <- sum(Y)
    pp <- sapply(1:max, function(j){
        prior <- log(lambda)*j - log(exp(lambda)-1) - lgamma(j+1)
        likelihood <- lchoose(j, m) + lgamma(m + 1) + lgamma(alpha*j) + sum(lgamma(Y + alpha)) - j*lgamma(alpha) - lgamma(sumy + j*alpha)
        if( j > m ) { likelihood <- likelihood + (j-m)*lgamma(alpha) } 
        else if( j < m ) { likelihood = -Inf }
        prior + likelihood
        })
    evidence <- log(sum(exp(pp))) # there's no check that this converges
    pp - evidence
}

Y_1 <- rep(10, 15)
pos_1 <- logPosteriorN_2(50, Y_1, 6, 1)
plot(1:50, exp(pos_1))