Een aselecte steekproef \(X_1,\dots ,X_n\) wordt genomen uit een verdeling \(F\) met dichtheid \(f(x) = (\theta^3/2)x^{-4}e^{-\theta/x}\) voor \(x > 0\), waar \(\theta\in\Theta=(0,\infty)\) een onbekende parameter is.
De volgende feiten zijn wellicht in het volgende van nut: als \(X\sim F\) een enkele trekking uit de de verdeling \(F\) is, geldt \(\textrm{E}(X) = \theta/2\) and \(\textrm{Var}(X) = \theta^2/4\). Bovendien, geldt dat \(1/X\sim \textrm{Gamma}(3,\theta)\): de gamma verdeling met vorm parameter \(\alpha=3\) en (inverse) schaal parameter \(\lambda=\theta\), in de terminologie van Rice.
Bepaal de meest aannemelijke schatter (MLE) \(\widehat\theta\) van \(\theta\).
De log likelihood is \(\text{const} + 3 n \log\theta -4 \sum \log X_i -\theta \sum 1/X_i\), wat (als functie van \(\theta\), bij gegeven data \(X_1,\dots, X_n\)) gelijk is aan een constante plus \(3 n \log\theta -\theta \sum 1/X_i\). Eerste afgeleide naar \(\theta\) (wat ook wel de score functie heet) is \(3 n / theta - \sum 1/X_i\); tweede afgeleide (wiens negatieve de observed information heet) is \(-3 n / \theta^2\). De tweede afgeleide is overal negatief dus de log likelihood is concaaf. De score functie is nul bij \(\theta = (\sum 1/X_i)/3n\). Dit is dus de meest aannemelijke schatter.
Laat zien zien dat \(T(X_1,\dots, X_n) = \sum_{i=1}^n X_i^{-1}\) een voldoende grootheid (Engels: “sufficient statistic”) is. (Hint: gebruik het factorizatie criterium).
We hadden al gezien dat log likelihood is gelijk aan functie van de data alleen, plus function van \(\theta\) en van \(\sum_{i=1}^n X_i^{-1}\). Dus likelihood is gelijk aan functie van data keer functie van \(\theta\) en van \(\sum_{i=1}^n X_i^{-1}\). Volgens factorisatie criterium is dus \(\theta\) en van \(\sum_{i=1}^n X_i^{-1}\) een voldoende grootheid voor \(\theta\).
Een alternative schatter van theta is \(\widetilde \theta = 2 \overline X\). Bepaal de onzuiverheid, variantie en verwachte kwadratische fout (Engels: “bias”, “variance”, “mean square error”) van deze schatter.
We weten al dat voor één waarneming \(X\) uit deze verdeling, \(\textrm{E}(X) = \theta/2\) and \(\textrm{Var}(X) = \theta^2/4\). Er geldt dus dat voor de gemiddelde van een steekproef ter grootte \(n\), \(\textrm{E}(\overline X) = \theta/2\) and \(\textrm{Var}(\overline X) = \theta^2/4n\). Dus, \(\textrm{E}(2 \overline X) = \theta\) and \(\textrm{Var}(2 \overline X) = \theta^2/n\).
De schatter \(\widetilde \theta\) is dus zuiver (heeft de goede verwachtingswaarde), en als gevolg hiervan is zijn “mean square error” gelijk aan zijn variantie, oftewel \(\theta^2/n\).
Wat zijn voor grote \(n\), bij benadering, de kansverdelingen van deze twee schatters? Welke is dan te preferen?
Beide schatters zijn voor grote \(n\) bij benadering normaal verdeeld: de mle wegens standaard asymptotische theorie mle; de andere schatter wegens de centrale limiet stelling. De verdeling van \(\widetilde \theta\) is, bij benadering (voor grote \(n\)), normaal met verwachting \(\theta\) en variantie \(\theta^2/n\)).
We zagen al bij onderdeel (a) dat de “observed information” gelijk is aan \(3 n / \theta^2\). De mle heeft dus bij benadering voor grote \(n\) de normale verdeling met verwachting \(\theta\) en variantie \(\theta^2/3n\).
De MLE is (zoals we van te voren konden verwachten) stukken beter.
Hoe zou U, voor grote \(n\), de variantie van de MLE \(\widehat \theta\) schatten?
Ik zou de variantie van de MLE schatten als \(\widehat\theta^2/3n\) waarbij \(\widehat\theta\) de MLE van \(\theta\) is.
Hoe zou U een benaderend 95% betrouwbaarheids interval voor \(\theta\) berekenen?
\(\widehat\theta \pm 2 \widehat\theta/\sqrt(3n)\) is een bij benadering 95% betrouwbaarheidsinterval voor \(\theta\). Hier staat de factor “2” eigenlijk voor 1.986:
qnorm(0.975)## [1] 1.959964qnorm(0.025)## [1] -1.959964De leden van een biologische populatie kunnen ingedeeld worden in drie genetische klassen A, B en C. Volgens een bepaald genetisch model komen de twee klassen B en C even vaak voor. Maar in het algemeen zouden de drie klassen volledig willekeurige kansen \(\theta_A, \theta_B, \theta_C\) kunnen hebben, waarbij natuurlijk wel moet gelden \(\theta_A+\theta_B+ \theta_C=1\). Een aselecte steekproef van 100 individuen uit de populatie levert de volgende aantallen op voor elk van de drie klassen, in dezelfde volgorde A, B, C: 63, 23, 14.
Wat is de meest aannemelijke schatter voor de vector parameter \(\theta=(\theta_A, \theta_B, \theta_C)\) onder de nul hypothese dat het eerder genoemde genetisch model geldt, en wat is de meest aannemelijke schatter voor \(\theta\) in het algemeen?
Zonder beperkingen op de drie kansen (dus “in het algemeen”) is de log likelihood gelijk aan \(63 \log \theta_A + 23 \log \theta_B + 14 \log \theta_C\), waarbij \(\theta_A+\theta_B+\theta_C=1\). Dit wordt gemaximaliseerd door \(\widehat \theta = (0.63, 0.23, 0.14)\).
Volgens het speciale genetisch model is \(\theta_B = \theta_C\). Definieer \(\theta_{B+C}:=\theta_B +\theta_C\). De log likelihood is nu gelijk aan \(63 \log\theta_A + 37 \log\theta_{B+C}\) waar \(\theta_A +\theta_{B+C} = 1\). Dit wordt gemaxialiseerd door \((\widehat \theta_A, \widehat \theta_{B+C}) = (0.63, 0.37)\). Dus in termen van de originele drie parametrisatie, \(\widehat\theta_0 = (0.63, 0.185, 0.185)\).
Geef een formule voor de toetsingsgrootheid \(-2 \log \Lambda\) in termen van de data en de eerder afgeleide schattingen. Hier is \(\Lambda\) de generalized likelihood ratio voor het toetsen van de eerder genoemde nul hypothese tegen het alternatief (binnen het algemene context van de opgave) dat de nul hypothese niet geldt. (U hoeft hierbij de numerieke waarde van \(- 2 \log \Lambda\) niet uit te rekenen als U geen zakrekenmachine bij de hand hebt: de opgave vraagt om een eenvoudige formule waarmee men dat in principe zou kunnen doen.)
We weten dat voor dit soort probleem, \(-2 \log \Lambda = 2 \sum O_i \log(O_i / E_i)\) is waar de index \(i\) de drie waardes \(A\), \(B\), \(C\) aanneemt; de \(O_i\) zijn de waargenomen aantallen 63, 23, 14; de \(E_i\) zijn de geschatte verwachte aantallen onder het model van de nulhypothese. Deze waardes zijn dus 63, 18.5, 18.5. Er volgt hieruit dat \(-2 \log\Lambda = 2\times(63 \log (63/63) + 23 \log(23/18.5) + 14 \log(14/18.5))\).
chisq <- 2*(63 * log (63/63) + 23 * log(23/18.5) + 14 * log(14/18.5))
chisq## [1] 2.211305pvalue <- pchisq(chisq, df = 1, lower.tail = FALSE)
pvalue## [1] 0.1370027Ik vond een waarde 2.211 voor \(- 2 \log \Lambda\). Als ik deze toets wil uitvoeren met significantie niveau \(\alpha = 10\%\), moet ik wel of niet de nul hypothese verwerpen? (Hint: appendix B van het boek van Rice bevat, o.m., tabellen met quantielen/percentielen van verschillende standaard verdelingen.)
Ik vind vandaag ook weer 2.211305. De bijhorende p-waarde is 0.1370027. Hier, vergelijk ik de waargenomen waarde van de toetsingsgrootheid met de chi-kwadraat verdeling met één vrijheidsgraad (verschil aantal parameters in nul hypothese, algemene situatie). De p-waarde (kans op zó’n extreme waarde als we hebben waargenomen) is groter dan 10% dus kunnen we de nul-hypothese niet verwerpen.
Hetzelfde conclusie vinden we, uiteraard, als we de waarde van de toetsingsgrootheid vergelijken met de 90% procentiel van de chi-kwadraat(1) verdeling:
qchisq(0.90, df = 1)## [1] 2.705543chisq > qchisq(0.90, df = 1)## [1] FALSEDefinieer \(f(x;\theta)=\theta x^{\theta-1}\) voor \(0\le x\le 1\), \(f(x)=0\) voor \(x<0\) en voor \(x>1\), waar \(\theta>0\) een vaste parameter is.
Laat zien dat de bijbehorende cumulatieve verdelingsfunctie gelijk is aan \(F(x;\theta)=x^\theta\) voor \(0\le x\le 1\); \(F(x;\theta)=0\) voor \(x\le 0\); \(F(x;\theta)=1\) voor \(x\ge 1\). Welke bekende kansverdeling stelt dit voor, in het geval dat \(\theta=1\)?
Voor \(x\in [0, 1]\) geldt \(F(x; \theta) = \int_0^x \theta x^{\theta-1}\) d x = x^$. Omdat \(\theta > 0\) is deze functie stijgend, gelijk aan 0 in \(x = 0\) en gelijk aan 1 in \(x = 1\). Als \(\theta = 1\) hebben we te maken met de uniforme verdeling op het interval \([0, 1]\).
Laat zien dat de quantiel functie \(Q(p;\theta)=p^{\frac 1 \theta}\) is (\(0\le p \le 1\)).
De oplossing van de vergelijking \(p = x^\theta\) is $Q(p; ) = x = \(p^{1/\theta}\).
Voor gegeven \(\theta\), hoe zou u een trekking uit deze verdeling op een computer kunnen nabootsen? Uw computer kan wel uniform verdeelde getallen op het interval \([0,1]\) nabootsen.
Zij \(U\sim \text{Unif}[0, 1]\). Dan heeft \(X = Q(U)\) verdelingsfunctie \(F\).
Hoe ziet een QQ plot uit van de logarithmes van de quantielen van 'e'en verdeling uit deze klasse van verdelingen tegen de logaritmes van de quantielen van een andere verdeling van dezelfde klasse? Noem de parameters van de twee verdelingen \(\theta_1\) en \(\theta_2\) en gebruik de \(x\)-as voor de log quantielen van de eerste verdeling, de \(y\)-as voor de log quantielen van de andere.
Voor een gegeven kans \(p\) zijn de twee quantielen \(x = p^{1/\theta_1}\), \(y = p^{1/\theta_2}\). Dus \(\log(x) = \log(p)/\theta_1\), \(\log(y) = \log(p)/\theta_2\). Als ik de variabele \(p\) uit deze twee vergelijkingen elimineer, vind ik \(\log(y) = (\theta_1/\theta_2)\log(x)\).
Dus de QQ plot van log quantielen van de ene verdeling tegen de andere is een rechte lijn door de oorsprong met helling \(\theta_1/\theta_2\).