Resolver analíticamente muestreo con o sin reemplazo después binomial Poisson/negativa

Question

Versión corta

Estoy tratando de resolver analíticamente/aproximada del compuesto probabilidad de que los resultados de independiente de Poisson empates y más muestreo con o sin reemplazo (no me importa que uno). Quiero usar la probabilidad con MCMC (Stan), así que tengo la solución sólo hasta un término constante. Finalmente, quiero un modelo en el cual el proceso inicial de los sorteos son de neg. distribución binomial, pero creo que voy a ser capaz de llegar con una solución para el caso de Poisson.

Es bien posible que la solución no es factible (no entiendo las matemáticas suficiente para ser capaz de decir si esto es un simple o un problema muy difícil). Yo soy así también interesado en las aproximaciones, los resultados negativos o de la intuición de por qué el problema es, probablemente, insolubles (por ejemplo, comparar a un conocido problema difícil). Enlaces a documentos útiles/teoremas/trucos que me ayudará a avanzar son buenas respuestas, incluso si su conexión con el problema en cuestión no está totalmente resuelto.

Declaración Formal

Más formalmente, la primera $Y = (y_1, ..., y_N), y_n \sim Pois(\lambda_n)$ es elegido de forma independiente y luego me muestra $k$ artículos al azar de entre todos los de $Y$ conseguir $Z = (z_1,...,z_N)$. I. e. Puedo dibujar $k$ bolas de colores a partir de una urna donde la cantidad de bolas de color $n$ se extrae de $Pois(\lambda_n)$. Aquí, $k$ se supone conocida y fija y nos condición en $\sum_n y_n \geq k$. Técnicamente, el muestreo se hace sin reemplazo, pero asumiendo muestreo con reemplazo no debe ser gran cosa.

He probado con dos enfoques para resolver para muestreo sin reemplazo (como este parecía ser el más fácil de caso debido a que algunos de los términos de cancelación fuera), pero se quedó atascado con ambos. La probabilidad cuando el muestreo sin reemplazo es:

$$P(Z = (z_1, ..., z_N) | \Lambda = (\lambda_1, ..., \lambda_N)) = \frac{ \sum_{Y;\forall n: y_n \geq z_n} \left( \frac{\prod_{n=1}^N{y_n \elegir z_n}}{ {\sum_{n=1}^N y_n} \elegir k} \prod_{n=1}^N de Poisson(y_n |\lambda_n) \right) }{ P(\sum_n y_n \geq k|\Lambda) }$$

EDIT: El "intento de la sección de soluciones fue removido como la solución en la respuesta no se puede construir sobre ellos (y es mejor)

Answer 1

El caso sin reemplazo

Si usted tiene $n$ independiente de Poisson variables de distribución

$$Y_i \sim \text{Pois$(\lambda_i)$}$$

and condition on

$$\sum_{j=i}^n Y_j = K$$

then

$$\lbrace Y_i \rbrace \sim \text{Multinom} \left(K,\left(\frac{\lambda_i}{\sum_{j=1}^n \lambda_j} \right)\right)$$

---

So you could fill your urn with the $n$ times $Y_i$ colored balls like first drawing the value for the total $K$ (which is Poisson distributed cutoff by the condition $K \geq k$) and then fill the urn with $K$ balls according to the multinomial distribution.

This filling of the urn with $K$ balls, according to a multinomial distribution, is equivalent to drawing for each ball independently the color from the categorical distribution. Then you can consider the first $k$ balls that have been added to the urn as defining the random sample $\lbrace Z_i \rbrace$ (when this sample is drawn without replacement) and the distribution for this is just another multinomial distributed vector:

$$\lbrace Z_i \rbrace \sim \text{Multinom} \left(k,\left(\frac{\lambda_i}{\sum_{j=1}^n \lambda_j} \right)\right)$$

simulación

##### simulating sample process for 3 variables #######


# settings
set.seed(1)
k = 10
lambda = c(4, 4, 4)
trials = 1000000

# observed counts table
Ocounts = array(0, dim = c(k+1, k+1, k+1))

for (i in c(1:trials)) {
  # draw poisson with limit sum(n) >= k
  repeat {
    Y = rpois(3,lambda)
    if (sum(Y) >= k) {break}
  }
  # setup urn
  urn <- c(rep(1, Y[1]), rep(2, Y[2]), rep(3, Y[3]))
  # draw from urn
  s <- sample(urn, k, replace=0)
  Z = c(sum(s==1),sum(s==2),sum(s==3))
  Ocounts[Z[1]+1, Z[2]+1, Z[3]+1] = Ocounts[Z[1]+1, Z[2]+1, Z[3]+1] + 1
}



# comparison
observed = rep(0, 0.5*k*(k+1))
expected = rep(0, 0.5*k*(k+1))   
t = 1

for (z1 in c(0:k)) {
  for (z2 in c(0:(k-z1))) {  
    z3 = k-z1-z2 
    observed[t] = Ocounts[z1+1, z2+1, z3+1]
    expected[t] = trials*dmultinom(c(z1, z2, z3), prob=lambda)
    t = t+1
  }
}

plot(observed,expected)
x2 <- sum((observed-expected)^2/expected)
pvalue <- 1-pchisq(x2, 66-1)

resultados

> # results from chi-sq test
> x2
[1] 75.49286
> pvalue
[1] 0.1754805

Binomial negativa

Los argumentos sería el mismo para el caso de una distribución binomial negativa que los resultados (bajo ciertas condiciones) en una de Dirichlet-distribución multinomial.

A continuación es un ejemplo de simulación

##### simulating sample process for 3 variables #######

# dirichlet multinomial for vectors of size 3
ddirmultinom =  function(x1,x2,x3,p1,p2,p3) {
  (factorial(x1+x2+x3)*gamma(p1+p2+p3)/gamma(x1+x2+x3+p1+p2+p3))/
  (factorial(x1)*gamma(p1)/gamma(x1+p1))/
  (factorial(x2)*gamma(p2)/gamma(x2+p2))/
  (factorial(x3)*gamma(p3)/gamma(x3+p3))
}

# settings
set.seed(1)
k = 10
theta = 1
lambda = c(4,4,4)
trials = 1000000

# calculating negative binomials pars
means = lambda
vars = lambda*(1+theta)

ps = (vars-means)/(vars)
rs = means^2/(vars-means)


# observed counts table
Ocounts = array(0, dim = c(k+1,k+1,k+1))

for (i in c(1:trials)) {
  # draw poisson with limit sum(n) >= k
  repeat {
    Y = rnbinom(3,rs,ps)
    if (sum(Y) >= k) {break}
  }
  # setup urn
  urn <- c(rep(1,Y[1]),rep(2,Y[2]),rep(3,Y[3]))
  # draw from urn
  s <- sample(urn,k,replace=0)
  Z = c(sum(s==1),sum(s==2),sum(s==3))
  Ocounts[Z[1]+1,Z[2]+1,Z[3]+1] = Ocounts[Z[1]+1,Z[2]+1,Z[3]+1] + 1
}



# comparison
observed = rep(0,0.5*k*(k+1))
expected = rep(0,0.5*k*(k+1))   
t = 1

for (z1 in c(0:k)) {
  for (z2 in c(0:(k-z1))) {  
    z3 = k-z1-z2 
    observed[t]=Ocounts[z1+1,z2+1,z3+1]
    expected[t]=trials*ddirmultinom(z1,z2,z3,lambda[1]/theta,lambda[2]/theta,lambda[3]/theta)
    t = t+1
  }
}

plot(observed,expected)
x2 <- sum((observed-expected)^2/expected)
pvalue <- 1-pchisq(x2,66-1)

# results from chi-sq test
x2
pvalue