# Statistik-1 ## Ü1. Statistische Grundbegriffe ```r ########################## #Sonderblatt R ########################## #Machen Sie sich mit grundlegenden Funktionen in R vertraut, #indem Sie die folgenden Befehle nacheinander ausführen. x <- c(3,4,3,5,6,3,5,6,7,8) y <- c(1,1,1,1,1,2,2,2,2,2) z <- c(5,2,3,6,1,3,4,5,4,1) x y z length(x) sum(x) cumsum(y) plot(x,z) x[7] x[y==2] m <- matrix(c(x,y,z), nrow=10) m2 <- matrix(c(x,y,z), nrow=10, byrow=TRUE) m3 <- cbind(x,y,z) d <- data.frame(x,y,z) m d str(m) str(d) t(m) #transpose m[7,2] m[2,7] t(m)[2,7] m[2,] d[y==2,] summary(d) plot(d) #Hilfestellung in R ############################### ?mean apropos(mean) #what? help.start() ``` ## **Ü2.** Datenerhebung, Häufigkeiten ![Übung 2 Aufgabe 4](/files/-LMP1NPUXQ1M0e_hGVvO) ```r ### Aufgabe 4 data<-c("Erde","Winter","Sonne","Winter","Winter","Winter","Sonne", "Erde","Frühling","Frühling","Sonne","Norden") stat <- table(data) stat fre<-prop.table(stat) fre absolut <- c(2,4,3,2,1) names(absolut) <- c("Erde","Winter","Sonne","Frühling","Norden") absolut anzahl <- sum(absolut) anzahl relativ <- absolut / anzahl relativ round(relativ, 2) prozent <- relativ * 100 round(prozent, 1) #prepare for the graphic par(mfrow=c(1,1), lwd=1.5, font.axis=2, bty="n", ps=11, cex.axis=1) pie(absolut) barplot(absolut, names.arg=c("E","W","S","F","N"),col="blue", density=c(5,10,15,20,25), ylim=c(0,6)) ``` ![Übung 2 Aufgabe 5](/files/-LMS0rrVaNX6G2AM1H6x) ```r ### Aufgabe 5 absolut <- c(8,20,30,40,22) names(absolut) <- c("ohne","Haupt","Mittlere","Abitur","Uni") absolut anzahl <- sum(absolut) anzahl relativ <- absolut / anzahl round(relativ, 2) prozent <- relativ * 100 round(prozent, 1) stufen <- c(rep(1, 8), rep(2, 20), rep(3, 30), rep(4, 40), rep(5, 22)) stufen st <- table(stufen) st fre<-prop.table(st) fre par(mfrow=c(1,1)) plot.ecdf(stufen) ecdfx<-ecdf(stufen) ecdfx plot.stepfun(ecdfx,lty=1,verticals=F,xlim=c(0,10),ylim=c(0,1),main="",ylab="F (x)",xlab="x") plot.stepfun(ecdfx,lty=2,verticals=T,add=T,main="") title(main="Empirische Verteilungsfunktion") ### Alternativ hi<-c(8,20,30,40,22) names(hi)<-c("ohne", "Hauptschule", "Mittlere Reife", "Abitur", "Universität") # Berechnung der Häufigkeitstabelle fi<-hi / sum(hi) # relative Häufigkeiten Fi<-cumsum(fi) # kumulierte Häufigkeiten HT<-rbind(hi, fi, Fi);HT # Darstellung Häufigkeitstabelle - #Rundung auf 3 Stellen und Vertauschung von Zeilen und Spalten t(round(HT, 3)) # Säulendiagramm (angepasste Abstände wegen der langen Bezeichner) barplot(hi, ylim=c(0,40), space=0.15, cex.names=0.9,col="blue", names.arg=c("ohne","Hauptschule","Mittlere\nReife","Abitur","Universität") ) title(main="Absolute Häufigkeitsverteilung" ``` ![Übung 2 Aufgabe 6](/files/-LMS45CnpOq00JqMDKHi) ```r #Aufgabe 6 v<-c(39, 41, 42, 37, 42, 39, 40, 40, 40, 40, 43, 44, 39, 39, 43, 42, 42, 38, 42, 46) names(v)<-v vsort <- sort(v) n <- length(vsort) table(vsort) table(vsort)/n quantile(vsort,type=2) quantile(v,type=2) #sortiert selber par(mfrow=c(1,1)) boxplot(vsort, col="blue") ``` ![Übung 2 Aufgabe 7](/files/-LMS6Cud6zOUvBGOAfU9) ```r #Aufgabe 7: Pareto-Diagramm defect <- c(18, 2, 7, 44, 24, 9) names(defect) <- c("aa", "bb", "cc", "dd", "ee", "ff") par(mfrow=c(1,1), lwd=2.5, font.axis=2, bty="u",ps=10, cex.axis=0.8, cex=1.3) install.packages("qcc") library(qcc) # spez. Funktion pareto.chart() inlibrary(qcc) pareto.chart(defect, ylab = "Fehler-Häufigkeit", main =" ", xlab="Fehler-Ursache", las=3, col="blue") #alternativ tabelle<-as.table(defect) barplot(tabelle) barplot(sort(tabelle, decreasing=TRUE)) #sortiert nach Häufigkeit barplot(sort(prop.table(tabelle), decreasing=TRUE)) cumsum(sort(tabelle, decreasing=TRUE))/sum(defect) ``` ![Übung 2 Aufgabe 8](/files/-LMS8hG3SgEjcUeUQdBc) ```r ##Aufgabe 8 absolut <- matrix(c(22, 242, 176, 294, 92, 14), nrow=2, byrow=T) colnames(absolut) <- c("weniger 3","3 bis 5", "mehr als 5") rownames(absolut) <- c("Jugendliche", "Erwachsene") names(dimnames(absolut)) <- c("Personenkreis","Spiele"); absolut gesamt<-sum(absolut) tab<-as.table(absolut) #Mosaikplot plot(tab) addmargins(tab) prop.table(tab) #absoluten Häufigkeiten für Personenkreis margin.table(absolut, 1) #absoluten Häufigkeiten für Spiele margin.table(absolut, 2) #relativen Häufigkeiten für Personenkreis = Zeilenprozente round(prop.table(absolut, 1), 3) ztab<-prop.table(tab,1) addmargins(ztab) #relativen Häufigkeiten für Spiele = Spaltenprozente round(prop.table(absolut, 2), 3) stab<-prop.table(tab,2) addmargins(stab) ``` ## Ü3. Verteilungsanalyse, Lagemaße ![Übung 3 Aufgabe 1](/files/-LMSSsDOj1aBWqXiSrS0) ```r #Aufgabe 1 u<-c(3, 7, 12, 18, 19, 20, 25, 25, 27, 28, 29, 31, 32, 34, 37, 38, 40, 41, 45, 47) u h<-hist(u,breaks=(0:5)*10, col="red") h par(mfrow=c(1,1), lwd=1.5, font.axis=2, bty="n", ps=11, cex.axis=1) x <- h$breaks y <- c(0, round(cumsum(h$counts)/length(u),3)) plot(x, y, type="b", ylim=c(0,1), xlim=c(0,50), col="blue", ylab="rel. Summenhäufigkeit" ) abline(v=x, lty=2, col="green") abline(h=y, lty=2, col="red") ``` ![Übung 3 Aufgabe 2](/files/-LMSZ9XWSD-6bqIh1vE-) ```r #Aufgabe 2 b<-c(190, 210, 195, 209, 199, 189, 215) sort(b) sum(b)/length(b) median(b) mean(b) bb<-c(b,149) median(bb) mean(bb) ``` ![Übung 3 Aufgabe 5](/files/-LMSbWZLF9VwA2IJrULd) ```r #Aufgabe 5: harmonischer Mittelwert si<-c(25,15,2) vi<-c(17,12,21) ti <- (1/vi)*si v<-sum(si) / sum(ti) v # mittlere Geschwindigkeit t<-sum(ti) t #Gesamtzeit ``` ![Übung 3 Aufgabe 6](/files/-LMSc_s4S4imsVxg0Fs1) ```r #Aufgabe 6: geometrischer Mittelwert r<-c(7.1 , 6.5 , 6.2 , 5.1 , 4.8 , 5.2) #Renditen x<-1+0.01*r #Wachstumsfaktoren R<-log(x) #Log-Renditen (prod(x))^(1/length(x)) #1.05813 mean(x) #1.05817 mean(r) #5.82 #Zusammenhang mean(R) #5.6 log((prod(x))^(1/length(x))) #5.6 ``` ## Ü4. Histogramm, Streuungsmaße ![Übung 4 Aufgabe 1](/files/-LMSe7tJnaEF2ifbuUB0) ```r #Aufgabe 1 stufen <- c(rep(1, 25), rep(3, 40), rep(6, 10), rep(10, 20), rep(16,5)) windows() par(mfrow=c(1,1), lwd=1.5, font.axis=2, bty="n", ps=11, cex.axis=1) hist(stufen, breaks=c(0, 2,4,8,12,20), col="blue", xlim=c(0,20), xlab="Histogramm der verlorenen Pfunde", main=" ") h<-hist(stufen, breaks=c(0, 2,4,8,12,20)) x <- h$breaks y <- c(0, round(cumsum(h$counts)/length(stufen),3)) plot(x, y, type="b", ylim=c(0,1), xlim=c(0,20), col="blue", xlab="Diät", ylab="rel. Summenhäufigkeit" ) abline(v=9, lty=2, col="green") abline(v=6, lty=2, col="green") abline(v=2, lty=2, col="green") abline(h=.8, lty=2, col="red") abline(h=.25, lty=2, col="red") abline(h=.7, lty=2, col="red") ``` ![Übung 4 Aufgabe 3](/files/-LMSekghMXAAu-mkeMwu) ```r # Aufgabe 3 #Stamm-Blatt-Darstellung time <- c(93, 87, 96, 77, 73, 91, 82, 71, 98, 74, 95, 89, 79, 88) stem(time,scale = 1) stem(time,scale = 2) ``` ![Übung 4 Aufgabe 4](/files/-LMSgNqSsAUzmYLEpoB7) ```r # Aufgabe 4 kiosk<-c(56, 2, 7, 0, 42, 118, 35, 29, 10, 21, 50, 92, 28, 14, 11, 0, 6, 25, 17, 64) windows() par(lwd=1.5, font.axis=2, ps=14, cex.axis=1, bty="o") op <- par(mar=c(0,0,0,0), oma=c(0,0,0,0)+.1) layout(matrix(c(1,1,1,2), nc=1)) br<-quantile(kiosk, c(0,0.25, 0.50, 0.75,1),type=2) hist(kiosk, breaks=br, col="blue", main=" ") boxplot(kiosk, horizontal = TRUE, col = "blue", lwd=2, cex=1.5, pch=16, axes = FALSE) ``` ![Übung 4 Aufgabe 7](/files/-LMShD7Nae4CQB2NGbne) ```r #Aufgabe 7 Partei_A<-c(5.6,6.3,6.6,6.9,7.1,7.6,6.1) Partei_B<-c(40.4,41.9,47.9,40.4,48.9,41.4,42.9) sqrt(sum((Partei_A - mean(Partei_A))^2)/(length(Partei_A)-1)) sqrt(sum((Partei_A - mean(Partei_A))^2)/(length(Partei_A)-1))/mean(Partei_A) sqrt(sum((Partei_B - mean(Partei_B))^2)/(length(Partei_B)-1)) sqrt(sum((Partei_B - mean(Partei_B))^2)/(length(Partei_B)-1))/mean(Partei_B) ``` ## Ü5. Bivariate Datenanalyse ![Übung 5 Aufgabe 1](/files/-LMSu6u94nXQkRVPv6QS) ```r #Aufgabe 1 b<-c(6, 8.5, 0, 11.5, 13, 20, 0, 7.5, 8, 15.5) muc<-c(8, 9.5, 12.5, 0, 9.5, 13, 17, 21, 19, 14.5, 14, 18) bQ<-quantile(b, c(0,0.25, 0.50, 0.75,1),type=2) mucQ<-quantile(muc, c(0,0.25, 0.50, 0.75,1),type=2) par( mfrow=c(1,1) ) plot(bQ, mucQ, xlab="Berlin", ylab="München") abline(0,1, col="red", lty=2) sqrt(sum((b - mean(b))^2)/(length(b))) sqrt(sum((muc - mean(muc))^2)/(length(muc))) 100*sqrt(sum((b - mean(b))^2)/(length(b)))/mean(b) 100*sqrt(sum((muc - mean(muc))^2)/(length(muc)))/mean(muc) ``` ![](/files/-LMckUQcXZFUnSxr9mGL) ```r #Aufgabe2 gini <- function(x, y) { # Funktion zum Gini-Index area <- 0 # Trapezregel for (i in 2:n+1) area <- area + 0.5*((x[i]-x[i-1])*(y[i]+y[i-1])) gini <- 1 - 2*area; round(gini, 3) # Gini-Index } x <- c(120,120,120,120,120,300,300,300,300,1200); n <- length(x) u <- c(0, (1:n)/n) # Abszisse - relativer Index v <- c(0, (cumsum(x) / sum(x))) # Odinate - kumulierte rel. Anteile gini(u, v) u1 <- c(0,0.5,0.9,1) v1 <- c(0,0.2,0.6,1) par(mfrow=c(1,1), lwd=2, font.axis=2, bty="n", ps=12) plot(u1, v1, type="b", cex=1.5, xlim=c(0,1), ylim=c(0,1), xlab="u", ylab="v") abline(0,1, col="red", lty=2) polygon(c(u1,0), c(v1,0), angle = -45, border=NA, density = 10) text(0.85, 0.30, paste("Gini-Index=",round(gini(u, v),3))) ``` ![](/files/-LMckx5yNEFIl4t2ii_J) ![](/files/-LMclDbKtQwPlzsfaeSU) ```r #Aufgabe 3 x_f<-c(1,2,2,3,3) y_f<-c(4,6,3,5,3) set1<-data.frame(S=x_f,E=y_f) x_m<-c(4,5,5,6,6) y_m<-c(6,8,4,7,5) set2<-data.frame(S=x_m,E=y_m) cor(set1) #-0.046 cor(set2) #0 x<-c(x_f,x_m) y<-c(y_f,y_m) set<-data.frame(S=x,E=y) cor(set) #0.503 ``` ![](/files/-LMclzEi9nCodNbd4zdX) ```r #Aufgabe 4 g<-c(9,1,10,6,5,8) w<-c(7,5,12,10,8,2) vin<-data.frame(Gutachter=g,Weinfreund=w) cor(vin, method="spearman") #0.37 cor(vin) g<-c(9,1,10,6,5,9) w<-c(7,5,12,10,8,5) vin<-data.frame(Gutachter=g,Weinfreund=w) cor(vin, method="spearman") #0.43 ``` ![](/files/-LMcmLuI23tiQD9rh1AB) ```r #Aufgabe 5 # Eingabe der Daten: x = Alter, y = Verkaufswert x<-c(0,1,3,5,6) y<-c(65,45,40,30,20) Auto<-data.frame(Alter=x,Verkaufswert=y) plot(Auto,lwd=8,main="Zusammenhang zwischen Alter und Verkaufswert",col="red", ylim=c(0,70)) cor(x, y, method="spearman") cor(x, y, method="pearson") # oder einfacher: cor(x, y) # oder direkt mit cor(Auto, method="spearman") cor(Auto) help(lm) # R Hilfe zur linearen Regression fm<-lm( Verkaufswert ~Alter , data = Auto) coef(fm) # Koeffizienten der Ausgleichsgeraden # Graphische Darstellung: Scatterplot und Ausgleichsgerade plot(Auto,lwd=8,main="Zusammenhang zwischen Alter und Verkaufswert",col="red", ylim=c(0,70)) abline(fm,lwd=2, col="blue") ``` ## Ü6. Zeitreihen, Parameterschätzer ![](/files/-LMcuXl_ryHjulLGWGEX) ```r ################# #Aufgabe 1 #################### t<-seq(from=1,to=8,by=1) y<-c(49,47,56,58,59,53,57,54) # als time series d.y <- ts(y, start=c(2005,1),frequency=1) #filter y_f2 <- filter(d.y, filter=c(1,1,1)/3) y_f1 <- filter(d.y, filter=c(0.5,1,1,1,0.5)/4) par(mfrow = c(2,1) ) plot(d.y, type = "b", main = "Zeitreihe", lwd = 2,col=2) plot(y_f2, type = "b", main = "3 -glatt", lwd = 2,col=2) plot(d.y, type = "b", main = "Zeitreihe", lwd = 2,col=2) plot(y_f1, type = "b", main = "4-glatt", lwd = 2,col=2) #oder y1<-(1/3)*(d.y + lag(d.y,-1)+ lag(d.y,+1)) #zentriert #y2<-0.25*(d.y + lag(d.y,-1)+ lag(d.y,-2)+ lag(d.y,-3)) y2<-0.25*(d.y + lag(d.y,-1)+ 0.5*lag(d.y,-2)+ lag(d.y,+1)+ 0.5* lag(d.y,+2)) #zentriert par(mfrow = c(1,1) ) plot(d.y, ylab="ZR", xlab="Zeit",main="ZR und Glättung", type = "b",col="blue", frame.plot = TRUE, xaxs = "i", pch=4) lines(y1, col = 2, lwd = 2, type = "b") lines(y2, col = 3, lwd = 2, type = "b") ``` ![](/files/-LMcvEPry8BcQoWWETSN) ![](/files/-LMcvHtF5QAaU47JgREj) ![](/files/-LMcvQEYVrZdWx7qV2nJ) ![](/files/-LMcvkO9Q9ahCRbvAnwj) ```r t<-seq(from=1,to=13,by=1) U<-c(11300,10631,9949,9814,9454,8758,8549,7792,7772,7503,6977,6842,6606) # als time series t.U <- ts(U, start=c(1991,1),frequency=1) plot(t.U) #a) #Geometrisches Mittel (6606/11300)^(1/length(U)) # 0.96 -> jährlich 4 Prozent weniger #Prognose für 2004 6606*(6606/11300)^(1/length(U)) #6339 #Prognose für 2005 6606*(6606/11300)^(2/length(U)) #6083 #b) Trendmodell linearisieren #ln U = ln b0+ lnb1 *t u<-log(U) Trend<- lm(log(U) ~ t); a0 <- Trend$coeff[1]; a1<- Trend$coeff[2] coef(Trend) summary(Trend) #en detail sum((t-mean(t))^2) sum((u-mean(u))^2) sum((t-mean(t))*(u-mean(u))) beta<-sum((t-mean(t))*(u-mean(u)))/sum((t-mean(t))^2) alpha<-mean(u)-beta*mean(t) b0<-exp(a0) b1<-exp(a1) # 0.956 #Prognose für 2004: U(14) b0*b1^(14) #6203 #Prognose für 2005: U(15) b0*b1^(14) #5931 ``` ![](/files/-LMcwFS3aYZwCwDDFa3u) ![](/files/-LMcxJ1Pp28dOB-tvEzU) ```r #Aufgabe 3 m<-seq(from=1,to=6,by=1) i<-c(100,150,230,340,500,750) n<-length(i) sum((m-mean(m))^2) sum((i-mean(i))^2) sum((i-mean(i))*(m-mean(m))) #Lineare Regression lm(i ~ m) reg <- lm(i ~ m); reg$coefficients summary(reg) #fit i.hat <- reg$coefficients[1] + reg$coefficients[2]*m #Residuen resid <-( i - i.hat)/sqrt((sum((i-i.hat)^2))/n) #resid <-i - i.hat #Lin Reg nach Trafo i2<-log(i) reg2 <- lm(i2 ~ m); reg2$coefficients summary(reg2) i2.hat <- reg2$coefficients[1] + reg2$coefficients[2]*m #resid2 <-( i2 - i2.hat)/sqrt((sum((i2-i2.hat)^2))/n) resid2 <- i2 - i2.hat #Plot observed + fitted par(mfrow=c(2,2), lwd=1.5, font.axis=2, bty="n", ps=11, cex.axis=1) plot(m, i, pch=16, cex=1.5, ylab="Index", xlab="Monat") abline(reg$coefficients[1], reg$coefficients[2], col="red") plot(m, resid, cex=1.5, ylim=c(-2,2), pch=16,ylab="Residuen", xlab="Monat") abline(h=0, col="red") plot(m, i2, pch=16, cex=1.5, ylab="ln-Index", xlab="Monat") abline(reg2$coefficients[1], reg2$coefficients[2], col="red") plot(m, resid2, cex=1.5, pch=16,ylab="ln-Residuen", ylim=c(-2,2), xlab="Monat") abline(h=0, col="red") #vgl. Bestimmtheitsmaß cor(i,m)^2 cor(i2,m)^2 #zu Fuss: n<-length(m) b1<-(n*sum(i*m)-sum(i)*sum(m))/(n*sum(m^2)-(sum(m))^2) b0<-(sum(i)-b1*sum(m))/n b0;b1 iM<-sum(i)/n R2<-1-sum(resid^2)/sum((i-iM)^2) ``` ![](/files/-LMcy0bpOB0tPNYc1cBy) ## Vorlesung R-skript ```r ################# #Vorlesung 2 ################# #Nominalskaliert ############################################ # Verkehrsmittel absolut <- c(15,193,11,24,10) names(absolut) <- c("DB","BVG","Pkw","velo","s") absolut anzahl <- sum(absolut) anzahl relativ <- absolut / anzahl round(relativ, 2) prozent <- relativ * 100 round(prozent, 1) # Torten- und Balkendiagramm par(mfrow=c(1,2), lwd=1.5, font.axis=2, bty="n", ps=11, cex.axis=1) pie(absolut, labels=c("DB","BVG","Pkw","velo","s")) barplot(absolut, names.arg=c("db","bvg","a","velo","s"),col="blue",density=c(5,10,15,20,25), ylim=c(0,200)) # Pareto- und Balkendiagramm par(mfrow=c(1,2), lwd=1.5, font.axis=2, bty="n", ps=11, cex.axis=1) barplot(sort(prop.table(absolut), decreasing=TRUE)) barplot(absolut, names.arg=c("db","bvg","a","velo","s"),col="blue",density=c(5,10,15,20,25), ylim=c(0,200)) ############################################ # #Ordinalskaliert #wichtig: 1. geordnet darstellen, 2. alle Merkmalsausprägungen mit aufnehmen (auch die, die nicht vorkommen) # ############################################ stufen <- c(rep(1, 21), rep(2, 70), rep(3, 87), rep(4, 67), rep(5, 37),rep(6, 0)) st <- table(stufen) st prop.table(st) plot(st) barplot(st) #Faktor bilden stufenF<-factor(stufen,levels=1:6) par(mfrow=c(1,3)) plot(table(stufenF),col="blue",density=c(5,10,15,20,25,30)) barplot(table(stufenF), col="blue",density=c(5,10,15,20,25,30)) barplot(st, col="blue") par(mfrow=c(1,1)) table(stufenF)/sum(stufen) prop.table(st) cumsum(prop.table(st)) par(mfrow=c(1,1)) plot.ecdf(stufen) plot(cumsum(table(stufen))) ############################## #Vorlesung 3 # ############################## # # Pareto-Diagramm # ################################## defect <- c(12, 2, 32, 4, 19, 9, 1) names(defect) <- c("A", "B", "C", "D", "E", "F", "G") par(mfrow=c(1,1), lwd=2.5, font.axis=2, bty="u",ps=10, cex.axis=0.8, cex=1.3) library(qcc) # spez. Funktion pareto.chart() inlibrary(qcc) pareto.chart(defect, ylab = "Fehler-Häufigkeit", main =" ",xlab="Fehler-Ursache", las=3, col="blue") #alternativ tabelle<-as.table(defect) barplot(tabelle) barplot(sort(tabelle, decreasing=TRUE)) #sortiert nach Häufigkeit barplot(sort(prop.table(tabelle), decreasing=TRUE)) ######################################################################### #Quantile c1<-c(8,6,4,1,3) c2<-sort(c1) quantile(sort(c1),probs=c(0.2,0.8),type=2) quantile(c2,probs=0.2,type=2) ?quantile # Quartile vor <- c(3, 4, 6, 4, 8, 9, 2, 7, 10, 7, 5, 6, 5 ) vsort <- sort(vor); n <- length(vsort) Q1 <- vsort[floor((n+1)*0.25)]; Q1 Q2 <- vsort[floor((n+1)*0.50)]; Q2 Q3 <- vsort[floor((n+1)*0.75)]; Q3 median(vor); quantile(vor, c(0.25, 0.50, 0.75),type=2) boxplot(vor,col="blue") ######################################################################### # #Bivariate # absolut <- matrix(c(30, 10, 5, 40, 39, 7, 2, 22), nrow=2, byrow=T) colnames(absolut) <- c("A","B", "AB","0") rownames(absolut) <- c("maennlich", "weiblich") names(dimnames(absolut)) <- c("Geschlecht","Blutgruppe"); absolut gesamt<-sum(absolut) tab<-as.table(absolut) #Mosaikplot plot(tab) addmargins(tab) prop.table(tab) #absoluten Häufigkeiten für Geschlecht margin.table(absolut, 1) #absoluten Häufigkeiten für Blutgruppe margin.table(absolut, 2) #relativen Häufigkeiten für Geschlecht = Zeilenprozente round(prop.table(absolut, 1), 3) ztab<-prop.table(tab,1) addmargins(ztab) #relativen Häufigkeiten für Blutgruppe = Spaltenprozente round(prop.table(absolut, 2), 3) ``` ```r ##################### # #Daten Mietspiegel (LMU) #R-Kurs Kap 4 # #Verzeichnis: Mathematik\Lehre\Stat 1\R ############################ library(MASS) library(lattice) daten <- read.table(file="miete03.asc", header=TRUE) attach(daten) summary(daten) sum(nmqm>0) #Balkendiagramm table(rooms) rooms barplot( table(rooms)) barplot(table(rooms),col="gray",border="blue", main="Balkendiagramm für die Variable rooms", xlab="Anzahl der Räume", ylab="Absolute Häufigkeiten", ylim=c(0,800) ) pdf("barplot2.pdf") par(cex=2) barplot(table(rooms),col="gray",border="blue", main="Balkendiagramm für die Variable rooms", xlab="Anzahl der Räume", ylab="Absolute Häufigkeiten", ylim=c(0,800) ) dev.off() #Kreisdiagramm pie(table(rooms), clockwise=T) #Boxplot boxplot(nmqm, main="Boxplot der Nettomiete pro qm", col="lightgray", border="red") # Gruppierter Boxplot boxplot(nmqm ~ rooms, col="blue", main="Nettomiete pro qm nach Anzahl der Zimmer", xlab="Anzahl der Zimmer", ylab="Nettomiete pro qm") #Histogramm par( mfrow=c(2,2) ) hist(nmqm) hist(nmqm, col = "blue", freq = FALSE, main = "", xlab = "Nettomiete") # vorgegebene Klassenzahl von 50 hist(nmqm, breaks = 50, col = "blue", freq = FALSE, main = "",xlab = "Nettomiete") # direkte Vorgabe Klassen hist(nmqm, breaks = c(0, 3, 6, 9, 12, 15, 18,25),col = "red", freq = FALSE, main = "", xlab = "Nettomiete") #Histogramm (!metrische Variablen) par( mfrow=c(2,2) ) hist(nmqm) truehist(nmqm) hist(rooms) truehist(rooms) ################################################ #Streudiagramm pairs( cbind(nm,wfl), labels=c("Nettomiete", "Wohnfläche") ) #Koplots # 2 Merkmale mit 1 Faktor coplot( nm~wfl | as.factor(rooms), panel= panel.smooth ) # 2 Merkmale mit 2 Faktor coplot( nm~wfl | as.factor(rooms) * as.factor(badextra), panel= panel.smooth ) ################################################ # #Trellis # #Bivariate # ################################################ xyplot(nm~wfl | as.factor(rooms),main="Nettomiete versus Wohnfläche gegen Anzahl Zimmer" ) xyplot( nm~wfl | as.factor(rooms) * as.factor(badextra) ) splom( cbind(nm,nmqm,wfl) ) ####### # #qq Plot # ############ # Daten n <- 200 x <- rnorm(n=n) y <- rt(n=n, df=6) z <- rt(n=n, df=2) w <- exp(rnorm(n=n)) # Referenzverteilung ist Normalverteilung par( mfrow=c(2,2) ) qqnorm(x) qqnorm(y) qqnorm(z) qqnorm(w) ################# #Klassen # ### kategorisieren nmqmk<-as.factor((nmqm>5)+(nmqm>10)) summary(nmqmk) table(nmqmk) ``` ```r ############ ##Vorlseung 5 #QQ Plot ############ # Daten n <- 200 x <- rnorm(n=n) y <- rt(n=n, df=6) z <- rt(n=n, df=2) w <- exp(rnorm(n=n)) # Referenzverteilung ist Normalverteilung par( mfrow=c(2,2) ) qqnorm(x) qqnorm(y) qqnorm(z) qqnorm(w) ## #kernel density estimate #N(0,1), N(3,1) -> bimodal data<-c(rnorm(500), rnorm(500, mean=4)) summary(data) par(mfrow=c(2,2)) hist(data) plot(density(data)) plot(density(data,bw=0.1)) plot(density(data,bw=2)) ``` ```r ########## ##Vorlseung 6 # #Lagemaße + Streuungsmaße #QQ Plot ############ # Daten n <- 200 x <- rnorm(n=n) y <- rt(n=n, df=6) z <- rt(n=n, df=2) w <- exp(rnorm(n=n)) # Referenzverteilung ist Normalverteilung par( mfrow=c(2,2) ) qqnorm(x) qqnorm(y) qqnorm(z) qqnorm(w) ### BSP: geometrischer Mittelwert dis<-c(25,-20,60,-37.5) #diskrete Renditen disR<-dis/100 mean(disR) # 7% (prod(1+disR))^(1/length(disR)) #1 kein Wachstum mean(log(1+disR)) #0 ################################################################ # Body-Mass-Index) bmi <- c(28.2,23.9,20.3,26.7,25.6,32.5,23.5,19.7,27.8,26.7,20.7,28.4,33.3) n <- length(bmi) Summe <- sum(bmi); Summe Summe/n # arithmetisches Mittel mean(bmi) # Funktion mean() ######################################################################### # modifizierte BMI-Werte; gestutzter Mittelwert bmi <- c(22.2,23.9,20.3,26.7,25.6,22.5,23.5,24.7,27.8,26.7,20.7,26.4,40.3) sort(bmi) mean(bmi) mean(bmi, trim=0.1) # geometrischer Mittelwert gehalt <- c(1.025, 1.10, 1.22) # Gehaltserhoehungen lg.gehalt <- log10(gehalt) 10^mean(lg.gehalt) # mittlere Gehaltserhoehung ######################################################################### # # harmonischer Mittelwert stueck <- c(10, 5, 8) # Kosten / Stueckzahl rez.stueck <- 1/stueck; n <- length(stueck) n / sum(rez.stueck) # mittlere Stueckzahl ######################################################################### #Streuungsmaße ######################################################################### # Standardabweichung bmi <- c(28.2,23.9,20.3,26.7,25.6,32.5,23.5,19.7,27.8,26.7,20.7,28.4,33.3) m <- mean(bmi) saq <- (bmi - m)^2 sqrt(sum(saq)/(n-1)) sd(bmi) # Funktion sd() ######################################################################### # x+-sigma: Errorbar-Plot; Abbildung 3.6 library(Hmisc) par(mfrow=c(1,1), lwd=2.5, font.axis=2, bty="n", ps=11, cex.axis=1) set.seed(1) set.seed(1) x <- 1:5 y <- c(5, 6, 9, 4, 6) + rnorm(5) delta <- rnorm(5, sd=1.5) errbar( x, y, y + delta, y - delta , xlim=c(0.5,5.5), ylim=c(2,10), cex=1.5, xlab=" ", ylab=expression(bar(x)%+-%s)) ``` ```r ####### # VL 9 # Konzentrationsmaße ##### gini <- function(x, y) { # Funktion zum Gini-Index area <- 0 # Trapezregel for (i in 2:n+1) area <- area + 0.5*((x[i]-x[i-1])*(y[i]+y[i-1])) gini <- 1 - 2*area; round(gini, 3) # Gini-Index } ######################################################################### x <- c(2, 8, 10, 15, 20, 45); n <- length(x) u <- c(0, (1:n)/n) # Abszisse - relativer Index v <- c(0, (cumsum(x) / sum(x))) # Odinate - kumulierte rel. Anteile gini(u, v) # Lorenzkurve und Gini-Index par(mfrow=c(1,1), lwd=2, font.axis=2, bty="n", ps=12) plot(u, v, type="b", cex=1.5, xlim=c(0,1), ylim=c(0,1), xlab="u", ylab="v") abline(0,1, col="red", lty=2) text(0.85, 0.30, paste("Gini-Index=",round(gini(u, v),3))) polygon(c(u,0), c(v,0), angle = -45, border=NA, density = 10) ##################################################### # BSP # x <- c(2, 2, 2, 2, 2,2); n <- length(x) u <- c(0, (1:n)/n) # Abszisse - relativer Index v <- c(0, (cumsum(x) / sum(x))) # Odinate - kumulierte rel. Anteile gini(u, v) # Lorenzkurve und Gini-Index par(mfrow=c(1,1), lwd=2, font.axis=2, bty="n", ps=12) plot(u, v, type="b", cex=1.5, xlim=c(0,1), ylim=c(0,1), xlab="u", ylab="v") abline(0,1, col="red", lty=2) text(0.85, 0.30, paste("Gini-Index=",round(gini(u, v),3))) polygon(c(u,0), c(v,0), angle = -45, border=NA, density = 10) ######################################################################### # Gini-Index; numerische Bestimmung gini_num <- function(x, unbiased =FALSE) { N <- length(x) mu <- mean(x) n <- if (unbiased) N * (N-1) else N * N ox <- x[order(x)] dsum <- drop(crossprod(2 * 1:N - N - 1, ox)) round(dsum / (mu * n), 3) } x <- c(2, 2,2,2,2,2); gini_num(x, unbiased=FALSE) x <- c(0, 0,0,0,1); gini_num(x, unbiased=FALSE) n <- length(x) u <- c(0, (1:n)/n) # Abszisse - relativer Index v <- c(0, (cumsum(x) / sum(x))) # Odinate - kumulierte rel. Anteile # Lorenzkurve und Gini-Index par(mfrow=c(1,1), lwd=2, font.axis=2, bty="n", ps=12) plot(u, v, type="b", cex=1.5, xlim=c(0,1), ylim=c(0,1), xlab="u", ylab="v") abline(0,1, col="red", lty=2) text(0.85, 0.30, paste("Gini-Index=",round(gini_num(x,unbiased=FALSE ),3))) polygon(c(u,0), c(v,0), angle = -45, border=NA, density = 10) ######################################################################### ######################################################################### # Scatterplot # Alter und Körpergröße in Kalama # ############### x <- seq(18, 29, by=1) y <- c(76.1, 77.0, 78.1, 78.2, 78.8, 79.7, 79.9, 81.1, 81.2, 81.8, 82.8, 83.5) var(x); var(y); cov(x, y) cor(x, y) tous<-data.frame(Alter=x,Größe=y) cor(tous, method="spearman") cor(tous) fm<-lm( Größe ~Alter , data = tous) par(mfrow=c(1,1), lwd=1.5, font.axis=2, bty="n", ps=10, cex.axis=1) plot(x, y, pch=16, cex=1.5, xlab="Alter (Monate)", ylab="Größe (cm)", xlim=c(17, 30), ylim=c(75, 85)) abline(fm, col="red", lty=2) ######################################################################### # x <- c(13, 17, 10, 17, 20, 11, 15) y <- c(12, 17, 11, 13, 16, 14, 15) tous<-data.frame(Alter=x,Freund=y) cor(tous, method="spearman") cor(tous) var(x); var(y); cov(x, y) cor(x, y) ######################################################################### # # Rangkorrelation # L <- c(1, 2, 2, 2, 3, 3, 4, 4); M <- c(2, 4, 1, 3, 4, 3, 4, 3) r.L <- rank(L, ties.method="average"); r.L # Rangzahlen zu x r.M <- rank(M, ties.method="average"); r.M # Rangzahlen zu y D <- r.L - r.M; n <- length(D) 1- 6*sum(D^2)/(n*(n^2-1)) # Rangkorrelationskoeffizient (Spearman) cor(r.L, r.M) # Korrelationskoeffizient aus Rangzahlen ######################################################################### # # Kendall Korrelationskoeffizient x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10) y <- c(2, 1, 5, 3, 4, 6, 7, 9, 8, 10) n <- length(x) inv <- 0; prov <- 0 for (i in 1:n) { for (j in i:n) { if (x[i]y[j]) inv <- inv + 1 if (x[i] ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.