# # Stochastik-Praktikums, WS 2020/2021, JGU Mainz # 25.1.2021 ## ############################################################# ## ## Themen: ## (*) zum t-Test ## (*) Tests fuer Kategorielle Daten, genauer ## zum Chiquadrat-Test ## zu Fishers exaktem Test ## ## ggfs.: etwas zu Pseudozufallszahlen (falls Zeit bleibt ;-) ## ############################################################# # zum t-Test # t-Statistik: t.stat <- function(x, mu0=0) mean(x-mu0)/(sd(x)/sqrt(length(x))) # ein "klassisches" Datenbeispiel: # aus Student (William S. Gosset), The Probable Error of a Mean, Biometrika 6:1–25 (1908) schlaf <- c(0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0) schlaf mean(schlaf) sd(schlaf) sd(schlaf)/sqrt(10) t.test(schlaf) # betrachten wir die empirische Verteilung von T unter # der Nullhypothese n <- 5 x <- rnorm(n, mean=0, sd=1) x t.stat(x) M <- 10000 # wir machen M simulierte Experimente t.emp <- replicate(M, t.stat(rnorm(n, mean=0, sd=1))) hist(t.emp, breaks=c(min(c(t.emp-1,-5)),seq(from=-4,to=4,by=0.2),max(c(t.emp+1,5))), xlim=c(-4,4), prob=T, xlab='Wert der T-Statistik', main=paste('empir. Verteilung von T_n, n=',n,sep='')) curve(dnorm, col='red', add=T) # wir fuegen die Normaldichte hinzu curve(dt(x,df=n-1), col='blue', add=T) # und die Student-(n-1)-Dichte # mit verschiedenen Werten von n wiederholen: n <- 7 n <- 10 n <- 20 n <- 100 # wir koennen die Guetefunktion und die "Robustheit" # des T-Tests per Simulation einschaetzen (s.a. 10. Uebungsblatt) ## ########################################### ## Zum Chi-Quadrat-Test ## ("feste", d.h. einpunktige Nullhypothese) ## G.Mendel, Versuche ueber Pflanzenhybriden, ## Verhandlungen des naturforschenden Vereines in Bruenn, ## Bd. IV fuer 1865, Abhandlungen: 3–47, (1866) phaenotypen <- c('rund-gelb', 'rund-gruen','kantig-gelb','kantig-gruen') erbsen <- c(315, 108, 101, 32) names(erbsen) <- phaenotypen # Gregor Mendels Befunde aus 556 Erbsen-Kreuzungsexperimenten: erbsen sum(erbsen) erbsen/sum(erbsen) # Anteile # gemaess Mendel'schen Regeln erwartete Anteile # (bei Kreuzung von "Doppelhybriden"), # 'rund' und 'gelb' sind dominant, 'kantig' und 'gruen' rezessiv theoret <- c(9/16,3/16,3/16,1/16) chisq.test(erbsen, p=theoret) # falls wir uns nicht auf die Asymptotik verlassen wollen: chisq.test(erbsen, p=theoret, simulate.p.value=TRUE) ################################################# # Aus Benson Rosen, Thomas H. Jerdee, # Influence of sex role stereotypes on personnel decisions, # J. Appl. Psych. 59, 9--14, 1974; siehe Table~1 dort # 48 Teilnehmern eines Management-Kurses wurde je eine (fiktive) # Personalakte vorgelegt, und sie sollten anhand der Aktenlage # entscheiden, ob sie die betreffende Person befoerdern oder die # Akte zunaechst ablegen und weitere Kandidaten begutachten wuerden. # Die Akten waren identisch bis auf die Geschlechtsangabe -- # 24 waren als "weiblich und 24 als "maennlich" gekennzeichnet -- # und wurden rein zufaellig an die Teilnehmer verteilt. jr <- matrix(c(14,10,21,3), nrow=2, dimnames=list(c("bef","abl"),c("w","m"))) jr plot(dhyper(0:24,35,13,24),main="Hypergeom. Gewichte (35 bef, 13 abl, ziehe 24 mal)") ## ## zu Fishers "exaktem Test" FisherTest <- function(x, alternative="zweiseitig") { # berechne p-Wert fuer Fishers exakten Test, # x sollte eine 2x2-Matrix sein, unsere # "Testgroesse" ist Eintrag x[1,1] zs <- c(x[1,1]+x[1,2],x[2,1]+x[2,2]) # Zeilensummen ss <- c(x[1,1]+x[2,1],x[1,2]+x[2,2]) # Spaltensummen ges <- sum(x) if (alternative=="zweiseitig") { beob.gew <- dhyper(x[1,1],zs[1], zs[2], ss[1]) alle.gew <- dhyper(0:(min(ss[1],zs[1])),zs[1], zs[2], ss[1]) pwert <- sum(alle.gew[alle.gew <= beob.gew]) #if (x[1,1] >= ss[1]*zs[1]/ges) { # pwert<-(phyper(x[1,1]-1, zs[1], zs[2], ss[1], lower.tail=FALSE) # +phyper(zs[1]-x[1,1],zs[1], zs[2], ss[1])) #} else { # pwert<-(phyper(x[1,1], zs[1], zs[2], ss[1]) # + phyper(zs[1]-x[1,1]-1, zs[1], zs[2], ss[1], lower.tail=FALSE)) # } else if (alternative=="groesser") { pwert <- phyper(x[1,1]-1, zs[1], zs[2], ss[1], lower.tail=FALSE) } else if (alternative=="kleiner") { pwert <- phyper(x[1,1], zs[1], zs[2], ss[1], lower.tail=TRUE) } pwert } FisherTest(jr) ## p.s.: "eingebauter" R-Befehl: fisher.test ## ######################################################### ## Peter B. Hjortrup, Nicolai Haase, Jorn Wetterslev, Anders Perner, ## Gone fishing in a fluid trial. ## Critical Care and Resuscitation Volume 18 Issue 1 (Mar 2016) fisch <- matrix(c(45,25,380,348),nrow=2,ncol=2, dimnames=list(c("ueberl","gest"),c("Fische","anderes"))) fisch sum(fisch) uews.hut.ges <- (45+380)/798 uews.hut.ges uews.hut.fische <- 45/(25+45) uews.hut.anderes <- 380/(380+348) uews.hut.fische; uews.hut.anderes uews.hut.fische/uews.hut.anderes FisherTest(fisch, alternative="groesser") # uebrigens: wenn man zweiseitig testet, sieht es # schon nicht mehr signifikant aus: FisherTest(fisch) hjortrup <- read.table("sepsis.txt", header=TRUE) hjortrup attach(hjortrup) sum(total) sum(died) # schauen wir die anderen Sternzeichen ebenfalls an anz.ges <- sum(total) anz.gest <- sum(died) anz.ueberl <- anz.ges-sum(died) pwerte <- matrix(0, nrow=12, ncol=3) for (i in 1:12) { gest <- hjortrup[i,3] ueberl <- hjortrup[i,2]-hjortrup[i,3] mat <- matrix(c(ueberl,gest,anz.ueberl-ueberl,anz.gest-gest),nrow=2) pwerte[i,1] <- FisherTest(mat) pwerte[i,2] <- FisherTest(mat,alternative="groesser") pwerte[i,3] <- FisherTest(mat,alternative="kleiner") } hjortrup round(pwerte, digits=3) # Generiere "zufaellige Permutation": neudaten <- function() { ueberl <- numeric(12) lebende <- anz.ueberl tote <- anz.gest for (i in 1:12) { ueberl[i] <- rhyper(1,lebende,tote,hjortrup$total[i]) lebende <- lebende-ueberl[i] tote <- tote-(hjortrup$total[i]-ueberl[i]) } ueberl } nd <- neudaten() nd hjortrup.neu <- hjortrup; hjortrup.neu[,3] <- nd hjortrup.neu ## wie saehen nun die p-Werte aus? pwerte.neu <- matrix(0, nrow=12, ncol=3) for (i in 1:12) { gest <- hjortrup.neu[i,3] ueberl <- hjortrup.neu[i,2]-hjortrup.neu[i,3] mat <- matrix(c(ueberl,gest,anz.ueberl-ueberl,anz.gest-gest),nrow=2) pwerte.neu[i,1] <- FisherTest(mat) pwerte.neu[i,2] <- FisherTest(mat,alternative="groesser") pwerte.neu[i,3] <- FisherTest(mat,alternative="kleiner") } hjortrup.neu round(pwerte.neu, digits=3) min(pwerte.neu) arrayInd(which.min(pwerte.neu),dim(pwerte.neu)) metaexperiment <- function() { nd <- neudaten() hjortrup.neu <- hjortrup; hjortrup.neu[,3] <- nd pwerte.neu <- matrix(0, nrow=12, ncol=3) for (i in 1:12) { gest <- hjortrup.neu[i,3] ueberl <- hjortrup.neu[i,2]-hjortrup.neu[i,3] mat <- matrix(c(ueberl,gest,anz.ueberl-ueberl,anz.gest-gest),nrow=2) pwerte.neu[i,1] <- FisherTest(mat) pwerte.neu[i,2] <- FisherTest(mat,alternative="groesser") pwerte.neu[i,3] <- FisherTest(mat,alternative="kleiner") } min(pwerte.neu) } metaexperiment() metaexperiment2 <- function() { nd <- neudaten() hjortrup.neu <- hjortrup; hjortrup.neu[,3] <- nd pwerte.neu <- matrix(0, nrow=12, ncol=1) for (i in 1:12) { gest <- hjortrup.neu[i,3] ueberl <- hjortrup.neu[i,2]-hjortrup.neu[i,3] mat <- matrix(c(ueberl,gest,anz.ueberl-ueberl,anz.gest-gest),nrow=2) pwerte.neu[i,1] <- FisherTest(mat) } min(pwerte.neu) } metaexperiment2() detach(hjortrup) ## ##################################################### ## ## Simulation von Pseudo-Zufallszahlen via lineare Kongruenzen-Generatoren # Fuer sehr viele Simulationsprobleme benoetigt man # sog. "Pseudozufallszahlen": # Beobachtungen oder Werte x[1], x[2], x[3], ... die keine # erkennbare Regelmäßigkeit haben und fuer die Zwecke der # Berechnung als unabhaengig und zufaellig generiert angenommen werden # (duerfen), d.h. wir tun beispielsweise so, als ob die x[1], x[2], ... # durch wiederholtes, unabhaengiges Drehen eines Gluecksrad gewonnen # worden waeren -- obwohl sie aus Praktikabilitaetsgruenden durch einen # (deterministischen) Algorithmus (einen sog. (Pseudo-)Zufallsgenerator) # im Computer erzeugt worden sind. # # Bemerkung: # Eine "Standardreferenz" zu Pseudozufallszahlen ist Chapter 3 in # Donald E. Knuth, The art of computer programming, # Vol. 2 / Seminumerical algorithms, 3rd Ed., Addison-Wesley, 1998. # Das dort behandelte Material geht weit über den hier vorgestellten # Stoff hinaus. # Ein einfaches Beispiel sind die # linearen Kongruenzengeneratoren: # Beginne mit einem (ganzzahligen) Startwert ("Zufallssame", "random seed"). # Wenn der aktuelle Wert x ist, so ist der naechste Wert # (a*x+c) mod M # z.B. M <- 2048; a <- 65; c <- 1 # "Selbstbau"-linearer Kongruenzengenerator: linKonGen <- function() { x.neu <- (a*x + c) %% M # folgendes setzt den Wert der (globalen) Variable x auf den von x.neu # (alternativ: x <<- x.neu) assign("x", x.neu, .GlobalEnv) # gebe "uniform" auf (0,1] verteilten Wert aus: (x.neu+1)/M } # Startwert: x <- 5 # betrachte N Aufrufe: N <- 5000 werte <- replicate(N, linKonGen()) # Anschauen: plot(werte, pch=20, cex=0.5) hist(werte, prob=TRUE) werte <- replicate(500, linKonGen()) # Zeichne Wert gegen Nachfolger-Wert: plot(werte, werte[c(2:500,1)]) # zum Vergleich: x11() werteR <- runif(500) plot(werteR, werteR[c(2:500,1)]) # Perioden können ein Problem sein (speziell bei kleinem M): M <- 16 a <- 5; c<-1 #a <- 4; c <-0 x <- 5 replicate(30, linKonGen()) # Zwei Beispiele: M <- 256 a <- 17; c <- 1 # gibt relativ gleichmaessige Verteilung N <- M werte <- numeric(N) for (i in 1:N) werte[i]<-linKonGen() plot(M*werte, M*werte[c(2:N,1)]) x11() M <- 256 a <- 129; c <- 1 # gibt wenig gleichmaessige Verteilung N <- M werte <- numeric(N) for (i in 1:N) werte[i]<-linKonGen() plot(M*werte, M*werte[c(2:N,1)]) # Rs Zufallsgenerator-Hilfe ?set.seed ############################################ # # Eine Illustration zu Knuths Warnung # "... random numbers should not be generated with a method chosen # at random. Some theory should be used." # (loc. cit., S. 6, siehe auch Exercise 11, S.8) m <- 100000 f <- sample(1:m, m, replace=T) # Definiere "Zufallsfolge" durch y[n+1]=f[y[n]] # mit zufaelligem Startwert y[0]. # Diese wird offenbar schliesslich zyklisch, wie lang ist der # resultierende Zyklus? gesehen <- logical(m) # initialisiert mit m-mal FALSE y <- sample(1:m,1) y # der Startwert anfangslaenge <- 0 while(gesehen[y]==FALSE) { anfangslaenge <- anfangslaenge+1 z<-y gesehen[z]<-TRUE y<-f[y] } anfangslaenge # Wie lang war das Stück vor dem Zyklus? # Wir sind auf einen Zyklus gestossen, die Werte von z und y sind # nun so, dass f[z]=y, d.h. wir haben einen Zyklus # y_0=y, y_1=f[y_0], y_2=f[y_1], ..., y_{n-1}=f[y_{n-2}]=z, y_n=f[z]=y_0. # Bestimme die Laenge des Zyklus: zyklaenge <- 1 yy <- f[y] while (yy != y) { zyklaenge <- zyklaenge+1 yy <- f[yy] # cat(yy, " ") # ggf. ent-kommentieren, um Zykel anzuschauen } zyklaenge # Ggf. einige Male wiederholen. # Beobachtung: auch bei grossem m laeuft die Folge typischerweise # auf recht kurze Zyklen. ############################################ # # (empirische) Tests mit Zufallsgeneratoren # hier ggf. eigenen Zufallsgenerator einsetzen zufgen <- function() runif(1) # zufgen <- linKonGen # Beispielwerte: # M <- 2048; a <- 65; c <- 1 # (das Mini-Beispiel vom Anfang, versagt kläglich) # M <- 10^10; a<-3141592621; c <- 1 # M <- 2^31; a <- 65539; c <- 0 # Dies ist "RANDU", s.a. Knuth, S. 107 # M <- 2^35; a <- 2^18+1; c <- 1 # "Generator F" aus Knuth, S. 47 # ggf. auch Startwert variieren: # x <- 5; x <- 314159265 # 1. Sind die Werte gleichmäßig verteilt? # Wir benutzen den chi^2-Test: # Zerlege (0,1] in kl Klassen (0,1/kl], (1/kl,2/kl],...,((kl-1)/k1,1] # wenn wir n Werte generieren lassen, erwarten wir in jeder Klasse n/kl # Werte (mit "Zufallsfluktuationen") kl <- 5 # variieren n <- 100000 beob <- rep(0, times=kl) for (i in 1:n) { w <- ceiling(kl*zufgen()) # verwandelt Wert aus (0,1] in Wert aus {1,...,kl} beob[w] <- beob[w]+1 } beob beob/n beob/n-rep(1/kl, times=kl) # Einschaetzung der Goesse der Abweichung vom "theoretischen Mittelwert" # via chi^2-Statistik chiquadratwert <- sum((beob-n/kl)^2/(n/kl)) chiquadratwert # Wie wahrscheinlich waere eine solche Abweichung für # "echte" Zufallszahlen (p-Wert des Tests)? pchisq(chiquadratwert, df=kl-1, lower.tail=F) # Dasselbe mit einem R-Befehl: chisq.test(beob) # # Analog für d-Tupel von Werten: # (sollten bei Vergroeberung uniform verteilt sein auf kl^d Klassen) d <- 3 # Andere Werte einsetzen, z.B. 2, 3, 5 kl <- 6 # Andere Werte einsetzen, z.B. 2, 3, 4, 5, 6 beob <- rep(0, times=kl^d) n <- 10000 for (i in 1:n) { # gewinne eine ganze Zahl aus {0,1,...,kl^d-1} # durch Darstellung mit d Ziffern im kl-System stelle <- 1; w <- 0 for (j in 1:d) { w <- w+(ceiling(kl*zufgen())-1)*stelle stelle <- stelle*kl } w <- w+1 # (Rs array-Indizes beginnen bei 1, nicht bei 0) beob[w]<-beob[w]+1 } beob beob-n/(kl^d) beob/n-1/(kl^d) # Wie gut passt die empirische Verteilung von d-Tupeln? chisq.test(beob) # # "Lueckentest": Wie lange muss man warten # (genauer: wieviele Fehlversuche), bis der # naechste Wert in [s,t] kommt? s <- 0.2; t <- 0.5 sim.wartezeit <- function() { i <- 0 repeat { u <- zufgen() if (u >=s && u <=t) break i<-i+1 } i } L <- 50000 wz <- replicate(L, sim.wartezeit()) # vergleiche solche Wartezeiten mit der geometrischen Verteilung # mit Erfolgsparameter t-s: # grafisch: if (max(wz)>=10) { klassen <- c(seq(from=-0.5, to=10.5,by=1), max(wz)+0.5) geom.gewichte <- c(dgeom(0:10, prob=t-s), pgeom(10, prob=t-s, lower.tail=F)) } else { klassen <- seq(from=-0.5, to=10.5,by=1) geom.gewichte <- dgeom(0:10, prob=t-s) } hi <- hist(wz, prob=T, breaks=klassen) points(hi$mids, geom.gewichte/(hi$breaks[-1]-hi$breaks[-length(hi$breaks)]), col="red") # und mit dem chi^2-Test: chisq.test(hi$counts, p=geom.gewichte) # # Maximalwert von t Beobachtungen: # Für U_1,..., U_t u.a. unif([0,1]) ist # P(max(U1,...,U_t)<=x)=x^t = int_0^t t*x^{t-1} dt, # d.h. max(U1,...,U_t) sollte Beta(3,1)-verteilt sein: t <- 3 # ggf. variieren max(replicate(t, zufgen())) L <- 10000 mwerte <- replicate(L, max(replicate(t, zufgen()))) # Vergleiche empirische und theoretische Verteilungsfkt. plot(ecdf(mwerte)) curve(pbeta(x, shape1=3, shape2=1), add=T, col="red") # # Befunde: # Rs eingebauter Zufallsgenerator "besteht" diese Tests, # die Beispiel-LKG haben z.T. Schwierigkeiten. ## ######################################################### ## ######################################################### ## ## Übrigens: den Zustand von Rs eingebauten Zufallsgenerator ## erfährt man mit .Random.seed ## in Aktion: runif(1); .Random.seed[1:6]; runif(1); .Random.seed[1:6] # Informationen zum Zufallsgenerator: RNGkind() # Startwert ("seed") setzen: set.seed(5)