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

8

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 nkY=(y1,...,yN),ynPois(λn)kYZ=(z1,...,zN)knPois(λn)knynk

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:

P(Z=(z1,...,zN)|Λ=(λ1,...,λN))=Y;n:ynzn(n=1N(ynzn)(n=1Nynk)n=1NPoisson(yn|λn))P(nynk|Λ)

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

Martin Modrák
quelle

Antworten:

3

Der Fall ohne Ersatz

Wenn Sie unabhängige Poisson-verteilte Variablen habenn

YiPois(λi)

und Bedingung auf

j=inYj=K

dann

{Yi}Multinom(K,(λij=1nλj))

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 .nYiKKkK

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: Kk{Zi}

{Zi}Multinom(k,(λij=1nλj))

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
Sextus Empiricus
quelle
1
Danke für die Antwort. Ich habe diesen Ansatz bereits ausprobiert und habe zwei Probleme: a) Ich verstehe nicht, warum die Konditionierung der Summe mit dem Resampling (mathematisch) identisch ist, und b) meine (zugegebenermaßen begrenzten) simulierten Daten aus Poisson + Resampling konnten nicht gut angepasst werden mit multinomialer Verteilung. Könnten Sie näher erläutern, warum Konditionierung und Resampling gleich sind?
Martin Modrák
Oh, bei weiterer Lektüre denke ich, dass ich Ihren Standpunkt verstehe. Vielleicht habe ich gerade einen blöden Fehler mit den Simulationen gemacht, werde nachsehen.
Martin Modrák
Ich benutze zwei Schritte. 1) Das Füllen der Urne mit farbigen Kugeln durch Ziehen von aus unabhängigen Poisson-Verteilungen entspricht dem Füllen der Urne durch Ziehen einer Gesamtsumme aus der Poisson-Verteilung und dann des aus einer Multinomialverteilung. 2) Das Zeichnen einer Teilmenge von Kugeln aus einer Urne, die mit farbigen Kugeln gefüllt ist, deren Anzahl durch eine Multinomialverteilung definiert ist, entspricht einer weiteren Multiniomialverteilung (Intuition kann man sehen, dass die Urne für die Farbe jeder Kugel, aus der eine Ziehung gezogen wird, gefüllt ist die kategoriale Verteilung). YiYiYi=KYi
Sextus Empiricus
Ja, es gab einen Fehler in meinem Simulationscode :-) Danke für die Hilfe! Sie haben Recht und ich verstehe auch den Schritt in Ihrer Argumentation, den ich vorher nicht gesehen habe. Ich versuche zu verstehen, warum die gleiche Beziehung nicht für neg gilt. Binomial- und Dirichlet-Multinomialverteilung. (Unter der Annahme, dass die NB-Variablen konstante Fano-Faktoren haben, abhängig von ihrer Summe erhalte ich DM) - weil "Abtastung ohne Ersatz durch Multinomial" multinomial ist, aber "Abtastung ohne Ersatz durch DM" nicht DM ist (ist diese Argumentation Ihrer Ansicht nach richtig? )
Martin Modrák
1
Es stellte sich heraus, dass ich auch die zweite Simulation durcheinander gebracht habe. Es scheint auch für den DM-Fall zu funktionieren.
Martin Modrák