Wie viele Seiten hat ein Würfel? Bayesianische Folgerung in JAGS

Problem

Ich möchte auf ein System schließen, das analog dazu ist, mit einer unbekannten Anzahl von Seiten zu sterben. Der Würfel wird mehrmals gewürfelt, wonach ich eine Wahrscheinlichkeitsverteilung über einen Parameter ableiten möchte, der der Anzahl der Seiten des Würfels entspricht, θ.

Intuition

Wenn Sie nach 40 Rollen 10 Rot-, 10 Blau-, 10 Grün- und 10 Gelbtöne beobachtet haben, scheint es, dass θ bei 4 seinen Höhepunkt erreichen sollte, und die Verzerrungen beim Rollen jeder Seite sind Verteilungen, die auf 1/4 zentriert sind.

θ hat eine triviale Untergrenze, dh die Anzahl der verschiedenen Seiten, die in den Daten beobachtet werden.

Die Obergrenze ist noch unbekannt. Es könnte eine fünfte Seite geben, die wahrscheinlich eine geringe Tendenz hätte. Je mehr Daten Sie beobachten, denen eine fünfte Kategorie fehlt, desto höher ist die hintere Wahrscheinlichkeit von θ = 4.

Ansatz

Ich habe JAGS für ähnliche Probleme (über R und rjags) verwendet, was hier angemessen erscheint.

Nehmen wir obs <- c(10, 10, 10, 10)in Bezug auf die Daten an, dass sie den Beobachtungen im obigen Beispiel entsprechen.

Ich denke, die Beobachtungen sollten mit einer multinomialen Verteilung modelliert werden obs ~ dmulti(p, n), wo p ~ ddirch(alpha)und n <- length(obs).

θ ist mit der Anzahl der Kategorien verknüpft, die durch impliziert werden alpha. Wie kann ich also modellieren alpha, um verschiedene mögliche Anzahlen von Kategorien zu erfassen?

Alternativen?

Ich bin ziemlich neu in Bayes'schen Analysen, könnte also den falschen Baum völlig bellen. Gibt es alternative Modelle, die unterschiedliche Einblicke in dieses Problem liefern könnten?

Danke vielmals! David

Dies ist ein interessantes Problem, das als "Artenprobenahme" bezeichnet wird und im Laufe der Jahre viel Aufmerksamkeit erhalten hat und viele andere Schätzprobleme umfasst (z. B. die Wiedererfassung von Markierungen). Es genügt zu sagen, dass JAGS Ihnen in diesem Fall nicht hilft - JAGS kann keine Markov-Ketten mit einer variablen Dimension über Iterationen hinweg verarbeiten. Man muss auf ein MCMC-Schema zurückgreifen, das für solche Probleme ausgelegt ist, wie beispielsweise das MCMC mit reversiblem Sprung.

Hier ist ein Ansatz, der für das spezifische Modell geeignet ist, das Sie beschreiben und das ich zum ersten Mal in der Arbeit von Jeff Miller ( arxived ) kennengelernt habe .

Teil I (Originalfrage)

Eine Annahme, die ich machen werde, ist, dass eine Beobachtung einer bestimmten Kategorie die Existenz von Kategorien mit einem niedrigeren Rang impliziert. Das heißt, das Beobachten eines Würfelwurfs auf Seite 9 impliziert die Existenz der Seiten 1-8. Es muss nicht so sein - die Kategorien könnten beliebig sein - aber ich gehe in meinem Beispiel davon aus. Dies bedeutet, dass im Gegensatz zu anderen Problemen bei der Artenschätzung 0-Werte beobachtbar sind.

Nehmen wir an, wir haben eine multinomiale Stichprobe.

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

Dabei ist die maximal beobachtete Kategorie, die (unbekannte) Anzahl von Kategorien und alle gleich 0. Der Parameter ist endlich, und wir brauchen ein Prior dafür. Jeder diskrete, ordnungsgemäße Vorgänger mit Unterstützung für wird funktionieren; Nehmen Sie zum Beispiel ein Poisson mit Null-Verkürzung:$m$$n$$\{y_{m+1},\dots,y_{n}\}$$n$$[1, \infty)$

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

Ein geeigneter Prior für die multinomialen Wahrscheinlichkeiten ist das Dirichlet,

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

Und nehmen Sie einfach an: .$\alpha_1 = \alpha_2 = \dots = \alpha_n = \tilde{\alpha}$

Um das Problem leichter handhabbar zu machen, werden die Gewichte marginalisiert:

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

Was in diesem Fall zur gut untersuchten Dirichlet-Multinomialverteilung führt . Das Ziel ist dann, den bedingten posterioren zu schätzen,

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

Wo ich explizit davon dass und feste Hyperparameter sind. Es ist leicht zu sehen, dass:$\tilde{\alpha}$$\lambda$

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

Wobei wobei . Diese unendliche Reihe sollte ziemlich schnell konvergieren (solange der Schwanz des Prior nicht zu schwer ist) und ist daher leicht zu approximieren. Für den abgeschnittenen Poisson hat es die Form:$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!}$$

Führend zu:

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

Welches hat Unterstützung auf . In diesem Fall ist MCMC nicht erforderlich, da die unendliche Reihe im Nenner der Bayes-Regel ohne großen Aufwand angenähert werden kann.$[m, \infty)$

Hier ist ein schlampiges Beispiel in 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")

Ihre Intuition ist richtig: Eine spärliche Auswahl über Kategorien hinweg führt zu einer größeren Unsicherheit über die Gesamtzahl der Kategorien. Wenn Sie als unbekannten Parameter behandeln möchten, müssen Sie MCMC und alternative Aktualisierungen von und .$\tilde{\alpha}$$n$$\tilde{\alpha}$

Dies ist natürlich ein Ansatz zur Schätzung. Mit ein wenig Suche finden Sie leicht andere (mit Bayes'schen und nicht-Bayes'schen Aromen).

Teil II (Antwort auf Kommentar)

$Y = \{y_1, \dots, y_m, y_{m+1}, \dots, y_n \}$ ist ein teilweise beobachteter multinomialer Vektor mit entsprechenden Wahrscheinlichkeiten : $\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}$$

Wobei , und aber ansonsten sind die Indizes willkürlich. Nach wie vor besteht das Problem darin, die wahre Anzahl der Kategorien , und wir beginnen mit einem Prior auf wie beispielsweise einem Poisson mit Null-Verkürzung: $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}^+$$

Ebenso wie zuvor behandeln wir die multinomialen Wahrscheinlichkeiten als Dirichlet, verteilt mit einem symmetrischen Hyperparameter , dh für ein gegebenes , $\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}$$

Das Integrieren (Marginalisieren) über den Wahrscheinlichkeitsvektor ergibt das multinomiale 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})$$

Hier weichen wir vom Modell in Teil I oben ab. In Teil I gab es eine implizite Reihenfolge nach Kategorien: Beispielsweise haben in einem seitigen Würfel die Kategorien (Seiten) eine implizite Reihenfolge, und die Beobachtung einer Kategorie impliziert die Existenz kleinerer Kategorien . In Teil II haben wir einen teilweise beobachteten multinomialen Zufallsvektor, der keine implizite Ordnung hat. Mit anderen Worten, die Daten stellen eine ungeordnete Aufteilung der Datenpunkte in beobachtete Kategorien dar. Ich werde die ungeordnete Partition, die sich aus ergibt, ergänzt durch nicht beobachtete Kategorien, als .$n$$i \in \{1 \dots n\}$$j < i$$m \leq n$$Y$$n-m$$\mathcal{P}[Y]$

Die Wahrscheinlichkeit der ungeordneten Partition, die von einer wahren Anzahl von Kategorien abhängig ist , kann unter Berücksichtigung der Anzahl von Permutationen von Kategorien ermittelt werden, die zu derselben Partition führen: $n$

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

Und dies kann über integriert werden, um zu ergeben: $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)$$

Verwenden der Bayes-Regel zum Abrufen des Seitenzahns:

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

Schließen Sie einfach die obigen Definitionen an. Auch hier ist der Nenner eine unendliche Reihe, die schnell konvergiert: In diesem einfachen Modell muss MCMC keine angemessene Annäherung geben.

Durch Ändern des R-Codes aus Teil 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))