Analytically solving sampling with or without replacement after Poisson/Negative binomial

Question

Short version

I am trying to analytically solve/approximate the composite likelihood that results from independent Poisson draws and further sampling with or without replacement (I don't really care which one). I want to use the likelihood with MCMC (Stan), so I need the solution only up to a constant term. Ultimately, I want to model a process where the initial draws are from neg. binomial distribution, but I think I will be able to get there with a solution for the Poisson case.

It is well possible that the solution is not feasible (I don't understand the math enough to be able to tell whether this is a simple or very difficult problem). I am thus also interested in approximations, negative results or intuition why the problem is probably intractable (e.g. comparing to a known hard problem). Links to useful papers/theorems/tricks that will help me move forward are good answers even if their connection to the problem at hand is not fully worked out.

Formal statement

More formally, first $Y = (y_1, ..., y_N), y_n \sim Pois(\lambda_n)$ is drawn independently and then I sample $k$ items at random from all of $Y$ to get $Z = (z_1,...,z_N)$. I.e. I draw $k$ coloured balls from an urn where the amount of balls of color $n$ is drawn from $Pois(\lambda_n)$. Here, $k$ is assumed known and fixed and we condition on $\sum_n y_n \geq k$. Technically the sampling is done without replacement, but assuming sampling with replacement should not be a big deal.

I have tried two approaches to solve for sampling without replacement (as this seemed to be the easier case due to some terms cancelling out), but got stuck with both. The likelihood when sampling without replacement is:

$$ 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 \choose z_n}}{ {\sum_{n=1}^N y_n} \choose k} \prod_{n=1}^N Poisson(y_n |\lambda_n) \right) }{ P(\sum_n y_n \geq k|\Lambda) } $$

EDIT: The "attempted solutions section was removed as the solution in the answer does not build on them (and is way better)

Answer 1

The case without replacement

If you have $n$ independent Poisson distributed variables

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

simulation

##### 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)

results

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

Negative binomial

The arguments would work the same for the case of a negative binomial distribution which results (under certain conditions) into a Dirichlet-multinomial distribution.

Below is an example simulation

##### 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