Afsnit 3.5: Opgave med besvarelse: levetid af overfladebobler

I artiklen How does the presence of stevia glycosides impact surface bubbles stability? undersøges betydningen af af erstatte sukker med sødestoffer i drikkevarer med hensyn til de bobler, der dannes på overfladen. Specielt laves et eksperiment, hvor man kan måle levetiden af en enkelt boble. Data, vi betragter, er levetiden for 575 bobler i en opløsning med sukker uden sødestoffer. Data er simulerede i overensstemmelse med det "fintmaskede" histogram i figur 6 i artiklen. I den samme figur kan I også se forskellen til situationer, hvor sukker erstattes af sødestoffer. Et tæthedshistogram for data med sukkeropløsning er vist i den følgende figur, og data er indskrevet i kodevinduet i det skjulte punkt 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 python med kommandoerne weibull_\text{\textunderscore}min.pdf(x,α\alpha,λ\lambda) og weibull_\text{\textunderscore}min.cdf(x,α\alpha,λ\lambda). Til at lave goodness of fit testet skal du benytte en intervalinddeling med intervaller af længde 2.5 startende i nul. Desuden skal der bruges, at maksimum likelihood skønnene baseret på antallene i de forskellige intervaller er α^=1.5113\hat\alpha=1.5113 og λ^=14.32885.\hat\lambda=14.32885.
Eksempel 3.5.1. (Besvarelse)
Idet den største værdi i data er 45.1, laver vi intervalinddelingen (0,2.5],(2.5,5.0],…,(42.5,45.0],(45.0,∞).(0,2.5],(2.5,5.0],\ldots,(42.5,45.0],(45.0,\infty). Antallene i de forskellige intervaller betegnes (a1,…,a19)(a_1,\ldots,a_{19}) og findes i python med kommandoen hist, se kodevinduet nedenfor. For de tilhørende stokastiske variable bruges Statistik Model 3.1.1,
(A1,…,A19)∼multinom(575,(π1,…,π19)),πj≥0,∑jπj=1. (A_1,\ldots,A_{19})\sim\text{multinom}(575,(\pi_1,\ldots,\pi_{19})), \quad \pi_j\geq 0,\enspace\sum_j\pi_j=1.
Vi ønsker at teste hypotesen
πj=Fw(2.5j,α,λ)−Fw(2.5(j−1),α,λ),j=1,…,18,π19=1−Fw(45.0,α,λ),α,λ>0,\begin{aligned} \pi_j& =F_w(2.5j,\alpha,\lambda)-F_w(2.5(j-1),\alpha,\lambda),\quad j=1,\ldots,18, \\ \pi_{19}&=1-F_w(45.0,\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 α^=1.5113\hat\alpha=1.5113 og λ^=14.3285.\hat\lambda=14.3285. De forventede kan derfor beregnes som
ej=575⋅(Fw(2.5j,1.5113,14.3285)−Fw(2.5(j−1),1.5113,14.3285)),j=1,…,18,e19=575⋅(1−Fw(45.0,1.5113,14.3285)).\begin{aligned} e_j&=575\cdot \big( F_w(2.5j,1.5113,14.3285)-F_w(2.5(j-1),1.5113,14.3285)\big), \quad j=1,\ldots,18, \\ e_{19}&=575\cdot (1-F_w(45.0,1.5113,14.3285)). \end{aligned}
Fra beregningen i kodevinduet får vi de observerede (første række) og forventede (anden række med 1 decimal):

Kasse   1    2    3    4    5    6    7    8    9   10   11   12   13  14  15  16  17  18  19
Obs    40   64   72   80   59   51   61   42   32   19    7   15   16   7   3   2   1   2   2
Forv 39.7 66.3 74.2 73.1 66.8 58.0 48.3 38.8 30.3 23.0 17.1 12.4  8.8 6.1 4.2 2.8 1.9 1.2 2.0

For at få alle de forventede større end eller lig med 5 slås kasse 15 og 16 sammen såvel som kasserne 17, 18 og 19. Dette giver de observerede antal 4 og 5, og de forventede antal 7.0 og 5.1. Efter denne sammenlægning er der 16 kasser, hvorfor antallet af frihedsgrader i χ2\chi^2-fordelingen bliver 16-1-2=13, 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 kodevinduet viser, at G=21.51,G=21.51, og den tilhørende pp-værdi er P(G≥13.72)=1−χcdf2(21.51,13)=0.063.P(G\geq 13.72)=1-\chi^2_{\text{cdf}}(21.51,13)=0.063. Da pp-værdien ligger lidt over 0.05, siger vi at data ikke strider mod hypotesen om, at levetiderne af overfladeboblerne er weibullfordelt.

3.5.2 Beregning i python af goodness of fit

Koden nedenfor indeholder en del linjer, der er specifikke for det konkrete eksempel, der betragtes. Alle steder, hvor weibullfordelingen bruges, er naturligvis specifikke, såvel som intervalinddeling og sammenlægning af intervaller for at opnå, at de forventede er store nok. Disse steder er markeret i kommentarerne med "eksempel". Udover beregningen af goodness of fit testet indtegnes også weibulltætheden i histogrammet. Bemærk i koden brug af np.delete og np.append.

Se opstartskoden (til/fra)

ForegåendeNæste