# Mein erstes R-Kommando 
4+6

# Benutzung der Hilfe
help("median")


###########################
# Zuweisungsoperator: <-  #
###########################
a <- 5
a

x <- c(1, 2, 3, 4, 5)
x

woerter <- c("eins", "zwei", "drei", "vier", "fuenf")

y <- c(1:100)
y

z <-seq(from=1, to=100, by=0.5)
z

w <- rep(c(1,2,3), 10)
w
length(w)



#############################
# Indexoperator: [...]      #
#############################
w[1:4]
w[seq(from=2, to=16, by=2)]



########################
# Vektoroperationen    #
########################
x<-c(1,2,3)
y<-c(2,4,6)

x + y
10*x
x*y

x < 3
sum(x < 3)


##########################
# Programmierkonstrukte  #
##########################

# Verzweigung
if(x[2] < 4) 
{
 r <- 1
} else
{
 r <- 0
} 
r

#for-Schleife
for (i in c(1:3))
{
 print(x[i])
}

for (i in c(1:5))
{
 print(woerter[i])
}


############################
# Mathematische Operatoren #
############################
log(exp(1))

a<-9
sqrt(a)

acos(-1)

##################################
#Verteilungsfunktion: P(X <=x )  #
##################################
#empirische Verteilungsfunktion: relative Haeufigkeit von Elementen
#in einem Vektor, die kleiner oder gleich x sind
#Aufgabe: empirische Verteilungsfunktion ermitteln 
#durch Verwendung von sum bzw. cumsum:

help("cumsum")
vektor <- seq(from=2, to=10, by=0.1)
plot(vektor)

#Plotte empirische Verteilungsfunktion von vektor:
gitter <- c(2:10)
n<-length(gitter)
rel_haeufigkeit<-rep(0, n)
for (j in c(1:length(gitter)))
{
 rel_haeufigkeit[j] <- sum(vektor <= gitter[j]) / length(vektor)
}
plot(gitter, rel_haeufigkeit)


#Erzeuge eine Zufallsstichprobe aus vektor mit Zuruecklegen
vektor2 <- sample(vektor, replace=TRUE)
vektor2<-sort(vektor2)
table(vektor2)

#Plotte empirische Verteilungsfunktion von vektor2:
plot(unique(vektor2), cumsum(table(vektor2)) / length(vektor2), 
     xlab='x', ylab='rel. Häufigkeit <= x');


############################
# Einige Grafikanwendungen #
############################
x <- seq(0, 10, 0.1)
x
x_sample <- x[sample(1:length(x))]
x_sample
y <- -1.5 + 4*x + rnorm(length(x), 0, 5)
plot(x, y)

# rnorm: erzeuge normalverteilte Zufallszahlen
# runif: erzeuge gleichverteilte Zufallszahlen

hist(y, probability=TRUE); #Histogramm: relative Haeufigkeit in Klassen


# Erzeuge normalverteilte Zufallszahlen 
x <- rnorm(100)

##########################
# Deskriptive Statistik  #
##########################
plot(x) #Werte von x gegen ihren Index plotten
hist(x, probability=TRUE) #relative Haeufigkeiten von Werten, die in Klassen fallen
mean(x) #Mittelwert von x
max(x)
min(x)
range(x) #Spannweite von x (min(x), max(x))
median(x) # Wert in der Mitte der geordneten Stichprobe
IQR(x) #empirischer Interquartilsabstand

# empirische Medianberechnung von Hand:
x_geordnet <- sort(x)
x_geordnet
length(x)
0.5*(x_geordnet[50] + x_geordnet[51])
median(x)


#Stichprobenziehung: sample()
sample(x, 20, replace=FALSE) #ohne Zuruecklegen
sort(sample(x, 20, replace=TRUE)) #mit Zuruecklegen

# Gebe Auswahlwahrscheinlichkeiten vor:
x_teil <- x[1:10]; x_teil;
prob_x <- c(0.004, 0.9, 0.3, 0.001, 0.6, 0.005, 0.12, 0.1, 0.0006, 1)
sample(x_teil, replace=TRUE, prob=prob_x)


#########################################
# Mehr zu statistischen Graphiken mit R #
#########################################

# Verteilungsfunktion: P(X <= x )
# Empirische Verteilungsfunktion: relative Haeufigkeit von 
# Elementen in einem Vektor, die kleiner oder gleich x sind
# Aufgabe: empirische Verteilungsfunktion ermitteln 
# durch Verwendung von sum bzw. cumsum:

#help("cumsum")
vektor <- rnorm(50)  
sort(vektor)
sort(vektor, decreasing=TRUE)
rev(sort(vektor))-sort(vektor, decreasing=TRUE)
order(vektor)
vektor[order(vektor)] - sort(vektor)

plot(vektor)

#Plotte empirische Verteilungsfunktion von vektor:
gitter <- seq(from=-4, to=4, by=0.2)
n<-length(gitter)
rel_haeufigkeit<-rep(0, n)
for (j in c(1:length(gitter)))
{
rel_haeufigkeit[j] <- sum(vektor <= gitter[j]) / length(vektor)
}
plot(gitter, rel_haeufigkeit)

#englischer Begriff fuer empirische Verteilungsfunktion:
#empricial cumulative distribution function, kurz: ecdf
help("ecdf")
plot(ecdf(vektor), add=TRUE, col="blue")

#Erzeuge eine Zufallsstichprobe aus vektor mit Zuruecklegen
vektor2 <- sample(vektor, replace=TRUE)
vektor2<-sort(vektor2)
table(vektor2)

#Plotte empirische Verteilungsfunktion von vektor2:
windows()    #Linux: X11()
plot(unique(vektor2), cumsum(table(vektor2)) / length(vektor2))
plot(ecdf(vektor2), add=TRUE, col="red")

###############################################
# Beispiel: Histogramm und Dichte in ein Bild #
###############################################
set.seed(123)    # Setze den Startwert des 
                 # Zufallszahlengenerators

x <- rnorm(100)  # 100 normalverteilte Zufallszahlen
par(las = 1)
#Histogramm, Aussehen angepasst
hist(x, main="Histogramm und Dichte", freq=FALSE, 
col="grey", ylab="Dichte", xlim=c(-5, 5), ylim=c(0, 0.6))

# Theoretische Dichte
curve(dnorm, from=-5, to=5, add=TRUE, lwd=3, lty=2)
legend(-4.5, 0.55, legend=c("emp. Histogramm", "theor. Dichte"),
       col=c("grey", "black"), lwd=5)


########################################
# Wahrscheinlichkeitsverteilungen in R #
########################################

# Zufallszahlen gemaess gegebener Verteilung erzeugen:
# Prefix "r" - randon number generation

# Verteilungsfunktion einer gegebenen Verteilung:
# Prefix "p" - probability F_X(x) = P(X <= x)

# Dichte (stetig) bzw. Wahrscheinlichkeitsfunktion (diskret)
# Prefix "d" - density

# Quantilsfunktion ((verallg.) Umkehrfunktion 
# der Verteilungsfunktion)
# Prefix "q" - quantile

# Kommandos werden zusammengesetzt aus Prefix + 
# Verteilungsname


####################################
# Beispiele: diskrete Verteilungen #
####################################
# Diskrete Gleichverteilung: sample(c(1:n), replace=TRUE)
obergrenze <- 1000
zufall_gleich <- sample(c(1:obergrenze), replace=TRUE)
plot(zufall_gleich)
plot(sort(zufall_gleich))

# Binomialverteilung: binom
p = 0.2
n = 30
stellen <- c(0:n)
help("pbinom")
plot(stellen, pbinom(stellen, size=n, prob=p),pch=16,xlab="x", ylab="W'keit", 
main="Binomialverteilung mit Parametern n=30 und p=0.2");
points(stellen, dbinom(stellen, size=n, prob=p), pch=16,col="blue");

pbinom(stellen, size=n, prob=p)
cumsum(dbinom(stellen, size=n, prob=p))

zufallszahlen <-rbinom(n=500, size=n, prob=p)
plot(ecdf(zufallszahlen))
points(stellen, pbinom(stellen, size=n, prob=p), col="blue")

# Poissonverteilung: pois
help("dpois")
my_lambda = 2.5

plot(stellen, ppois(stellen, lambda=my_lambda),pch=16,ylim=c(-0.01, 1), 
xlab="x", ylab="W'keit", main="Poisson-Verteilung mit Intensitätsparameter lambda=2.5");
points(stellen, dpois(stellen, lambda=my_lambda), pch=16,col="blue");

zufallszahlen <-rpois(n=500, lambda=my_lambda)
plot(ecdf(zufallszahlen))
points(stellen, ppois(stellen, lambda=my_lambda), col="blue")

####################################
# Beispiele: stetige Verteilungen  #
####################################

# Exponentialverteilung: exp

help("dexp") 

lambda = 0.2 # Erwartungswert ist 1/lambda = 5

stellen <- seq(from=0, to=5, by=0.1)
plot(dexp(stellen, lambda), type="l", xlab='t', ylab='Dichte(t)');
curve(dexp(x, lambda), from=0, to=5, n=length(stellen), xlab='t', ylab='Dichte(t)')
zufallszahlen <- rexp(n=10000, rate=lambda)
mean(zufallszahlen)

argument <- 4.1
pexp(argument, rate=lambda) 
f <- function(x) dexp(x, rate=lambda)
integrate(f, lower=0, upper=argument)

# Normalverteilung: norm
mu <- 3.0
sigma <- 0.2

stellen <- seq(from=mu-3*sigma, to=mu+3*sigma, by=0.01)
plot(stellen, dnorm(stellen, mean=mu, sd=sigma), type="l",
xlab='x', ylab='Dichte(x)');

plot(stellen, pnorm(stellen, mean=mu, sd=sigma), type="l",
xlab='x', ylab='Vtfkt(x)');

qnorm(0.5, mean=mu, sd=sigma)
qnorm(0.8, mean=mu, sd=sigma)

# Stetige Gleichverteilung: unif
help("runif")

a <- 2
b <- 5
stellen <- seq(from=a, to=b, by=0.01)
1/(b-a) # = 1/3 = konstanter Wert der Dichtefunktion

plot(stellen, dunif(stellen, min=a, max=b), col="red", type="l",ylim=c(0, 1),
xlab='x', ylab='Dichte(x) und Vtfkt(x)');
curve(punif(x, min=a, max=b), from=a, to=b, col="blue", type="l", add=TRUE);

g <- function(x) punif(x, min=a, max=b)

plot(stellen, dunif(stellen, min=a, max=b), col="red", type="l",
     xlim=c(a-1, b+1), ylim=c(-0.1, 1), xlab='x', ylab='Dichte(x) und Vtfkt(x)')
curve(g, from=a-1, to=b+1, n=length(stellen), add=TRUE, col="blue")

punif(1.0, min=a, max=b) # = 0