Afsnit 3.5: Opgave med besvarelse: vindhastigheder

Data i denne opgave består af den daglige middelvind i Tirstrup gennem hele 2019. Data er hentet hos Iowa Environmental Mesonet, og de daglige middelvinde er givet i kilometer per time. Et tæthedshistogram er vist i nedenstående figur, og data er indskrevet i kodevinduet nedenfor.
Data af denne type beskrives ofte med weibullfordelingen, og i opgaven her skal der laves et goodness of fit test for, om weibullfordelingen beskriver data. Hvis den stokastiske variabel XX er weibullfordelt, gælder der
P(X>x)=e−(x/λ)α,x≥0, P(X>x)=e^{-(x/\lambda)^\alpha},\quad x\geq 0,
hvor α>0\alpha>0 kaldes en formparameter og λ>0\lambda>0 en skalaparameter. Tæthedsfunktionen og fordelingsfunktionen for en weibullfordeling beregnes i R med kommandoerne dweibull(x,α\alpha,λ\lambda) og pweibull(x,α\alpha,λ\lambda). Til at lave goodness of fit testet skal der benyttes en intervalinddeling med intervaller af længde 3 startende i nul. Desuden skal der bruges, at maksimum likelihood skønnene baseret på antallene i de forskellige intervaller er α^=3.8851\hat\alpha=3.8851 og λ^=36.1116.\hat\lambda=36.1116.
Eksempel 3.5.1. (Besvarelse)
Idet den største værdi i data er 54, laver vi intervalinddelingen (0,3],(3,6],…,(48,51],(51,∞).(0,3],(3,6],\ldots,(48,51],(51,\infty). Antallene i de forskellige intervaller betegnes (a1,…,a18)(a_1,\ldots,a_{18}) og findes i R med kommandoen hist(vind,breaks=c(0:18)*3)$\text{\textdollar}counts, hvor vind er en vektor med de målte værdier. For de tilhørende stokastiske variable bruges Statistik Model 3.1.1,
(A1,…,A18)∼multinom(365,(π1,…,π18)),πj≥0,∑jπj=1. (A_1,\ldots,A_{18})\sim\text{multinom}(365,(\pi_1,\ldots,\pi_{18})), \quad \pi_j\geq 0,\enspace\sum_j\pi_j=1.
Vi ønsker at teste hypotesen
πj=Fw(3j,α,λ)−Fw(3(j−1),α,λ),j=1,…,17,π18=1−Fw(51,α,λ),α,λ>0,\begin{aligned} \pi_j& =F_w(3j,\alpha,\lambda)-F_w(3(j-1),\alpha,\lambda),\quad j=1,\ldots,17, \\ \pi_{18}&=1-F_w(51,\alpha,\lambda),\quad \alpha,\lambda>0, \end{aligned}
hvor Fw(x,α,λ)F_w(x,\alpha,\lambda) er fordelingsfunktionen for en weibullfordeling. Fra opgaveformuleringen vides, at skønnene over de ukendte parametre er α^=3.8851\hat\alpha=3.8851 og λ^=36.1116.\hat\lambda=36.1116. De forventede kan derfor beregnes som
ej=365⋅(Fw(3j,3.8851,36.1116)−Fw(3(j−1),3.8851,36.1116)),j=1,…,17,e18=1−Fw(51,3.8851,36.1116).\begin{aligned} e_j&=365\cdot \big( F_w(3j,3.8851,36.1116)-F_w(3(j-1),3.8851,36.1116)\big), \quad j=1,\ldots,17, \\ e_{18}&=1-F_w(51,3.8851,36.1116). \end{aligned}
Fra R-beregningen får vi de observerede (første række) og forventede (anden række):

  0   0   0   1   9   12   22   26   37   45   37   38   36   35   33   15  13   6
0.0 0.3 1.3 3.4 6.8 11.8 18.2 25.7 33.3 39.8 44.0 44.5 41.1 34.3 25.7 17.0 9.8 8.0

For at få alle de forventede større end eller lig med 5 slås de fire første kasser sammen. Dette giver det observerede antal 1 og det forventede antal 5.02. Efter denne sammenlægning er der 15 kasser, hvorfor antallet af frihedsgrader i χ2\chi^2-fordelingen bliver 15-1-2=12, idet vi under hypotesen har to frie parametre (α\alpha og λ\lambda). GG-teststørrelsen for vores hypotese beregnes fra formlen G=2∑ja~jlog⁡(a~j/e~j),G=2\sum_j\tilde a_j\log(\tilde a_j/\tilde e_j), hvor a~j\tilde a_j og e~j\tilde e_j er de observerede og forventede, efter at kasser er slået sammen. Beregningen i R viser, at G=13.72,G=13.72, og den tilhørende pp-værdi er P(G≥13.72)=1−χcdf2(13.72,12)=0.32.P(G\geq 13.72)=1-\chi^2_{\text{cdf}}(13.72,12)=0.32. Da pp-værdien ligger langt over 0.05, strider data ikke mod hypotesen om, at de daglige middelvinde er weibullfordelt.

3.5.2 Beregning i R af goodness of fit

I R-kørslen nedenfor har jeg også indtegnet weibulltætheden i histogrammet. Desuden binder jeg de forskellige dele af output sammen ved at bruge R-kommandoen list.

Se opstartskoden (til/fra)

ForegåendeNæste