# # Stochastik-Praktikum, WS 2020/2021, JGU Mainz # 11.1.2021 ## ################################################## ## heute: ## 1) Reprise zur Brownschen Bewegung ("Pfad-Version des ZGWS") ## ## 2) Zur linearen Regression ## ## ################################################## # #################################################### # # Nochmakl zur Brownschen Bewegung (vgl. auch Sitzung von 4.1.21): # Normalverteilte, "stationäre" (in Verteilung), # unabhängige Zuwächse # Simuliere längs einem Gitter von Zeitwerten mit Abstand delta: delta <- 0.001 # variieren: delta <- 0.002; delta <- 0.001 N <- ceiling(1/delta) zeitgitter <- (0:N)/N # Einen Brownschen Pfad auf die offensichtliche Weise simulieren: BBpfad <- seq(0, along=zeitgitter) for (i in 1:N) { BBpfad[i+1] <- BBpfad[i]+rnorm(1, sd=sqrt(delta)) } # und anschauen plot(zeitgitter, BBpfad, type="l", xlab="Zeit",ylab="", ylim=c(-2,2), main="Pfad der BB") # für ein festes Zeitintervall geht es schneller wie oben mittels BBpfadsim <- function(n, delta) c(0,cumsum(rnorm(n, sd=sqrt(delta)))) BBpfad2 <- BBpfadsim(N, delta) points(zeitgitter, BBpfad2, type="l") # # Betrachten wir wieder die gemeinsame Verteilung der Pfade zur # Zeit 0.5 und zur Zeit 1: wdh <- 1000 beob.BB <- matrix(0, nrow=wdh, ncol=2) for (i in 1:wdh) { zw1 <- rnorm(1, sd=sqrt(0.5)); zw2 <- rnorm(1, sd=sqrt(0.5)) beob.BB[i,]<-c(zw1, zw1+zw2) } # Vert. zur (reskalierten) Zeit 0.5: mean(beob.BB[,1]); var(beob.BB[,1]) qqnorm(beob.BB[,1]) # Vert. zur (reskalierten) Zeit 1: mean(beob.BB[,2]); var(beob.BB[,2]) qqnorm(beob.BB[,2]) # gemeinsam: plot(beob.BB[,1], beob.BB[,2], xlab="Pos. zur Zeit 0.5", ylab="Pos. zur Zeit 1") # Gemeinsame Verteilung: # Welcher Anteil hat Pos. zur Zeit 0.5 <= x1 und # Pos. zur Zeit 1.0 <= x2 ? x1 <- 0.5; x2 <- -0.1 sum(beob.BB[,1] <= x1 & beob.BB[,2] <=x2)/N ####################################################### # # Betrachten wir zur Zeit 1, wie weit der Pfad nun oberhalb # seines bisherigen Minimalwerts liegt # (beachte: dies hat nicht die Eigenschaft, # von nur endlich vielen Zeitpunkten abzuhängen) # Wir betrachten wie oben drei Verteilungen von Inkrementen: N <- 1000 # auch variieren, N <- 10000 delta <- 1/N p <- 0.4 mu <- p*1+(1-p)*(-1) # sigma=sqrt(1-mu^2) pfadsim.gewIrrf <- function(n) c(0,(cumsum(2*rbinom(n, size=1, prob=p)-1)-(1:n)*mu)/(sqrt(1-mu^2)*sqrt(n))) lambda <- 1.8 # ggfs. variieren, z.B. 2.0 # hier: mu= 1/lambda # sigma= 1/lambda pfadsim.ExpIrrf <- function(n) c(0,(cumsum(rexp(n, rate=lambda))-(1:n)*(1/lambda))/((1/lambda)*sqrt(n))) # hier delta=1/n pfadsim.BB <- function(n) c(0,cumsum(rnorm(n, sd=sqrt(1/n)))) # wdh <- 1000 Beob <- matrix(0, nrow=3, ncol=wdh) for (i in 1:wdh) { pf <- pfadsim.gewIrrf(N); Beob[1,i] <- pf[N]-min(pf) pf <- pfadsim.ExpIrrf(N); Beob[2,i] <- pf[N]-min(pf) pf <- pfadsim.BB(N); Beob[3,i] <- pf[N]-min(pf) } # hist(Beob[1,], prob=TRUE, main=c("Reskalierte Gew. Irrf.","empir. Vert von Endpunkt - Minimum")) hist(Beob[2,], prob=TRUE, main=c("Reskalierte Irrf. mit exp-verteilten Inkrementen","empir. Vert von Endpunkt - Minimum")) hist(Beob[3,], prob=TRUE, main=c("Reskalierte Irrf. mit normalverteilten Inkrementen","empir. Vert von Endpunkt - Minimum")) qqplot(Beob[1,],Beob[3,]) qqplot(Beob[2,],Beob[3,]) ## ############################################################################## ## ## [Zur] Lineare[n] Regression ## ## ######################################################### # Lineare Regression, ein Beispiel: # Rs cats-Datensatz library("MASS") data(cats) ?cats attach(cats) koerper <- Bwt[Sex=="M"]; herz <- Hwt[Sex=="M"] # (Wir betrachten nur die 97 Kater) plot(koerper, herz, xlab='Körpergewicht [kg]',ylab='Herzgewicht [g]') # wir passen die Ausgleichsgerade # Herzgewicht = a + b * Koerpergewicht # an: modell <- lm(herz ~ koerper) modell # die Koeffizienten: modell$coefficients # die Residuen: modell$residuals # zeichnen wir die Regressionsgerade ein: abline(modell, col='red') # (aequivalent ist hier:) abline(modell$coefficients[1], modell$coefficients[2]) # hist(modell$residuals) qqnorm(modell$residuals) # Praediktion: a <- modell$coefficients[1]; b <- modell$coefficients[2] neues.kg <- 3.75 predicted.hg <- b * neues.kg + a predicted.hg points(neues.kg, predicted.hg, pch=22, col='blue') # Bem.: Extrapolation kann natuerlich auch in der Anwendungssituation # "sinnlose" Werte ergeben, z.B. neues.kg <- 0.25 # ein sehr leichter Kater predicted.hg <- b * neues.kg + a predicted.hg # (das waere eine negatives vorhergesagtes Herzgewicht, # was natuerlich keinen Sinn hat ...) detach(cats) ## ###################################################### # # Woher kommt der Name "Regression"? # Francis Galtons Beobachtung der "Rueckkehr zum Mittelwert" ## Beispiel: K. Pearsons Daten über Größen von Vätern und Söhnen pearson <- read.table("pearson.dat", header=T) pearson[1:10,] attach(pearson) length(fheight) plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)") # Mittelwerte einzeichnen mean(fheight) # 67.6871 mean(sheight) # 68.68407 abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) mo <- lm(sheight ~ fheight) mo # und die Regeressionsgerade: abline(mo, lwd=2, col="red") # (empirische) Korrelation cor(fheight, sheight) # 0.5013383 sq.f <- sum((fheight-mean(fheight))^2)/length(fheight) sq.f # 7.527313 sqrt(sq.f) # 2.743595 sq.s <- sum((sheight-mean(sheight))^2)/length(sheight) sq.s # 7.915196 sqrt(sq.s) # 2.813396 cov.fs <- sum((fheight-mean(fheight))*(sheight-mean(sheight)))/length(fheight) cov.fs # 3.869739 cov.fs/(sqrt(sq.f)*sqrt(sq.s)) # Korrelationskoeff: 0.5013383 # uebrigens: die Daten (zumindest die Marginalien) sind nicht unplausibel # durch die Normalverteilung beschrieben: hist(fheight) qqnorm(fheight) hist(sheight) qqnorm(sheight) ## ######################################### ## grosse und kleine Vaeter anschauen: ## grosse: gr1 <- 71.5; gr2 <- 72.5 ausw.gr <- (fheight >= gr1) & (fheight <= gr2) sum(ausw.gr) farbe.gr <- "blue" plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="n") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!ausw.gr], sheight[!ausw.gr], pch=20, cex=0.5) points(fheight[ausw.gr], sheight[ausw.gr], pch=20, cex=0.5, col=farbe.gr) points((gr1+gr2)/2, mean(sheight[ausw.gr]), pch=19, cex=1.3, col=farbe.gr) ##points((gr1+gr2)/2, mean(sheight[ausw.gr]), pch=19, cex=0.7, col="white") ## points(fheight, sheight, pch=20, cex=0.5) # Histogramme vergleichen kgr <- 57:80 kgr2 <- seq(from=58, to=80, by=2) hi1 <- hist(sheight, breaks=kgr, plot=F) hi2 <- hist(sheight[ausw.gr], breaks=kgr2, plot=F) plot(hi1$mids, hi1$density, xlim=c(55, 80), ylim=c(0,1.1*max(hi1$density, hi2$density)), xlab="", ylab="Dichte", type="l") abline(h=0, lty=1, col="darkgray") points(hi2$mids, hi2$density, type="l", col=farbe.gr) abline(v=mean(sheight), lty=2) abline(v=mean(sheight[ausw.gr]), lty=2, col=farbe.gr) ## zusaetzlich mit Rand-Histogrammen: ## ## gr-Vaetern werden markiert: x11(width=9) layout(matrix(c(1,2),nrow=1), widths=c(3,1), heights=c(1,1)) # Scatter: par(mar=c(4.1,4.1,4.1,1)) plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="n") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!ausw.gr], sheight[!ausw.gr], pch=20, cex=0.5) points(fheight[ausw.gr], sheight[ausw.gr], pch=20, cex=0.5, col=farbe.gr) points((gr1+gr2)/2, mean(sheight[ausw.gr]), pch=19, cex=1.3, col=farbe.gr) abline(v=gr1, lty=3, col=farbe.gr); abline(v=gr2, lty=3, col=farbe.gr) lines(c((gr1+gr2)/2,82),rep(mean(sheight[ausw.gr]),2),col=farbe.gr,lty=2) # dann Histogramm: par(mar=c(4.1,0,4.1,0.2)) plot(hi1$density, hi1$mids, ylim=c(55, 80), xlim=c(0,1.1*max(hi1$density, hi2$density)), ylab="", xlab="", #xlab="Dichte", main=c("Dichte","(Gr. d. Sohns)"), type="l") abline(v=0, lty=1, col="darkgray") points(hi2$density, hi2$mids, type="l", col=farbe.gr) abline(h=mean(sheight), lty=2) abline(h=mean(sheight[ausw.gr]), lty=2, col=farbe.gr) ## ############################################### ## kleine: kl1 <- 64.5; kl2 <- 65.5 ausw.kl <- (fheight >= kl1) & (fheight <= kl2) sum(ausw.kl) farbe.kl <- "darkgreen" plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="n") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!ausw.kl], sheight[!ausw.kl], pch=20, cex=0.5) points(fheight[ausw.kl], sheight[ausw.kl], pch=20, cex=0.5, col=farbe.kl) points((kl1+kl2)/2, mean(sheight[ausw.kl]), pch=19, cex=1.3, col=farbe.kl) ##points((kl1+kl2)/2, mean(sheight[ausw.kl]), pch=19, cex=0.7, col="white") ## points(fheight, sheight, pch=20, cex=0.5) # Histogramme vergleichen kkl <- 57:80 kkl2 <- seq(from=58, to=80, by=2) hi1 <- hist(sheight, breaks=kkl, plot=F) hi3 <- hist(sheight[ausw.kl], breaks=kkl2, plot=F) plot(hi1$mids, hi1$density, xlim=c(55, 80), ylim=c(0,1.1*max(hi1$density, hi3$density)), xlab="", ylab="Dichte", type="l") abline(h=0, lty=1, col="darkgray") points(hi3$mids, hi3$density, type="l", col=farbe.kl) abline(v=mean(sheight), lty=2) abline(v=mean(sheight[ausw.kl]), lty=2, col=farbe.kl) ## mit Rand-Histogrammen, # kl-Vaetern markiert: x11(width=9) layout(matrix(c(1,2),nrow=1), widths=c(3,1), heights=c(1,1)) # Scatter: par(mar=c(4.1,4.1,4.1,1)) plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="p") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!ausw.kl], sheight[!ausw.kl], pch=20, cex=0.5) points(fheight[ausw.kl], sheight[ausw.kl], pch=20, cex=0.5, col=farbe.kl) points((kl1+kl2)/2, mean(sheight[ausw.kl]), pch=19, cex=1.3, col=farbe.kl) abline(v=kl1, lty=3, col=farbe.kl); abline(v=kl2, lty=3, col=farbe.kl) lines(c((kl1+kl2)/2,82),rep(mean(sheight[ausw.kl]),2),col=farbe.kl,lty=2) # dann Histogramm: par(mar=c(4.1,0,4.1,0.2)) plot(hi1$density, hi1$mids, ylim=c(55, 80), xlim=c(0,1.1*max(hi1$density, hi3$density)), ylab="", xlab="", #xlab="Dichte", main=c("Dichte","(Gr. d. Sohns)"), type="l") abline(v=0, lty=1, col="darkgray") points(hi3$density, hi3$mids, type="l", col=farbe.kl) abline(h=mean(sheight), lty=2) abline(h=mean(sheight[ausw.kl]), lty=2, col=farbe.kl) ## ############################################### ## alles zusammen: kl und gr x11(width=9) layout(matrix(c(1,2),nrow=1), widths=c(3,1), heights=c(1,1)) # Scatter: par(mar=c(4.1,4.1,4.1,1)) plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="n") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!(ausw.gr | ausw.kl)], sheight[!(ausw.gr | ausw.kl)], pch=20, cex=0.5) points(fheight[ausw.gr], sheight[ausw.gr], pch=20, cex=0.5, col=farbe.gr) points((gr1+gr2)/2, mean(sheight[ausw.gr]), pch=19, cex=1.3, col=farbe.gr) abline(v=gr1, lty=3, col=farbe.gr); abline(v=gr2, lty=3, col=farbe.gr) lines(c((gr1+gr2)/2,82),rep(mean(sheight[ausw.gr]),2),col=farbe.gr,lty=2) points(fheight[ausw.kl], sheight[ausw.kl], pch=20, cex=0.5, col=farbe.kl) points((kl1+kl2)/2, mean(sheight[ausw.kl]), pch=19, cex=1.3, col=farbe.kl) abline(v=kl1, lty=3, col=farbe.kl); abline(v=kl2, lty=3, col=farbe.kl) lines(c((kl1+kl2)/2,82),rep(mean(sheight[ausw.kl]),2),col=farbe.kl,lty=2) # dann Histogramme: par(mar=c(4.1,0,4.1,0.2)) plot(hi1$density, hi1$mids, ylim=c(55, 80), xlim=c(0,1.1*max(hi1$density, hi2$density)), ylab="", xlab="", #xlab="Dichte", main=c("Dichte","(Gr. d. Sohns)"), type="l") abline(v=0, lty=1, col="darkgray") points(hi2$density, hi2$mids, type="l", col=farbe.gr) abline(h=mean(sheight), lty=2) abline(h=mean(sheight[ausw.gr]), lty=2, col=farbe.gr) points(hi3$density, hi3$mids, type="l", col=farbe.kl) abline(h=mean(sheight[ausw.kl]), lty=2, col=farbe.kl) ## ################################# ## Andersherum: Vater als Funktion des Sohns mo.andersherum <- lm(fheight ~ sheight) mo.andersherum plot(sheight, fheight, main="1078 Größen von Sohn und Vater", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, ylab="x[Vater] : Größe des Vaters (in Inches)", xlab="x[Sohn] : Größe des Sohns (in Inches)") abline(v=mean(sheight), lty=2); abline(h=mean(fheight), lty=2) abline(mo.andersherum, lwd=2, col="violet") abline(0, 1) # Diagonale einzeichnen # detach(pearson)