Wie kann man feststellen, ob Daten einer Poisson-Verteilung in R folgen?

Ich bin Student und habe ein Projekt für meine Wahrscheinlichkeitsklasse. Grundsätzlich habe ich einen Datensatz über die Wirbelstürme, die mein Land mehrere Jahre lang heimgesucht haben.

In meinem Wahrscheinlichkeitsbuch (Wahrscheinlichkeit und Statistik mit R) gibt es ein (nicht vollständiges) Beispiel dafür, wie überprüft werden kann, ob die Daten einer Poisson-Verteilung folgen. Sie versuchen zu beweisen, dass diese drei Kriterien eingehalten werden: (Aus meinem Buch, Seite 120 (Kriterien) Seite 122-123 Beispiel)

1- Die Anzahl der Ergebnisse in nicht überlappenden Intervallen ist unabhängig. Mit anderen Worten ist die Anzahl der Ergebnisse im Zeitintervall (0, t] unabhängig von der Anzahl der Ergebnisse im Zeitintervall (t, t + h], h> 0

2- Die Wahrscheinlichkeit von zwei oder mehr Ergebnissen in einem ausreichend kurzen Intervall ist praktisch Null. Mit anderen Worten ist, vorausgesetzt, h ist ausreichend klein, die Wahrscheinlichkeit, zwei oder mehr Ergebnisse in dem Intervall (t, t + h) zu erhalten, vernachlässigbar, verglichen mit der Wahrscheinlichkeit, ein oder null Ergebnisse in dem gleichen Zeitintervall zu erhalten.

3- Die Wahrscheinlichkeit, dass genau ein Ergebnis in einem ausreichend kurzen Intervall oder einer kleinen Region auftritt, ist proportional zur Länge des Intervalls oder der Region. Mit anderen Worten ist die Wahrscheinlichkeit eines Ergebnisses in einem Intervall der Länge h Lambda * h.

Kriterium 3 bleibt jedoch "als Übung".

A- Kann mir jemand sagen, ob es eine "einfachere" Möglichkeit gibt, festzustellen, ob mein Datensatz einer Poisson-Verteilung folgt?

B- Kann mir jemand Kriterium 1 und 3 mit einem Beispiel erklären (wenn es mit R ist, fantastisch)?

Vielen Dank!

Hinweis: Entschuldigung für den langen Beitrag. Außerdem muss ich die Daten konvertieren, damit ich eine Tabelle wie folgt habe:

  number of hurricanes       | 0 | 1 | 2  etc.
  -----------------------------------------
  total years that have      |   |   |
  that number of hurricanes  |   |   |

Es gibt unendlich viele Möglichkeiten, wie sich eine Verteilung geringfügig von einer Poisson-Verteilung unterscheidet. Sie können nicht erkennen, dass es sich um einen Datensatz handelt aus einer Poisson - Verteilung gezogen. Was Sie tun können, ist nach Inkonsistenzen mit dem zu suchen, was Sie mit einem Poisson sehen sollten, aber ein Mangel an offensichtlicher Inkonsistenz macht es nicht zu einem Poisson.

Wenn Sie jedoch diese drei Kriterien überprüfen, überprüfen Sie nicht, ob die Daten aus einer Poisson-Verteilung stammen (dh, indem Sie Daten betrachten), sondern indem Sie bewerten, ob der Prozess, durch den die Daten generiert werden, den Anforderungen entspricht Bedingungen eines Poisson-Prozesses; Wenn die Bedingungen alle oder fast alle eingehalten werden (und das ist eine Überlegung des Datenerzeugungsprozesses), könnten Sie etwas von einem Poisson-Prozess haben oder einem sehr nahe kommen, was wiederum eine Möglichkeit wäre, Daten zu erhalten, die von einem nahe liegenden Prozess stammen Poisson-Verteilung.

Aber die Bedingungen sind in mehrfacher Hinsicht nicht zutreffend ... und das am weitesten von der Wahrheit entfernte ist Nummer 3. Auf dieser Grundlage gibt es keinen besonderen Grund, einen Poisson-Prozess geltend zu machen, obwohl die Verstöße möglicherweise nicht so schlimm sind, dass die resultierenden Daten weit entfernt sind von Poisson.

Wir kehren also zu statistischen Argumenten zurück, die sich aus der Untersuchung der Daten selbst ergeben. Wie würden Daten zeigen, dass es sich bei der Verteilung eher um Poisson als um etwas Ähnliches handelt?

Wie eingangs erwähnt, können Sie überprüfen, ob die Daten offensichtlich nicht mit der zugrunde liegenden Verteilung von Poisson inkonsistent sind. Das bedeutet jedoch nicht, dass sie aus einem Poisson stammen (Sie können sich bereits darauf verlassen, dass dies der Fall ist) nicht).

Sie können diese Prüfung über die Prüfung der Passgenauigkeit durchführen.

Das erwähnte Chi-Quadrat ist eines davon, aber ich würde den Chi-Quadrat-Test für diese Situation nicht selbst empfehlen **; es hat eine geringe Leistung gegen interessante Abweichungen. Wenn Sie gute Leistung anstreben, werden Sie das nicht so sehen (wenn Sie sich nicht für Leistung interessieren, warum sollten Sie dann testen?). Ihr Hauptwert liegt in der Einfachheit und sie hat einen pädagogischen Wert. Abgesehen davon ist es kein konkurrenzfähiger Fitnesstest.

** In der späteren Bearbeitung hinzugefügt: Da nun klar ist, dass dies Hausaufgaben sind, steigt die Wahrscheinlichkeit, dass Sie einen Chi-Quadrat-Test durchführen , um die Daten zu überprüfen, nicht im Widerspruch zu einem Poisson. Sehen Sie sich mein Beispiel für einen Chi-Quadrat-Fit-Test an, der unterhalb des ersten Poissonness-Diagramms durchgeführt wurde

Die Leute machen diese Tests oft aus dem falschen Grund (zB weil sie sagen wollen, dass es in Ordnung ist, andere statistische Dinge mit den Daten zu machen, die davon ausgehen, dass die Daten Poisson sind). Die eigentliche Frage ist, wie schlimm das gehen könnte. ... und die Güte von Fit-Tests helfen bei dieser Frage nicht wirklich weiter. Oft ist die Antwort auf diese Frage bestenfalls eine, die unabhängig (/ fast unabhängig) von der Stichprobengröße ist - und in einigen Fällen mit Konsequenzen, die tendenziell mit der Stichprobengröße verschwinden ... während ein Test der Anpassungsgüte unbrauchbar ist kleine Stichproben (bei denen das Risiko von Verstößen gegen Annahmen häufig am größten ist).

Wenn Sie eine Poisson-Distribution testen müssen, gibt es ein paar sinnvolle Alternativen. Eine Möglichkeit wäre, einen Anderson-Darling-Test durchzuführen, der auf der AD-Statistik basiert, aber eine simulierte Verteilung unter der Null verwendet (um die Doppelprobleme einer diskreten Verteilung zu berücksichtigen und um Parameter abzuschätzen).

Eine einfachere Alternative könnte ein Glättungstest für die Anpassungsgüte sein. Hierbei handelt es sich um eine Sammlung von Tests, die für individuelle Verteilungen entwickelt wurden, indem die Daten unter Verwendung einer Familie von Polynomen modelliert wurden, die in Bezug auf die Wahrscheinlichkeitsfunktion in der Null orthogonal sind. Alternativen niedriger (dh interessanter) Ordnung werden getestet, indem geprüft wird, ob die Koeffizienten der Polynome über der Basis Eins von Null verschieden sind, und diese können normalerweise die Parameterschätzung behandeln, indem die Terme niedrigster Ordnung aus dem Test weggelassen werden. Es gibt so einen Test für den Poisson. Ich kann eine Referenz ausgraben, wenn Sie sie brauchen.

$n(1-r^2)$$\log(x_k)+\log(k!)$$k$

Hier ist ein Beispiel für diese Berechnung (und Darstellung) in R:

y=rpois(100,5)
n=length(y)
(x=table(y))
y
 0  1  2  3  4  5  6  7  8  9 10 
 1  2  7 15 19 25 14  7  5  1  4 

k=as.numeric(names(x))
plot(k,log(x)+lfactorial(k))

Hier ist die Statistik, die ich vorgeschlagen habe, um den Fitnesstest eines Poisson durchzuführen:

n*(1-cor(k,log(x)+lfactorial(k))^2)
[1] 1.0599

Um den p-Wert zu berechnen, müssten Sie natürlich auch die Verteilung der Teststatistik unter der Null simulieren (und ich habe nicht besprochen, wie man mit Nullzählungen innerhalb des Wertebereichs umgehen könnte). Dies sollte einen einigermaßen leistungsfähigen Test ergeben. Es gibt zahlreiche andere alternative Tests.

Hier ist ein Beispiel für ein Poissonness-Diagramm mit einer Stichprobe der Größe 50 aus einer geometrischen Verteilung (p = 0,3):

Wie Sie sehen, wird ein deutlicher Knick angezeigt, der auf Nichtlinearität hinweist

Referenzen für das Poissonness-Diagramm wären:

David C. Hoaglin (1980),
"A Poissonness Plot",
Der amerikanische Statistiker
Vol. 34, Nr. 3 (August), S. 146-149

und

Hoaglin, D. und J. Tukey (1985),
"9. Überprüfen der Form diskreter Verteilungen",
Untersuchen von
Datentabellen, Trends und Formen , (Hoaglin, Mosteller & Tukey eds)
John Wiley & Sons

Die zweite Referenz enthält eine Anpassung des Diagramms für kleine Zählungen. Sie würden es wahrscheinlich einbauen wollen (aber ich habe den Hinweis nicht zur Hand).

Beispiel für einen Chi-Quadrat-Anpassungstest:

Abgesehen von der Ausführung der Chi-Quadrat-Anpassungsgüte, wie es normalerweise in vielen Klassen zu erwarten ist (wenn auch nicht so, wie ich es tun würde):

1: Ausgehend von Ihren Daten (die ich als die Daten ansehen werde, die ich zufällig in 'y' oben generiert habe), generieren Sie die Zählungstabelle:

(x=table(y))
y
 0  1  2  3  4  5  6  7  8  9 10 
 1  2  7 15 19 25 14  7  5  1  4 

2: Berechne den erwarteten Wert in jeder Zelle unter der Annahme eines von ML angepassten Poisson:

 (expec=dpois(0:10,lambda=mean(y))*length(y))
 [1]  0.7907054  3.8270142  9.2613743 14.9416838 18.0794374 17.5008954 14.1173890  9.7611661
 [9]  5.9055055  3.1758496  1.5371112

3: Beachten Sie, dass die Endkategorien klein sind; Dies macht die Chi-Quadrat-Verteilung weniger gut als eine Annäherung an die Verteilung der Teststatistik nah, aber der allgemeine Ansatz kann an eine strengere Regel angepasst werden). Reduzieren Sie benachbarte Kategorien, sodass die erwarteten Mindestwerte nicht zu weit unter 5 liegen (eine Kategorie mit einem erwarteten Countdown bei 1 von mehr als 10 Kategorien ist nicht schlecht, zwei ist ziemlich grenzwertig). Beachten Sie auch, dass wir die Wahrscheinlichkeit über "10" hinaus noch nicht berücksichtigt haben. Deshalb müssen wir auch Folgendes berücksichtigen:

expec[1]=sum(expec[1:2])
expec[2:8]=expec[3:9]
expec[9]=length(y)-sum(expec[1:8])
expec=expec[1:9]
expec
sum(expec) # now adds to n

4: in ähnlicher Weise zuklappen Kategorien auf dem beobachteten:

(obs=table(y))
obs[1]=sum(obs[1:2])
obs[2:8]=obs[3:9]
obs[9]=sum(obs[10:11])
obs=obs[1:9]

$(O_i-E_i)^2/E_i$

print(cbind(obs,expec,PearsonRes=(obs-expec)/sqrt(expec),ContribToChisq=(obs-expec)^2/expec),d=4)
  obs  expec PearsonRes ContribToChisq
0   3  4.618   -0.75282      0.5667335
1   7  9.261   -0.74308      0.5521657
2  15 14.942    0.01509      0.0002276
3  19 18.079    0.21650      0.0468729
4  25 17.501    1.79258      3.2133538
5  14 14.117   -0.03124      0.0009761
6   7  9.761   -0.88377      0.7810581
7   5  5.906   -0.37262      0.1388434
8   5  5.815   -0.33791      0.1141816

$X^2 = \sum_i (E_i-O_i)^2/E_i$

(chisq = sum((obs-expec)^2/expec))
[1] 5.414413
(df = length(obs)-1-1) # lose an additional df for parameter estimate
[1] 7
(pvalue=pchisq(chisq,df))
[1] 0.3904736

Sowohl die Diagnose als auch der p-Wert zeigen hier keinen Mangel an Übereinstimmung ... was wir erwarten würden, da die von uns generierten Daten tatsächlich Poisson waren.

Bearbeiten: Hier ist ein Link zu Rick Wicklins Blog, der die Poissonness-Darstellung und Implementierungen in SAS und Matlab behandelt

http://blogs.sas.com/content/iml/2012/04/12/the-poissonness-plot-a-goodness-of-fit-diagnostic/

Edit2: Wenn ich es richtig habe, wäre das modifizierte Poissonness-Diagramm aus der Referenz von 1985 *:

y=rpois(100,5)
n=length(y)
(x=table(y))
k=as.numeric(names(x))
x=as.vector(x)
x1 = ifelse(x==0,NA,ifelse(x>1,x-.8*x/n-.67,exp(-1)))
plot(k,log(x1)+lfactorial(k))

* Sie passen auch den Achsenabschnitt an, aber ich habe das hier nicht getan. Dies hat keine Auswirkungen auf das Erscheinungsbild des Diagramms. Sie müssen jedoch aufpassen, wenn Sie etwas anderes aus der Referenz implementieren (z. B. die Konfidenzintervalle).

(Im obigen Beispiel ändert sich das Aussehen kaum von der ersten Poissonness-Darstellung.)

Führen Sie den Chi-Quadrat-Fit-Test durch. Im Falle von goodfit()Zähldaten können wir die im vcd-Paket enthaltenen verwenden. Beachten Sie, dass wenn der p-Wert größer als 0,05 ist, wir h0 nicht zurückweisen können: Der Prozess ist ein Poisson-Prozess. Oder es ist kein Poisson-Prozess.

# load the vcd package
library(vcd) ## loading vcd package

# generate two processes for test
set.seed(2014);y=rpois(200,5)
set.seed(2014);y=rnorm(100, 5, 0.3) # goodfit asks for non-negative values
# output the results
gf = goodfit(y,type= "poisson",method= "ML")
plot(gf,main="Count data vs Poisson distribution")
summary(gf)

# to automatically get the pvalue
gf.summary = capture.output(summary(gf))[[5]]
pvalue = unlist(strsplit(gf.summary, split = " "))
pvalue = as.numeric(pvalue[length(pvalue)]); pvalue

# to mannualy compute the pvalue
chisq = sum(  (gf$observed-gf$fitted)^2/gf$fitted )

df = length(gf$observed)-1-1
pvalue = pchisq(chisq,df)
pvalue