## Stochastik-Praktikum, WS 2010/2011 ## Sitzung vom 14.2.2011 ## ## Häufigkeiten bei Lottozahlen ## und multiples Testen ## # Beim Lotto 6 aus 49 ist für eine gegebene Zahl die Wahrscheinlichkeit, # dass sie in einer Ziehung vorkommt, gleich # choose(48,5)/choose(49,6)=(48!/(5!*43!))/(49!/(6!*43!))=(48!/49!)*(6!/5!) # =6/49 theor.wert <- 6/49 theor.wert # Wir verfügen über folgende Beobachtungen: # Für die Mittwoch- und Samstagziehungen des Lottos 6 aus 49 # (ohne Zusatzzahl) von 1955 bis 2011 (letzte erfasste Ziehung 9.2.2011), # insgesamt 4933 Ziehungen wurde für jede Zahl i=1,2,...,49 # gezählt, wie häufig sie bei diesen Ziehungen vorkam. # Diese Häufigkeitstabelle wurde der Website # http://www.dielottozahlende.net/ entnommen, die Veranstalter # des Stochastik-Praktikums stimmen mit den dort # geäußerten Ansichten (über das bloße Zahlenmaterial hinaus) # aber keineswegs überein. gesamt.ziehungen <- 4933 # Wie oft erwarten wir, eine gegebene Zahl in den Ziehungen zu sehen? gesamt.ziehungen*theor.wert # Die Daten: lz <- read.table("Lottozahlen_1955_2011.txt", header=TRUE) lz # ein Trick, eine etwas übersichtlichere Darstellung der # Häufigkeiten zu erzwingen: Erzeuge einen # Häufigkeitsvektor mit Namen hmn <- lz$Haeufigk; names(hmn) <- 1:49 hmn # Wir sehen beispielsweise: Die 13 ist am seltensten gezogen worden. # Die Häufigkeiten grafisch inspizieren: plot(lz$Zahl, lz$Haeufigk, xlab="Zahl", ylab="Häufigkeit", col="red", main="Häufigkeit der Zahlen 1-49 in 4933 Lottoziehungen") for (i in 1:49) lines(c(i,i),c(0,lz$Haeufigk[i])) # Theoret. Erwartungswert einzeichnen: abline(h=gesamt.ziehungen*theor.wert, lty=2) ########################################################### # Frage: Sind die beobachteten Abweichungen vom # theoretischen Erwartungswert signifikant, oder können sie # genausogut durch reine Zufallsfluktuationen erklärt werden? # Sei B[i] die Anzahl Ziehungen, in denen die Zahl i beobachtet wurde. # Das Nullmodell "reiner Zufall" besagt: # Für jede Zahl i ist die Wahrscheinlichkeit p[i], dass i in einer # gewissen Ziehung vorkommt, gleich 6/49, und die Ziehungen sind unabhängig, # d.h. ist B[i] Binomial(4933,6/49)-verteilt. # Erster Vorschlag zur Prüfung des Nullmodels: # Wir konstruieren für jede Zahl i ein Konfidenzintervall für p[i], # sagen wir, zum Irrtumsniveau alpha <- 0.05 # und prüfen, ob es den theoretischen Wert 6/49 enthält # (dies entspricht einem Z-Test, wenn wir ein "Standard"-Konf.int. benutzen): n <- gesamt.ziehungen # Standard-Konfidenzintervall (basierend auf asymptot. Normalität) # für den Erfolgsparameter im Binomialmodell # (b Erfolge unter n Versuchen beobachtet): # phut = b/n, Konf.int = phut +/- qnorm(1-alpha/2)*sqrt(phut*(1-phut))/sqrt(n) phut <- lz$Haeufigk/n KI.untergr <- phut-qnorm(1-alpha/2)*sqrt(phut*(1-phut))/sqrt(n) KI.obergr <- phut+qnorm(1-alpha/2)*sqrt(phut*(1-phut))/sqrt(n) # Ergebnis grafisch inspizieren: plot(lz$Zahl, phut, xlab="Zahl", ylab="relative Häufigkeit", ylim=c(min(KI.untergr),max(KI.obergr)), col="red", main="Relative Häufigkeit der Zahlen 1-49 in 4933 Lottoziehungen") abline(h=theor.wert, lty=2) for (i in 1:49) lines(c(i,i), c(KI.untergr[i],KI.obergr[i]), lwd=2, col="blue") # Wir sehen: Das Konfidenzintervall für p[13] überdeckt 6/49 nicht: KI.untergr[13]; KI.obergr[13]; 6/49 # und tatsächlich ist die Wahrscheinlichkeit, das eine # Bin(4933,6/49)-verteilte Zufallsgröße einen Wert <= 545 annimmt, # ziemlich klein: pbinom(lz$Haeufigk[13], size=gesamt.ziehungen, prob=6/49) # und auch wenn man die entsprechende symmetrische Abweichung nach oben # mit betrachtet, wird es nicht viel besser: pbinom(lz$Haeufigk[13], size=gesamt.ziehungen, prob=6/49)+ pbinom(2*gesamt.ziehungen*theor.wert-lz$Haeufigk[13], size=gesamt.ziehungen, prob=6/49, lower.tail=FALSE) # so gesehen könnten wir die Nullhypothese # "die Anz. Ziehungen, in denen die 13 vorkommt, ist Bin(4933,6/49)-verteilt" # also beispielsweise zum Signifikanzniveau 5% ablehnen. ########################################################### # # Haben wir also das Nullmodell "reiner Zufall" (statistisch) widerlegt? # Sollten wir also Alternativen ins Auge fassen, z.B. # "Betrug beim Lotto" oder "die 13 hat systematisch kleinere # Ziehungswahrscheinlichkeit"? # Wir haben es mit einem Problem des multiplen Testens zu tun: # Eigentlich haben wir 49 Tests (nämlich auf die 49 Nullhypothesen # "p[i]=6/49", i=1,...,49) zugleich ausgeführt, und für jeden # dieser Tests gilt, dass er nur mit W'keit <= 0.05 die # Nullhypothese ablehnt, obwohl sie zutrifft (d.h. dass wir einen sog. # Fehler 1. Art begehen). # Andererseits: Nehmen wir an, die 49 Nullhypothesen sind alle wahr # (und der Einfachheit halber: die Tests wären unabhängig), wie # wahrscheinlich wäre es dann, dass mindestens einer "aus Versehen" # anschlägt? # Zunächst: Antwort per Simulation, ggf. einige Male wiederholen # Erzeuge eine zufällige Pseudo-Ziehung der Lottozahlen von 1955-2011 # (wo wir so tun, als ob für jede Zahl unabhängig entschieden würde, # wie oft sie vorkommt): dummy.beob <- rbinom(49, size=gesamt.ziehungen, prob=6/49) KI.ueberdeckt.nicht <- function(b) { ph <- b/gesamt.ziehungen if (ph-qnorm(1-alpha/2)*sqrt(ph*(1-ph))/sqrt(gesamt.ziehungen) > 6/49) return(TRUE) if (ph+qnorm(1-alpha/2)*sqrt(ph*(1-ph))/sqrt(gesamt.ziehungen) < 6/49) return(TRUE) FALSE } ablehnungen <- sapply(dummy.beob, KI.ueberdeckt.nicht) sum(ablehnungen) # falls obiges >= 1 spaßeshalber: Wer war's, und wie groß? (1:49)[ablehnungen]; dummy.beob[ablehnungen] # (ggf. Simulation wiederholen) # (Eigentlich bräuchten wir für obige Frage nicht zu simulieren: # In dieser Bilderbuchsituation ist die Anzahl Tests, die # "aus Versehen" anschlägt (als die Anz. Fehler 1. Art, die wir # insgesamt begingen), Bin(49,0.05)-verteilt: 1-0.05^0*(1-0.05)^49 # Allerdings: Selbst unter der Nullhypothese "reiner Zufall beim Lotto" # sind die Tests (im Gegensatz zu obiger "Bilderbuchsituation") nicht # unabhängig -- das Auftreten verschiedener Zahlen i != j in derselben # Ziehung ist nämlich (negativ) korreliert. ################################################################### # # Ein sehr allgemeiner Ansatz, bei mehrfachen Tests das # Signifikanzniveau richtig einzustellen, ist die sogenannte # Bonferroni-Korrektur # (nach Carlo Emilio Bonferroni, 1892-1960, ital. Mathematiker): # # Bei gewünschtem Gesamt-Signifikanzniveau alpha ersetzt man # bei k Tests in jedem individuellen Test das Signifikanzniveau durch # alpha/k. # Sei also für i=1,...,k ein Test T_i gegeben, für den die # Wahrscheinlichkeit eines Fehlers 1. Art (<)= alpha/k ist, d.h. # unter der Nullhypothese gilt # P(T_i verwirft) <= alpha/k. # Dann hält der kombinierte Test # T_{komb} : verwerfe die Nullhypothese, wenn mindestens einer der Tests # T_1,...,T_k verwirft # das Signifikanzniveau alpha ein, denn unter der Nullhypothese gilt # P(T_{komb} verwirft) # = P(T_1 verwirft oder T_2 verwirft ... oder T_k verwirft) # <= P(T_1 verwirft) + ... + P(T_k verwirft) <= k*(alpha/k)=alpha. (*) # # Das Argument gilt auch, wenn die Tests T_1,...,T_k nicht # unabhängig sind. # # Ein analoger Ansatz greift für Konfidenzintervalle, dies # verfolgen wir nun: alpha.bonf <- 0.05/49 n <- gesamt.ziehungen # Standard-Konfidenzintervall (mit Bonferroni-korrigiertem Niveau) phut <- lz$Haeufigk/n KI.untergr <- phut-qnorm(1-alpha.bonf/2)*sqrt(phut*(1-phut))/sqrt(n) KI.obergr <- phut+qnorm(1-alpha.bonf/2)*sqrt(phut*(1-phut))/sqrt(n) # Ergebnis grafisch inspizieren: x11() plot(lz$Zahl, phut, xlab="Zahl", ylab="relative Häufigkeit", ylim=c(min(KI.untergr),max(KI.obergr)), col="red", main=c("Relative Häufigkeit der Zahlen 1-49 in 4933 Lottoziehungen","und Konfidenzintervalle mit Bonferroni-Korrektur")) abline(h=theor.wert, lty=2) for (i in 1:49) lines(c(i,i), c(KI.untergr[i],KI.obergr[i]), lwd=2, col="blue") # Wir sehen: Jetzt überdecken alle Konfidenzintervalle den Wert 6/49. # # Andererseits: Die Bonferroni-Korrektur ist zwar korrekt, aber auch # sehr "konservativ". Obige Ungleichung (*) ist i.A. (überhaupt) # nicht scharf. Ein Bonferroni-korrigierter Gesamttest könnte # daher nur eine geringe Macht haben (d.h. eine recht große # Wahrscheinlichkeit für einen Fehler 2. Art). # Ein Skeptiker könnte also entgegnen, dass es systematische # Verzerrungen (vom theoretischen Wert 6/49) der # Erfolgswahrscheinlichkeiten verschiedener Zahlen beim Lotto gibt, # die wir nur wegen der relativ langen Bonferroni-Konfidenzintervalle # nicht mehr entdecken können. # Diesem Einwand können wir folgendermaßen begegnen: # Zwar ist die Verteilung des größten und des kleinsten # beobachteten Werts in den Häufigkeiten der Zahlen 1,...,49 # in 4933 Ziehungen von 6 aus 49 unter der Nullhypothese "reiner # Zufall" schwierig exakt zu bestimmen, aber wir können sehr leicht mit # R (unter unserer Nullhypothese) die # Geschichte der 4933 Ziehungen nachsimulieren. # Wir tun dies oft, sagen wir 100 mal, und zählen mit, wie oft # die seltenste Zahl nur 545 oder seltener auftrat # (und symmetrisch auch, wie oft die häufigste Zahl 650 mal # oder häufiger auftrat). ################################################################### # # Ansatz per Simulationsstudie # eine Lottoziehung sample(1:49, 6, replace=FALSE) # 4933 Lottoziehungen: simHaeufigk <- numeric(49) for (t in 1:gesamt.ziehungen) { ziehung <- sample(1:49, 6, replace=FALSE) simHaeufigk[ziehung]<-simHaeufigk[ziehung]+1 } # Häufigkeiten der Zahlen 1,...,49 in den 4933 Ziehungen: simHaeufigk # kleinster und größter Wert darunter: min(simHaeufigk); max(simHaeufigk) # Simuliere die Geschichte der Lottoziehungen 1955-2011 # wdh mal, speichere jeweils Häufigkeit der # seltensten und der häufigsten Zahl wdh <- 100 minima <- numeric(wdh); maxima <- numeric(wdh) for (w in 1:wdh) { simHaeufigk <- numeric(49) for (t in 1:gesamt.ziehungen) { ziehung <- sample(1:49, 6, replace=FALSE) simHaeufigk[ziehung]<-simHaeufigk[ziehung]+1 } minima[w] <- min(simHaeufigk); maxima[w] <- max(simHaeufigk) } # die beobachteten Häufigkeiten der jeweils seltensten Zahl hist(minima) abline(v=lz$Haeufigk[13], col="red") # und die Häufigkeit der # seltensten Zahl (13) in den Daten sum(minima <= lz$Haeufigk[13])/wdh # die beobachteten Häufigkeiten der jeweils häufigsten Zahl hist(maxima) lz$Haeufigk[43] abline(v=lz$Haeufigk[43], col="red") # und die Häufigkeit der # häufigsten Zahl (43) in den Daten sum(maxima >= lz$Haeufigk[43])/wdh 2*gesamt.ziehungen*theor.wert-lz$Haeufigk[13] # alternativ könnte man # betrachten: analoge Abweichung nach oben, die die 13 in den # Daten nach unten zeigt. sum(maxima >= lz$Haeufigk[43] & minima <= lz$Haeufigk[13])/wdh #################################################### # # Unser Befund ist also: # Die beobachtete Spannweite von Häufigkeiten ist durchaus # verträglich mit der Nullhypothese "reiner Zufall beim Lotto" # (d.h. jeweils unabhängige, rein zufällige Ziehungen von # 6 Zahlen aus 49 Möglichen ohne Zurücklegen), daher ist die # beobachtete Seltenheit der 13 eigentlich statistisch gar nicht # auffällig.