Schätzung der Populationsgröße anhand der Häufigkeit der Stichproben von Duplikaten und Unikaten

Es gibt einen Webservice, über den ich Informationen zu einem zufälligen Artikel anfordern kann. Für jede Anfrage hat jeder Artikel die gleiche Chance, zurückgeschickt zu werden.

Ich kann weiterhin Artikel anfordern und die Anzahl der Duplikate und Unikate aufzeichnen. Wie kann ich diese Daten verwenden, um die Gesamtzahl der Artikel zu schätzen?

probability population coupon-collector-problem hoju
quelle

Was Sie schätzen möchten, ist nicht eine Stichprobengröße, sondern die Größe einer Grundgesamtheit (Gesamtzahl der vom Webservice zurückgegebenen eindeutigen Elemente).

GaBorgulya

Antworten:

Dies ist im Wesentlichen eine Variante des Problems des Gutscheinsammlers.

Wenn insgesamt $n$ Elemente vorhanden sind und Sie eine Stichprobengröße $s$ mit Ersatz genommen haben, ist die Wahrscheinlichkeit, $u$ eindeutige Elemente identifiziert zu haben ,

P r (U = u | n, s) = \frac{S_{2} (s, u) n!}{(n - u)! n^{s}}

$Pr(U=u|n,s) = \frac{S_2(s,u) n! }{ (n-u)! n^s }$ wogibtStirling Zahlen der zweiten Art

S_{2} (s, u)

$S_2(s,u)$

Alles was Sie jetzt brauchen , ist eine vorherige Verteilung für $Pr(N=n)$ , gelten Bayes - Theorem und eine hintere Verteilung erhalten $N$ .

Henry
quelle

Dies scheint einige Informationen zu verlieren, da die Häufigkeit, mit der Elemente 2, 3, 4, ... Mal beobachtet wurden, nicht berücksichtigt wird .

whuber

@whuber: Möglicherweise werden die Informationen nicht verwendet. Wenn Sie jedoch weitere Untersuchungen durchführen, sollten Sie feststellen, dass die Anzahl der eindeutigen Elemente eine ausreichende Statistik darstellt. Wenn Sie beispielsweise eine Stichprobe mit dem Ersatz von 4 Elementen aus einer Grundgesamtheit von

, beträgt die Wahrscheinlichkeit, dass 3 von einem Element und 1 von einem anderen erhalten werden,

n

$n$

, dass 2 jeweils zwei Punkte zu bekommen, egalwas

ist, so die detaillierten Frequenzenwissenüber die Bevölkerung nicht mehr nützliche Informationen gibtals nurwissenin der Probe gefunden zwei Unikate waren.

\frac{4}{3}

$\frac{4}{3}$

n

$n$

Henry

Interessanter Punkt zur ausreichenden Anzahl von Unikaten. Die Frequenzen können also zur Überprüfung der Annahmen (von Unabhängigkeit und gleicher Wahrscheinlichkeit) dienen, sind aber ansonsten nicht erforderlich.

whuber

Ich habe bereits einen Vorschlag gemacht, der auf Stirling-Zahlen der zweiten Art und Bayes'schen Methoden basiert.

Für diejenigen, die Stirling-Zahlen zu groß oder Bayes'sche Methoden zu schwierig finden, könnte eine gröbere Methode sein

E [U | n, s] = n (1 - {(1 - \frac{1}{n})}^{s})

$E[U|n,s] = n\left( 1- \left(1-\frac{1}{n}\right)^s\right)$

v a r [U | n, s] = n {(1 - \frac{1}{n})}^{s} + n^{2} (1 - \frac{1}{n}) {(1 - \frac{2}{n})}^{s} - n^{2} {(1 - \frac{1}{n})}^{2 s}

$var[U|n,s] = n\left(1-\frac{1}{n}\right)^s + n^2 \left(1-\frac{1}{n}\right)\left(1-\frac{2}{n}\right)^s - n^2\left(1-\frac{1}{n}\right)^{2s}$

und mit numerischen Methoden zurückrechnen.

Zum Beispiel mit GaBorgulya des am Beispiel und ein beobachtetes , dies könnte uns eine Schätzung geben für die Bevölkerung. $s=300$ $U = 265$ $\hat{n} \approx 1180$

Wenn das die Population gewesen wäre, hätte es uns eine Varianz für von ungefähr 25 gegeben, und zwei willkürliche Standardabweichungen auf beiden Seiten von 265 wären ungefähr 255 und 275 (wie gesagt, dies ist eine grobe Methode). 255 hätte uns eine Schätzung für ungefähr 895 gegeben, während 275 ungefähr 1692 gegeben hätte. Die 1000 des Beispiels liegt bequem innerhalb dieses Intervalls. $U$ $n$

Henry
quelle

(+1) Es ist interessant zu bemerken, dass wenn das Verhältnis

sehr klein ist, im Wesentlichen keine Informationen über

vorliegen können und man daher nicht erwarten kann,

sehr gut abzuschätzen . Wenn

sehr groß ist, ist

ein guter Schätzer. Wir brauchen also etwas, das im mittleren Bereich funktioniert.

s / n

$s/n$

n

$n$

n

$n$

s / n

$s/n$

U

$U$

Kardinal

Auch

wobei

ist die Taylorreihennäherung

ter Ordnung an

. Mit

1 - (1 - 1 / n)^{s} \geq (1 - f_{k} (s / n)) / f_{k} (s / n)

$1 - (1-1/n)^s \geq (1-f_k(s/n)) / f_k(s/n)$

f_{k} (x) = \sum_{i = 0}^{k} x^{i} / i!

$f_k(x) = \sum_{i=0}^k x^i/i!$

k

$k$

e^{x}

$e^x$

gibt einen Schätzer

k = 1

$k=1$

. Eine Kontinuitätskorrektur für kleine

kann durch Hinzufügen einer Konstanten (wie 1) im Nenner erhalten werden. Dieser Schätzer ist nicht so gut für das Beispiel als numerisch zur Lösung

, wie Sie getan haben, though.

\tilde{n} = \frac{s}{s - U} U

$\tilde{n} = \frac{s}{s-U} U$

s

$s$

\hat{n}

$\hat{n}$

Kardinal

Sie können die Capture-Recapture-Methode verwenden , die auch als Rcapture R-Paket implementiert ist .

Hier ist ein Beispiel, das in R codiert ist. Nehmen wir an, dass der Webdienst N = 1000 Elemente enthält. Wir werden n = 300 Anfragen stellen. Erzeugen Sie eine Zufallsstichprobe, wobei die Elemente von 1 bis k nummeriert werden, wobei k die Anzahl der verschiedenen Elemente ist, die wir gesehen haben.

N = 1000; population = 1:N # create a population of the integers from 1 to 1000
n = 300 # number of requests
set.seed(20110406)
observation = as.numeric(factor(sample(population, size=n,
  replace=TRUE))) # a random sample from the population, renumbered
table(observation) # a table useful to see, not discussed
k = length(unique(observation)) # number of unique items seen
(t = table(table(observation)))

Das Ergebnis der Simulation ist

  1   2   3 
234  27   4

Unter den 300 Anfragen gab es also 4 Artikel, die dreimal gesehen wurden, 27 Artikel, die zweimal gesehen wurden und 234 Artikel, die nur einmal gesehen wurden.

Schätzen Sie nun N aus dieser Stichprobe:

require(Rcapture)
X = data.frame(t)
X[,1]=as.numeric(X[,1])
desc=descriptive(X, dfreq=TRUE, dtype="nbcap", t=300)
desc # useful to see, not discussed
plot(desc) # useful to see, not discussed
cp=closedp.0(X, dfreq=TRUE, dtype="nbcap", t=300, trace=TRUE)
cp

Das Ergebnis:

Number of captured units: 265 

Abundance estimations and model fits:
                  abundance       stderr      deviance   df           AIC
M0**                  265.0          0.0  2.297787e+39  298  2.297787e+39
Mh Chao              1262.7        232.5  7.840000e-01    9  5.984840e+02
Mh Poisson2**         265.0          0.0  2.977883e+38  297  2.977883e+38
Mh Darroch**          553.9         37.1  7.299900e+01  297  9.469900e+01
Mh Gamma3.5**  5644623606.6  375581044.0  5.821861e+05  297  5.822078e+05

 ** : The M0 model did not converge
 ** : The Mh Poisson2 model did not converge
 ** : The Mh Darroch model did not converge
 ** : The Mh Gamma3.5 model did not converge
Note: 9 eta parameters has been set to zero in the Mh Chao model

So wird nur das Mh Chao Modell konvergiert, es wird geschätzt 1262,7 =. $\hat{N}$

BEARBEITEN: Um die Zuverlässigkeit der obigen Methode zu überprüfen, habe ich den obigen Code für 10000 generierte Samples ausgeführt. Das Mh Chao-Modell konvergierte jedes Mal. Hier ist die Zusammenfassung:

> round(quantile(Nhat, c(0, 0.025, 0.25, 0.50, 0.75, 0.975, 1)), 1)
    0%   2.5%    25%    50%    75%  97.5%   100% 
 657.2  794.6  941.1 1034.0 1144.8 1445.2 2162.0 
> mean(Nhat)
[1] 1055.855
> sd(Nhat)
[1] 166.8352

GaBorgulya
quelle

Es scheint eine Rechtfertigung für die Verwendung von Capture-Recapture-Modellen erforderlich zu sein, da dies kein Standard-Capture-Recapture-Experiment ist. (Möglicherweise kann es als 300 Capture-Ereignisse angesehen werden, aber der Aufruf von closedp scheint dies nicht anzuzeigen.)

whuber

@whuber Ja, ich habe das Beispiel als 300 Capture-Ereignisse angesehen. Wie meinst du das "der Anruf bei closedp scheint das nicht anzuzeigen"? Ich schätze (konstruktive) Kritik und bin gerne bereit, meine Antwort zu korrigieren (oder gegebenenfalls zu löschen), wenn sich herausstellt, dass sie falsch ist.

GaBorgulya

Dies scheint ein vernünftiger Ansatz zu sein. Allerdings werde ich R nicht verwenden, daher muss ich die Mathematik dahinter verstehen. Die Wiki-Seite behandelt eine Situation mit zwei Ereignissen - wie wende ich sie auf diesen Fall an?

Hoju

@Ga Ich verstehe: Sie haben eine 300 x 300-Matrix für die Daten erstellt! Die Ineffizienz dieses Codes hat mich getäuscht: Es wäre einfacher und direkter, `closedp.0 (Y, dfreq = TRUE, dtype =" nbcap ", t = 300) 'zu verwenden, wobei Y die Frequenzmatrix ist {{1,234}, {2,27}, {3,4}} (die Sie zweimal berechnet und tatsächlich angezeigt haben!) Genauer gesagt sind die Konvergenzfehler alarmierend, was darauf hindeutet, dass es Probleme mit dem zugrunde liegenden Code oder den zugrunde liegenden Modellen gibt. (Eine erschöpfende Suche in den Dokumenten nach "M0"

ergibt

@whuber Ich habe den Code nach Ihrem Vorschlag vereinfacht (dfreq = TRUE, dtype = "nbcap", t = 300). Danke noch einmal.

GaBorgulya