Statistik-1
My super prof. teach me this stuff. Many Thank for Prof. Dr. Beate Bergter HTW Berlin
Ü1. Statistische Grundbegriffe
##########################
#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
### 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))### 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"#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")#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)##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
#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") #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) #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#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
#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") # 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) # 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)#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
#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)
#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)))#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#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#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
#################
#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")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#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)
Vorlesung R-skript
#################
#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)#####################
#
#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)############
##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))##########
##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))#######
# 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]<x[j] & y[i]>y[j]) inv <- inv + 1
if (x[i]<x[j] & y[i]<y[j]) prov <- prov + 1
}
}
inv; prov
r.tau <- 1 - 4*inv/(n*(n-1)); r.tau
#########################################################################
cor(x, y , method = "kendall") # spez. Funktion cor()
cor.test(x, y, method="kendall") # spez. Funktion cor.test()
#########################################################################
#
# partielle Rangkorrelation
z <- c(1,2,3,4,5,6) # einfuehrendes Beispiel
x <- c(3,1,4,2,6,5)
y <- c(4,2,1,6,3,5)
tauXZ <- round(cor(x, z, method="kendall"), 4)
tauYZ <- round(cor(y, z, method="kendall"), 4)
tauXY <- round(cor(x, y, method="kendall"), 4)
tauXY.Z <- (tauXY - tauXZ*tauYZ) / sqrt((1-tauXZ^2)*(1-tauYZ^2))
cbind(tauXZ, tauYZ, tauXY, tauXY.Z)
# partielle Rangkorrelation
I <- c(1,2,3,4,5,6,7,8,9,10) # Intelligenz
A <- c(1,4,5,6,2,7,3,9,8,10) # Mathematik
B <- c(4,1,3,5,2,6,7,10,9,8) # Musik
tauAI <- round(cor(A, I, method="kendall"), 4)
tauBI <- round(cor(B, I, method="kendall"), 4)
tauAB <- round(cor(A, B, method="kendall"), 4)
tauAB.I <- (tauAB - tauAI*tauBI) / sqrt((1-tauAI^2)*(1-tauBI^2))
cbind(tauAI, tauBI, tauAB, tauAB.I)Last updated
