Jak sprawdzić, czy dane są zgodne z rozkładem Poissona w R?

Jestem studentem studiów licencjackich i mam projekt do mojej klasy prawdopodobieństwa. Zasadniczo mam zbiór danych o huraganach, które nawiedziły mój kraj przez szereg lat.

W mojej Księdze prawdopodobieństwa (Prawdopodobieństwo i statystyka z R) znajduje się (niekompletny) przykład, jak sprawdzić, czy dane są zgodne z rozkładem Poissona, zaczynają próbować udowodnić, że te 3 kryteria są spełnione: (Z mojej książki, strony 120 (kryteria), strona 122-123, przykład)

1- Liczba wyników w nie nakładających się przedziałach jest niezależna. Innymi słowy, liczba wyników w przedziale czasu (0, t] jest niezależna od liczby wyników w przedziale czasu (t, t + h], h> 0

2- Prawdopodobieństwo dwóch lub więcej wyników w wystarczająco krótkim odstępie wynosi praktycznie zero. Innymi słowy, pod warunkiem, że h jest wystarczająco mały, prawdopodobieństwo uzyskania dwóch lub więcej wyników w przedziale (t, t + h] jest znikome w porównaniu z prawdopodobieństwem uzyskania jednego lub zerowego wyniku w tym samym przedziale czasu.

3- Prawdopodobieństwo dokładnie jednego wyniku w wystarczająco krótkim odstępie lub małym obszarze jest proporcjonalne do długości przedziału lub regionu. Innymi słowy, prawdopodobieństwo jednego wyniku w przedziale długości h wynosi lambda * h.

Ale kryterium 3 pozostało „jako ćwiczenie”.

Odp .: Czy ktoś może mi powiedzieć, czy istnieje łatwiejszy sposób sprawdzenia, czy mój zestaw danych jest zgodny z rozkładem Poissona?

B- Czy ktoś może wyjaśnić mi kryterium 1 i 3 jakimś przykładem (jeśli jest to R, fantastycznie)?

Dzięki!

Uwaga: Przepraszamy za długi post. Ponadto muszę przekonwertować dane, aby mieć tabelę taką jak:

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

Istnieje nieskończona liczba sposobów, aby rozkład mógł się nieco różnić od rozkładu Poissona; nie można stwierdzić, że zestaw danych jest pobierany z rozkładu Poissona. Możesz poszukać niekonsekwencji z tym, co powinieneś zobaczyć z Poissonem, ale brak oczywistej niekonsekwencji nie czyni z niego Poissona.

Jednak to, o czym tu mówisz, sprawdzając te trzy kryteria, nie polega na sprawdzeniu, czy dane pochodzą z rozkładu Poissona metodami statystycznymi (tj. Na podstawie danych), ale poprzez ocenę, czy proces generowania danych spełnia warunki procesu Poissona; jeśli wszystkie warunki są utrzymywane lub prawie utrzymywane (a to bierze pod uwagę proces generowania danych), możesz mieć coś z procesu Poissona lub bardzo blisko niego, co z kolei byłoby sposobem na uzyskanie danych, które są pobierane z czegoś zbliżonego do Rozkład Poissona.

Ale warunki nie utrzymują się na kilka sposobów ... a najdalej od prawdziwości jest numer 3. Nie ma na tej podstawie szczególnego powodu, aby twierdzić, że proces Poissona jest poważny, chociaż naruszenia mogą nie być tak złe, że uzyskane dane są dalekie z Poisson.

Wracamy więc do argumentów statystycznych pochodzących z badania samych danych. W jaki sposób dane pokazują, że rozkład był Poissona, a nie coś podobnego?

Jak wspomniano na początku, możesz sprawdzić, czy dane nie są oczywiście niespójne z podstawową dystrybucją Poissona, ale to nie znaczy, że są pobierane z Poissona (możesz już mieć pewność, że są one nie).

Możesz to sprawdzić za pomocą testów dopasowania.

Wspomniany kwadrat chi jest jednym z takich, ale sam nie poleciłbym testu chi-kwadrat dla tej sytuacji **; ma niską moc w stosunku do interesujących odchyleń. Jeśli twoim celem jest mieć dobrą moc, nie zdobędziesz jej w ten sposób (jeśli nie zależy ci na mocy, dlaczego miałbyś ją testować?). Jego główna wartość polega na prostocie i ma wartość pedagogiczną; poza tym nie jest konkurencyjny jako test dopasowania.

** Dodano w późniejszej edycji: Teraz, gdy stało się jasne, że to zadanie domowe, szanse, że powinieneś wykonać test chi-kwadrat celu sprawdzenia danych nie są niespójne z Poissonem znacznie wzrosły. Zobacz mój przykład testu dobroci dopasowania chi-kwadrat wykonanego poniżej pierwszego wykresu Poissona

Ludzie często wykonują te testy z niewłaściwego powodu (np. Dlatego, że chcą powiedzieć „dlatego można zrobić coś innego statystycznego z danymi, które zakładają, że są to dane Poissona”). Prawdziwe pytanie brzmi: „jak bardzo źle mogło to pójść?” ... a trafność testów dopasowania niewiele pomaga w tym pytaniu. Często odpowiedź na to pytanie jest w najlepszym razie niezależna (/ prawie niezależna) od wielkości próbki - aw niektórych przypadkach jedna z konsekwencjami, które zwykle odchodzą od wielkości próbki ... podczas gdy test dobroci dopasowania jest bezużyteczny w przypadku małe próbki (w których ryzyko naruszenia założeń jest często największe).

Jeśli musisz przetestować rozkład Poissona, istnieje kilka rozsądnych alternatyw. Jednym z nich byłoby zrobienie czegoś podobnego do testu Andersona-Darlinga, opartego na statystyce AD, ale z wykorzystaniem symulowanego rozkładu poniżej wartości zerowej (aby uwzględnić bliźniacze problemy rozkładu dyskretnego i że musisz oszacować parametry).

Prostszą alternatywą może być Płynny Test na dobroć dopasowania - jest to zbiór testów zaprojektowanych dla poszczególnych rozkładów poprzez modelowanie danych przy użyciu rodziny wielomianów, które są ortogonalne względem funkcji prawdopodobieństwa w wartości zerowej. Testowane są alternatywy niskiego rzędu (tj. Interesujące), sprawdzając, czy współczynniki wielomianów powyżej podstawy są różne od zera, i zwykle mogą one poradzić sobie z estymacją parametrów, pomijając terminy najniższego rzędu w teście. Jest taki test dla Poissona. Mogę wykopać referencję, jeśli jej potrzebujesz.

$n(1-r^2)$$\log(x_k)+\log(k!)$ vs $k$ (patrz Hoaglin, 1980) - jako statystyka testowa.

Oto przykład tego obliczenia (i wykresu) wykonanego w 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))

Oto statystyka, którą zasugerowałem, może być wykorzystana do testu dopasowania Poissona:

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

Oczywiście, aby obliczyć wartość p, należy również zasymulować rozkład statystyki testowej poniżej wartości zerowej (i nie dyskutowałem, jak można sobie poradzić z zerami w zakresie wartości). To powinno dać dość mocny test. Istnieje wiele innych alternatywnych testów.

Oto przykład wykonania wykresu Poissona na próbce o wielkości 50 z rozkładu geometrycznego (p = 0,3):

Jak widać, wyświetla wyraźne „załamanie”, co wskazuje na nieliniowość

Odniesienia do wykresu Poissona byłyby następujące:

David C. Hoaglin (1980),
„A Poissonness Plot”,
The American Statistician
obj. 34, nr 3 (sierpień), s. 146–149

i

Hoaglin, D. i J. Tukey (1985),
„9. Sprawdzanie kształtu dyskretnych rozkładów ”,
badanie tabel danych, trendów i kształtów
rozkładów , (red. Hoaglin, Mosteller i Tukey)
John Wiley & Sons

Drugie odniesienie zawiera korektę wykresu dla małych liczb; prawdopodobnie chciałbyś to włączyć (ale nie mam odniesienia do ręki).

Przykład wykonania testu dopasowania chi-kwadrat:

Poza wykonaniem dobroci dopasowania chi-kwadrat, sposób, w jaki zwykle można się tego spodziewać w wielu klasach (choć nie tak, jakbym to zrobił):

1: zaczynając od twoich danych (które wezmę za dane, które losowo wygenerowałem w 'y' powyżej, wygeneruj tabelę zliczeń:

(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: oblicz oczekiwaną wartość w każdej komórce, przyjmując Poissona dopasowanego przez ML:

 (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: zwróć uwagę, że kategorie końcowe są małe; powoduje to, że rozkład chi-kwadrat jest mniej dobry jako przybliżenie rozkładu statystyki testowej (powszechną regułą jest oczekiwanie oczekiwanych wartości co najmniej 5, chociaż wiele artykułów wykazało, że zasada ta jest niepotrzebnie restrykcyjna; wezmę ją blisko, ale ogólne podejście można dostosować do surowszej zasady). Zwiń sąsiednie kategorie, tak aby minimalne oczekiwane wartości były co najmniej nie znacznie poniżej 5 (jedna kategoria z oczekiwanym odliczaniem w pobliżu 1 z więcej niż 10 kategorii nie jest taka zła, dwie są dość graniczne). Pamiętaj też, że nie uwzględniamy jeszcze prawdopodobieństwa przekraczającego „10”, dlatego musimy również uwzględnić:

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: podobnie, zwiń kategorie obserwowanych:

(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]

5: Wstaw do stołu (opcjonalnie) wraz z wkładem do kwadratu chi $(O_i-E_i)^2/E_i$ i resztkowy Pearson (podpisany pierwiastek kwadratowy wkładu), mogą być przydatne, gdy próbujemy zobaczyć, gdzie nie pasuje tak dobrze:

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

6: Oblicz $X^2 = \sum_i (E_i-O_i)^2/E_i$, z utratą 1df dla oczekiwanej sumy pasującej do obserwowanej sumy i 1 dodatkowej dla oszacowania parametru:

(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

Zarówno diagnostyka, jak i wartość p pokazują tutaj brak dopasowania ... czego się spodziewalibyśmy, ponieważ dane, które wygenerowaliśmy, to Poissona.

Edycja: oto link do bloga Ricka Wicklina, który omawia fabułę Poissonnessa i mówi o implementacjach w SAS i Matlabie

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

Edycja2: Jeśli mam rację, zmodyfikowany wykres Poissonnessa z referencji z 1985 roku byłby *:

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

* Właściwie dostosowują również przechwytywanie, ale nie zrobiłem tego tutaj; nie wpływa to na wygląd fabuły, ale musisz zachować ostrożność, jeśli zaimplementujesz cokolwiek innego z referencji (np. przedziały ufności), jeśli zrobisz to inaczej niż ich podejście.

(W powyższym przykładzie wygląd prawie nie zmienia się od pierwszego wykresu Poissona.)

Wykonaj test dobroci dopasowania chi-kwadrat. W przypadku danych zliczających możemy użyć goodfit()zawartych w pakiecie vcd. Zauważ, że jeśli wartość p jest większa niż 0,05, nie możemy odrzucić h0: proces jest procesem Poissona. W przeciwnym razie nie jest to proces Poissona.

# 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