Integración numérica / aproximada sobre variables gamma independientes

Question

Declaración Del Problema

Por un problema en la biología, la estoy probando una distribución conjunta de la forma:

$$ X \sim Multinomial(\frac{\theta_1}{\sum \theta_i}, ...,\frac{\theta_n}{\sum{\theta_i}}) \\ \theta_i \sim Gamma(\alpha_i, \beta_i) $$ (donde yo uso la forma de tasa de configuración de parámetros de la distribución Gamma).

Estoy interesado en el (logaritmo del) distribución marginal $$ P(X | \alpha_{1..n}, \beta_{1..n}) = \int Multinomial(X|\Theta) \prod_i Gamma(\theta_i | \alpha_i, \beta_i) d \Theta $$

Sé que, en el caso especial de $\beta_1 = \beta_2 = ... = \bar\beta$, la distribución marginal es de Dirichlet Multinomial (DM). La DM de distribución, sin embargo no se ajusta a los datos que estoy encontrando, los datos están más dispersos que lo que el DM de distribución puede acomodar. ¿Cómo puedo evaluar/aproximado de la general integral, al menos, un poco de manera eficiente? Yo necesitaría para esto, dentro de una gran inferencia, por lo que para muchas combinaciones de los parámetros.

El $\alpha_i$ e $\beta_i$ provienen de un nivel más profundo de la modelo y son limitados por lo que $\sum E(log(\theta_i)) = 0$.

Es posible que esta realidad es un problema muy difícil, y el puntero de por qué es difícil que iba a ser una respuesta suficiente.

Cosas que he probado

Dejando $\beta_{1..n}$ a variar y el tratamiento de la $\theta_{1..n}$ como explícita de variables latentes, este modelo se ajusta a los datos razonablemente bien, pero es muy lento, incluso para las pequeñas $n$ debido a la gran cantidad de variables latentes (estoy usando Stan).

Yo no creo que existe una solución analítica para $P(X | \alpha_{1..n}, \beta_{1..n})$ en el caso general, pero sería posible escribir una aproximación o tener un esquema de integración numérica que me dejaría para el cálculo de la (log) de esta densidad de una manera eficiente? Se siente como la estructura de independiente Gammas podría ser de alguna manera explotable (de hecho, la integral puede ser modificado para ser más yo.yo.d Gammas de tomar ventaja de la Gamma de propiedades de escala).

Yo he probado un ingenuo Monte-carlo (esquema de la muestra de las gammas, calcular el multinomial densidad promedio de más muestras) que estaba perdidamente lento para convergen incluso con unas dimensiones y un sencillo importancia esquema de muestreo, a partir de muestras de la distribución Dirichlet como las propuestas de $\theta_1 / \sum \theta_i$ (a de la muestra cerca de la máxima verosimilitud área de la multinomial), pero eso sólo empeoró la situación.

He tratado de encontrar la aproximación de Laplace a $P(log(\Theta)|X, \alpha_i, \beta_i)$ que se ve razonablemente de Gauss-como para muchas combinaciones de parámetros, pero resulta ser problemático, como cuando $x_i = 0$ e $\alpha_i \leq 1$, el modo no está definido. $P(\Theta)|X, \alpha_i, \beta_i)$ no parecerse a los de Gauss.

A partir de la aproximación del punto de vista de la pregunta ¿Cuál es el valor esperado de modificaciones de la distribución Dirichlet? (problema de integración) está relacionado. Pero todos mis intentos de composición como de muy baja aproximaciones.

---

Ligeramente philosphical nota: creo que DM no se adapta a mis datos porque DM surge también de muestreo multinomial de la binomial negativa variables con un factor de Fano (media/varianza de la proporción de) - véase, por ejemplo, Analíticamente la solución de muestreo con o sin reemplazo después de Poisson/binomial Negativa. Neg. binomio es el Gamma-Poisson, donde el $\beta$ parámetro determina el factor de Fano. Pero mis datos corresponden a multinomial de muestreo de neg. binom variables en las que el factor de Fano varía

Answer 1

Aquí es un boceto de una solución. La estrategia general es introducir el escalar suma $\theta_T=\sum\theta_i, $ , a continuación, marginar sobre él, a continuación, utilizarlo para activar el complicado multidimenional suma de las $\theta_{1..n}$ a un problema de programación dinámica agresiva de los puntos de corte para podar el espacio de solución como la dinámica del programa se ejecuta.

Sabemos que $\theta_T$ es una suma de gamma RVs con la heterogeneidad de los parámetros. Desde $\beta_i>0$ para todos los $i$, las desviaciones son finitos, y si $n$ es lo suficientemente grande, este va a sucumbir a la CLT, por lo $\theta_T \rightarrow_d ~~ N(\mu_T,\sigma^2_T)$ donde $\mu_T =\sum \alpha_i/\beta_i$ e $\sigma^2_T=\sum \alpha_i/\beta_i^2.$ Esto nos da una creíble subconjunto de $\theta_T$ más que centrar nuestra atención en $\theta \in [\theta_{T_0},\theta_{T_1}]$. Incluso si $n$ no es muy grande, somos capaces de aproximar $p(\theta_T)$ a través de métodos numéricos y otros métodos, que se basan en los datos.

Para cada una de las $\theta$ seleccionados a partir de este intervalo, se elige un número entero $M$ , de modo que $\theta/M$ es la cantidad más pequeña de $\theta$ a ser asignados a $\theta_i$ en el reparto de $\theta$ entre el $\theta_{1..n}$. Esto convierte el problema en un simple problema de programación dinámica en la que $\sum \theta_{1..n}=\theta_T$ es equivalente a $\sum m_i=M$ donde $0\le m_i \le M, m_i \in \mathbb{N}. $ El número de soluciones es $\binom{M+n-1}{n-1}$. Para los números como $M=20, n=10, \binom{29}{9} \sim 10^8$. Afortunadamente, podemos implementar la programación dinámica con un ojo a la selección de las soluciones para que $P(\theta_i|\alpha_i,\beta_i)>\epsilon$ que se reducirá el número enormemente. Esto es equivalente a $\sum m_i=M$ requiere que,por ejemplo, $1 \le m_i \le M$, y ahora el número de soluciones es $\binom{(M-n)+n-1}{n-1}$. Para $M=20, n=10, \binom{19}{9}\sim 10^5,$ o 3 órdenes de magnitud más pequeña con una simple, obviamente-correcto requisito.

La inteligencia entra al memo-ize la dinámica del programa para la reutilización de las soluciones, a continuación, utilizar una regla para la variación de $\alpha_{1..n},\beta_{1..n}$ quizás para mantener constantes los valores de $\mu_T$ e $\sigma^2_T$, y que también maximizar la reutilización de la programación dinámica solución. Esta constancia de $\mu_T$ e $\sigma^2_T$ puede ser perfectamente una dentro de la esfera de sus necesidades.

Otros trucos para hacer que el problema manejable es lo que yo pondría en una más-que-solo-a-sketch respuesta a esta pregunta.

En resumen, $$ \begin{align} P(X|\alpha_{1..n},\beta_{1..n}) &= \sum_{\theta_T} P(X|\theta_T,\alpha_{1..n},\beta_{1..n})P(\theta_T| \alpha_{1..n},\beta_{1..n} )\\ &= \sum_{\theta_T} \sum_{\{\theta_{1..n}|\sum \theta_i=\theta_T\}} P(X|\alpha_{1..n},\beta_{1..n})N(\theta_T|\mu_T,\sigma^2_T). \end{align} $$

A continuación, $P(X|\alpha_{1..n},\beta_{1..n})= \binom{X}{x_1,x_2..x_n}\theta_1^{x_1} \cdots \theta_n^{x_1}$ sólo necesitan ser evaluados sólo para los $\theta_{1..n}$ valores que son de alta probabilidad en el primer lugar. Tenga en cuenta que el $\theta_{1..n}$ son implícitamente funciones de $\alpha_{1..n},\beta_{1..n}$.

Algunas observaciones

1. Encontrar un problema de programación dinámica en la naturaleza es siempre un montón de diversión. Sería bastante alto perfil de resultado de investigación, que debo pensar.

2. Hay números de $M$ e $n$ aquí donde este enfoque se vuelve computacionalmente muy difícil, incluso con programación dinámica. Si $n$ tiene que ser grande, se puede compensar mediante el uso de pequeñas $M$ o más $\epsilon.$ en última instancia, el problema sólo puede ser solucionable, si usted está dispuesto a a) el uso de algunos restrictiva de los priores, y/o b) aceptar una baja resolución de la solución.

3. La elección para marginar a los más de $\theta_T$ significa que el uso de la log-verosimilitud tendrá que ser término a término en la suma, en lugar de todas a la vez, que no es tan incómodo en el gran esquema de las cosas.

Es esto a lo que está buscando?

---

Anexo

Un 20 línea dinámico de programación de secuencia de comandos de python es capaz de enumerar las $\binom{M+n-1}{n-1} M=100, ~n=20$ con una solución de $4.910371215196107e+21$ en $0.18s$. Esto implica contar soluciones de visita en la programación dinámica de ordenación de la manera en que redundante soluciones no necesita ser revisada, sólo representaron. En esta enumeración, $\epsilon=0$. Visita cada estado enumeración es mucho más lento, por supuesto. Para M=20 y n=10, la enumeración completa que visita cada uno de los $10^8$ soluciones actualmente se lleva a 42s, que es donde se $\epsilon$ es valiosa, y los requisitos que permiten la resolución baja de presentación de informes de resultados.

Answer 2

Su pregunta no especifica los ensayos parámetro de la distribución multinomial, así que voy a usar el estándar genérico caso de que usted tiene $n$ juicio y $m$ categorías. Dejaré $\tilde{\boldsymbol{\theta}} = (\tilde{\theta}_1,...,\tilde{\theta}_m)$ denotar la normalizado de los valores de $\tilde{\theta}_k = \theta_k / \sum_i \theta_i$, por lo que el modelo es de interés:

$$\mathbf{X} | \boldsymbol{\theta} \sim \text{Mu}(n,\tilde{\boldsymbol{\theta}}) \quad \quad \quad \theta_k \sim \text{Ga}(\alpha_k, \beta_k).$$

---

The required integral: The marginal mass function for $\mathbf{X}$ is given by:

$$\begin{equation} \begin{aligned} p(\mathbf{x}|\boldsymbol{\alpha}, \boldsymbol{\beta}) &= \int \limits_\Theta \text{Mu}(\mathbf{x} | \tilde{\boldsymbol{\theta}}) \prod_{i=1}^m \text{Ga}(\theta_i | \alpha_i, \beta_i) \ d \boldsymbol{\theta} \\[6pt] &\overset{\mathbf{x}}{\propto} {n \choose \mathbf{x}} \int \limits_\Theta \frac{\prod_{i=1}^m \theta^{x_i + \alpha_i - 1} \exp( - \beta_i \theta_i )}{(\sum_{i=1}^m \theta_i)^n} \ d \boldsymbol{\theta}. \\[6pt] \end{aligned} \end{equation}$$

If $\beta_1 = ... = \beta_m = \beta$ then $\tilde{\boldsymbol{\theta}}$ has a Dirichlet distribution which means that $\mathbf{X}$ then has a multinomial-Dirichlet distribution. (In your question you incorrectly state that $\mathbf{X}$ has a Dirichlet distribution in this case. That can't be right because the Dirichlet distribution generates non-integer vectors of values.) There are various ways to approximate the above integral, and I will show you one of these methods below. (If time permits, I may come back to this answer and add other methods.)

---

Approximation by direct sampling (useful if hyperparameters are fixed): Since the object of interest is a mass function obtained via a known conditional distribution, the simplest way to approximate the integral is simply to generate a large number of conditioning values and approximate the integral by an analogous sum. To do this, choose some large number of simulations $H \in \mathbb{N}$ and generate independent values:

$$\theta_k^{(h)} \sim \text{Ga}(\alpha_k, \beta_k) \quad \quad \quad \text{para todo } k=1,...,m \text{ y } h =1,...,H.$$

You can then use the approximation:

$$\hat{p}_\text{DS}(\mathbf{x}|\boldsymbol{\alpha}, \boldsymbol{\beta}) = \frac{1}{H} \sum_{h=1}^H {n \elegir \mathbf{x}} \frac{\prod_{k=1}^m \theta_k^{(h) \ x_k}}{(\sum_{k=1}^m \theta_k^{(h)})^n} .$$

This method requires you to generate $H \times m$ scalar parameter values and then evaluate the mass function as an average of $H$ multinomial mass functions. Note that this method evaluates the mass function for a fixed set of parameter values $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$, so it is not useful if you would like your approximation to serve as a function of those parameter values.

---

Approximation by importance sampling (useful if hyperparameters are variables): The method of direct sampling builds the parameters $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ into the simulation, so it is useful when these are fixed (i.e., when you don't need the function to vary over these values). If you would like your approximating function to also be able to vary over these parameters then you can use importance sampling instead. To do this, choose some mid-ranged hyperparameter values $\bar{\alpha}$ and $\bar{\beta}$, choose some large number of simulations $H \in \mathbb{N}$ and generate independent values:

$$\theta_k^{(h)} \sim \text{Ga}(\bar{\alpha}, \bar{\beta}) \quad \quad \quad \text{para todo } k=1,...,m \text{ y } h =1,...,H.$$

You can then use the approximation:

$$\hat{p}_\text{ES}(\mathbf{x}|\boldsymbol{\alpha}, \boldsymbol{\beta}) = \frac{1}{H} \sum_{h=1}^H {n \elegir \mathbf{x}} \frac{\prod_{k=1}^m \theta_k^{(h) \ x_k}}{(\sum_{k=1}^m \theta_k^{(h)})^n} \cdot \prod_{k=1}^m \frac{\text{Ga}(\theta_k^{(h)} | \alpha_k, \beta_k)}{\text{Ga}(\theta_k^{(h)} | \bar{\alpha}, \bar{\beta})}.$$

This method requires you to generate $H \times m$ scalar parameter values and then evaluate the mass function as an average of $H$ ponderado multinomial masa de funciones, donde la ponderación se utiliza para ajustar la hyperparameters. Tenga en cuenta que este método trata el hyperparameters como variables que pueden ser ajustados en el cálculo de la aproximación de la función (sin generar nuevos valores simulados). Aquí hay algunos R código para implementar este último método:

#Load required libraries and set seed
library(stats);
library(matrixStats);
set.seed(1);

#Set parameter values for simulation
m <- 3;
A <- 2;
B <- 3;
H <- 10^4;

#Simulate values of theta
THETA <- array(rgamma(H*m, shape = A, rate = B), dim = c(H, m));

#Create function to generate approximation using importance sampling 
LOGMASS  <- function(x, alpha, beta) { 

    LOGWEIGHT <- array(0, dim = c(H, m));
    for (h in 1:H) {
    for (k in 1:m) {
        T1 <- dgamma(THETA[h,k], shape = alpha[k], rate = beta[k], log = TRUE);
        T2 <- dgamma(THETA[h,k], shape = A,        rate = B,       log = TRUE);
        LOGWEIGHT[h,k] <- (T1 - T2); } }
    LOGW <- rowSums(LOGWEIGHT);

    LOGTERMS <- rep(0, H);
    for (h in 1:H) {
    LOGTERMS[h] <- (stats::dmultinom(x, size = sum(x), prob = THETA[h, ],
                                     log = TRUE) + LOGW[h]); }

    matrixStats::logSumExp(LOGTERMS) - log(H); } 

#Find approximate mass function at an example point
x     <- c(3,3,4);
alpha <- c(4,3,2);
beta  <- c(2,1,2);
LOGMASS(x, alpha, beta);

[1] -4.132541

Este código le permitirá generar el registro de masa para su distribución en cualquier valor de entrada para la observación del vector y el hyperparameters. Estableciendo H razonablemente grande, usted deberá obtener una buena aproximación a la verdad de registro-función de masa. Algunos teoría estadística y la precisión de los límites para este tipo de simulación se puede encontrar en O'Neill (2009).

---