Analytische Lösung der Probenahme mit oder ohne Ersatz nach Poisson / Negativ-Binomial

Kurzfassung

Ich versuche, die zusammengesetzte Wahrscheinlichkeit, die sich aus unabhängigen Poisson-Ziehungen und weiteren Stichproben mit oder ohne Ersatz ergibt, analytisch zu lösen / zu approximieren (es ist mir eigentlich egal, welche). Ich möchte die Wahrscheinlichkeit mit MCMC (Stan) verwenden, daher benötige ich die Lösung nur bis zu einer konstanten Laufzeit. Letztendlich möchte ich einen Prozess modellieren, bei dem die ersten Ziehungen von neg stammen. Binomialverteilung, aber ich denke, ich werde mit einer Lösung für den Poisson-Fall dorthin gelangen können.

Es ist gut möglich, dass die Lösung nicht durchführbar ist (ich verstehe die Mathematik nicht genug, um feststellen zu können, ob dies ein einfaches oder ein sehr schwieriges Problem ist). Ich interessiere mich daher auch für Annäherungen, negative Ergebnisse oder Intuition, warum das Problem wahrscheinlich unlösbar ist (z. B. im Vergleich zu einem bekannten schwierigen Problem). Links zu nützlichen Artikeln / Theoremen / Tricks, die mir helfen, voranzukommen, sind gute Antworten, auch wenn ihre Verbindung zum vorliegenden Problem nicht vollständig geklärt ist.

Formelle Stellungnahme

Formal wird zuerst unabhängig gezeichnet, und dann probiere ich zufällig Elemente aus ganz , um . Dh ich ziehe farbige Kugeln aus einer Urne, wobei die Anzahl der Kugeln der Farbe aus . Hier wird als bekannt und fest angenommen und wir bedingen . Technisch gesehen erfolgt die Probenahme ersatzlos, aber die Annahme, dass die Probenahme durch Ersatz erfolgt, sollte keine große Sache sein.k Y Z = ( z 1 , . . . , Z N ) k n P o i s ( λ n ) k ∑ n y n ≥ $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$

Ich habe zwei Ansätze ausprobiert, um die ersatzlose Probenahme zu lösen (da dies aufgrund der Aufhebung einiger Begriffe der einfachere Fall zu sein schien), blieb aber bei beiden hängen. Die Wahrscheinlichkeit bei ersatzloser Probenahme ist:

BEARBEITEN: Der Abschnitt "Lösungsversuche" wurde entfernt, da die Lösung in der Antwort nicht darauf aufbaut (und viel besser ist).

Der Fall ohne Ersatz

Wenn Sie unabhängige Poisson-verteilte Variablen haben$n$

und Bedingung auf

dann

Sie können also Ihre Urne mit den fachen farbigen Kugeln füllen, zuerst den Wert für das gesamte (das ist der Poisson-verteilte Cutoff durch die Bedingung ) zeichnen und dann die Urne mit Kugeln gemäß der Multinomialverteilung füllen .$n$$Y_i$$K$$K \geq k$$K$

Diese Füllung der Urne mit Kugeln gemäß einer multinomialen Verteilung entspricht dem Zeichnen für jede Kugel unabhängig von der Farbe von der kategorialen Verteilung. Dann können Sie die ersten betrachten Kugeln , die die Stichprobe als Definition in die Urne hinzugefügt wurden (wenn diese Probe gezogen wird , ohne Ersatz) und die Verteilung für diese ist nur ein weiteres multinomial verteilt Vektor: $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)

Ergebnisse

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

Negatives Binomial

Die Argumente würden für den Fall einer negativen Binomialverteilung, die ( unter bestimmten Bedingungen ) zu einer Dirichlet-Multinomialverteilung führt, gleich funktionieren .

Unten finden Sie eine Beispielsimulation

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