Statistik-1

My super prof. teach me this stuff. Many Thank for Prof. Dr. Beate Bergter HTW Berlin

Ü1. Statistische Grundbegriﬀe

##########################
#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äuﬁgkeiten

### 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)

PreviousR packages NextStatistik-2

Last updated 6 years ago

Statistik-1

My super prof. teach me this stuff. Many Thank for Prof. Dr. Beate Bergter HTW Berlin

Ü1. Statistische Grundbegriﬀe

##########################
#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äuﬁgkeiten

### 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)

PreviousR packages NextStatistik-2

Last updated 6 years ago