Сколько сторон у кубика? Байесовский вывод в JAGS

проблема

Я хотел бы сделать некоторые выводы о системе, аналогичной смерти с неизвестным числом сторон. Матрица бросается несколько раз, после чего я хотел бы вывести распределение вероятностей по параметру, соответствующему количеству сторон, которые имеет матрица, θ.

Интуиция

Если после 40 бросков вы наблюдали 10 красных, 10 синих, 10 зеленых и 10 желтых, кажется, что θ должен достигать пика в 4, а уклоны прокатки с каждой стороны - это распределения с центром в 1/4.

θ имеет тривиальную нижнюю границу, являющуюся количеством различных сторон, наблюдаемых в данных.

Верхняя граница пока неизвестна. Там может быть пятая сторона, которая, вероятно, будет иметь низкий уклон. Чем больше данных вы наблюдаете, не имея пятой категории, тем выше апостериорная вероятность θ = 4.

Подходить

Я использовал JAGS для подобных проблем (через R и RJAGS), которые здесь уместны.

Что касается данных, допустим, obs <- c(10, 10, 10, 10)чтобы соответствовать наблюдениям в примере выше.

Я думаю, что наблюдения должны быть смоделированы с многочленным распределением obs ~ dmulti(p, n), где p ~ ddirch(alpha)и n <- length(obs).

θ связано с количеством категорий, подразумеваемых alpha, так как я могу моделировать, alphaчтобы охватить различное возможное количество категорий?

Альтернативы?

Я довольно плохо знаком с байесовским анализом, так что, может быть, я полностью ошибаюсь, существуют ли альтернативные модели, которые могут по-разному понять эту проблему?

Большое спасибо! Дэвид

Эта интересная проблема называлась «выборка видов», которой много лет уделялось много внимания, и она включает в себя множество других проблем оценки (таких как повторная поимка). Достаточно сказать, что JAGS не поможет вам в этом случае - JAGS не может обрабатывать цепочки Маркова с переменным измерением на протяжении итераций. Нужно прибегнуть к схеме MCMC, разработанной для таких задач, как, например, обратимый скачок MCMC.

Вот один подход, который подходит для конкретной модели, которую вы описываете, с которой я впервые столкнулся в работе Джеффа Миллера ( arxived ).

Часть I (Оригинальный вопрос)

Одно из предположений, которое я сделаю, состоит в том, что наблюдение данной категории подразумевает существование категорий меньшего ранга. То есть, наблюдение броска кубика на стороне 9 подразумевает существование сторон 1-8. Так не должно быть - категории могут быть произвольными - но я предполагаю, что так в моем примере. Это означает, что наблюдаются 0 значений, в отличие от других проблем оценки видов.

Допустим, у нас есть полиномиальный образец,

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

где - максимальная наблюдаемая категория, - (неизвестное) количество категорий, и все равны 0. Параметр конечен, и нам нужно априор для этого. Любой дискретный, правильный предшествующий с поддержкой будет работать; Возьмем для примера усеченный до нуля Пуассон:n { y m + 1 , … , y n } n [ 1 , ∞ $m$$n$$\{y_{m+1},\dots,y_{n}\}$$n$$[1, \infty)$

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

Удобным приоритетом для полиномиальных вероятностей является Дирихле,

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

А для упрощенно предположим, что .$\alpha_1 = \alpha_2 = \dots = \alpha_n = \tilde{\alpha}$

Чтобы сделать задачу более удобной, мы выделим веса:

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

Что в данном случае приводит к хорошо изученному дирихле-полиномиальному распределению. Цель состоит в том, чтобы затем оценить условный апостериор,

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

Где я явно предполагаю, что и являются фиксированными гиперпараметрами. Легко увидеть, что:$\tilde{\alpha}$$\lambda$

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

Где где . Этот бесконечный ряд должен сходиться довольно быстро (пока хвост предыдущего не слишком тяжелый), и поэтому его легко аппроксимировать. Для усеченного Пуассона он имеет вид:п < $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!}$$

Ведущий к:

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

Который имеет поддержку на . В этом случае нет необходимости в MCMC, поскольку бесконечные ряды в знаменателе правила Байеса могут быть аппроксимированы без особых усилий.$[m, \infty)$

Вот небрежный пример в 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")

Ваша интуиция верна: редкая выборка по категориям приводит к большей неопределенности относительно общего количества категорий. Если вы хотите использовать как неизвестный параметр, вам нужно будет использовать MCMC и альтернативные обновления и . n ˜ $\tilde{\alpha}$$n$$\tilde{\alpha}$

Конечно, это один из подходов к оценке. Вы легко найдете другие (с байесовским и не байесовским вкусами) с небольшим поиском.

Часть II (Ответ на комментарий)

$Y = \{y_1, \dots, y_m, y_{m+1}, \dots, y_n \}$ - частично наблюдаемый полиномиальный вектор с соответствующими вероятностями : $\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}$$

Где , и но в остальном индексы являются случайными. Как и прежде, проблема заключается в том, чтобы вывести истинное количество категорий , и мы начнем с априора такого как усеченный до нуля Пуассон: $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}^+$$

Также как и прежде, мы рассматриваем полиномиальные вероятности как Dirichlet, распределенный с симметричным гиперпараметром , т.е. для данного , $\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}$$

Интегрирование (маргинализация) по вектору вероятностей дает многочлен Дирихле:

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

Вот где мы расходимся с моделью в части I выше. В первой части было неявное упорядочение по категориям: например, в стороннем кристалле категории (стороны) имеют неявное упорядочение, и наблюдение любой категории подразумевает существование меньших категорий . Во второй части мы имеем частично наблюдаемый полиномиальный случайный вектор, который не имеет неявного упорядочения. Другими словами, данные представляют собой неупорядоченное разбиение точек данных на наблюдаемых категорий. Я буду обозначать неупорядоченное разбиение, полученное в результате увеличенное на ненаблюдаемые категории , как .$n$$i \in \{1 \dots n\}$$j < i$$m \leq n$$Y$$n-m$$\mathcal{P}[Y]$

Вероятность неупорядоченного разбиения, обусловленного истинным числом категорий , можно определить, учитывая количество перестановок категорий, которые приводят к одному и тому же разбиению: $n$

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

И это можно интегрировать по чтобы получить: $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)$$

Использование правила Байеса для получения апостериорных:

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

Просто подключитесь из приведенных выше определений. Опять же, знаменатель представляет собой бесконечный ряд, который будет быстро сходиться: в этой простой модели MCMC не нужно давать адекватное приближение.

Изменяя код R из части 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))