Résolution analytique de l'échantillonnage avec ou sans remplacement après binôme Poisson / Négatif

Version courte

J'essaie de résoudre / approximer analytiquement la probabilité composite qui résulte de tirages de Poisson indépendants et d'un échantillonnage supplémentaire avec ou sans remplacement (je ne me soucie pas vraiment lequel). Je veux utiliser la vraisemblance avec MCMC (Stan), donc je n'ai besoin de la solution que jusqu'à un terme constant. En fin de compte, je veux modéliser un processus où les tirages initiaux proviennent de neg. distribution binomiale, mais je pense que je pourrai y arriver avec une solution pour le cas de Poisson.

Il est fort possible que la solution ne soit pas réalisable (je ne comprends pas suffisamment les mathématiques pour pouvoir dire s'il s'agit d'un problème simple ou très difficile). Je m'intéresse donc aussi aux approximations, aux résultats négatifs ou à l'intuition pour laquelle le problème est probablement insoluble (par exemple en le comparant à un problème difficile connu). Les liens vers des articles / théorèmes / astuces utiles qui m'aideront à avancer sont de bonnes réponses même si leur lien avec le problème en question n'est pas entièrement établi.

Déclaration officielle

Plus formellement, d'abord est dessiné indépendamment, puis éléments au hasard de tout pour obtenir . C'est-à-dire que je tire boules colorées d'une urne où la quantité de boules de couleur est tirée de . Ici, est supposé connu et fixe et nous conditionnons sur . Techniquement, l'échantillonnage se fait sans remplacement, mais en supposant que l'échantillonnage avec remplacement ne devrait pas être un gros problème.$Y = (y_1, ..., y_N), y_n \sim Pois(\lambda_n)$$k$$Y$$Z = (z_1,...,z_N)$$k$$n$$Pois(\lambda_n)$$k$$\sum_n y_n \geq k$

J'ai essayé deux approches pour résoudre l'échantillonnage sans remplacement (car cela semblait être le cas le plus facile en raison de l'annulation de certains termes), mais je suis resté coincé avec les deux. La probabilité lors de l'échantillonnage sans remplacement est:

EDIT: la section "tentatives de solutions a été supprimée car la solution dans la réponse ne s'appuie pas sur elles (et est bien meilleure)

Le boîtier sans remplacement

Si vous avez variables distribuées de Poisson indépendantes$n$

et condition sur

puis

Vous pouvez donc remplir votre urne avec les fois les boules colorées , comme dessiner d'abord la valeur du total (qui est la distribution de Poisson coupée par la condition ), puis remplir l'urne avec boules selon la distribution multinomiale.$n$$Y_i$$K$$K \geq k$$K$

Ce remplissage de l'urne avec des billes , selon une distribution multinomiale, équivaut à dessiner pour chaque balle indépendamment la couleur de la distribution catégorielle. Ensuite, vous pouvez considérer les premières boules qui ont été ajoutées à l'urne comme définissant l'échantillon aléatoire (lorsque cet échantillon est tiré sans remplacement) et la distribution pour cela n'est qu'un autre vecteur distribué multinomial: $K$$k$$\lbrace Z_i \rbrace$

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)

résultats

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

Binôme négatif

Les arguments fonctionneraient de la même manière dans le cas d'une distribution binomiale négative qui se traduit ( sous certaines conditions ) en une distribution Dirichlet-multinomiale.

Voici un exemple de 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