R Ergänzungen zur Vorlesung

Irrfahrt 

### R examples

#

# Example 0 - Binomial distribution

k=0:10

plot(k,dbinom(k,10,0.5))

points(k,dbinom(k,10,0.3),col="red")

#DIY: dhyper !

#

# Example 1 - Irrfahrt.

N=10000

omega1 = sample(c(-1,1),size=N,replace=TRUE)

omega2 = sample(c(-1,1),size=N,replace=TRUE)

# Der Pfad ist die kumulative Summe

S1 = cumsum(omega1)

S2 = cumsum(omega2)

# Für den Plot ist es nützlich min und max zu kennen

ymin = min(c(S1,S2))

ymax = max(c(S1,S2))

# Plot

plot(S1,type="l",main="Irrfahrt",ylim=c(ymin,ymax),xlab="n",ylab="S_n")

lines(S2,col="red")

# ---- dann setze N auf 10000

#

# Example 2 - Letzter Besuch der Null.

# Wir wiederholen obiges Beispiel und bestimmen den letzten

# Zeitpunkt an welchem S die Null erreicht hat

T1 = max(which ( S1 == 0))

T2 = max(which ( S2 == 0))

lines(c(0,N),c(0,0),lty=2)

lines (c(T1,T1),c(ymin,ymax),lty=2,col="green")

lines (c(T2,T2),c(ymin,ymax),lty=2,col="green")

# Nun versuchen wir die Verteilung dieser Zeiten empirisch zu schätzen

N=200

N2=10000

T = rep(0,N2)

for (i in (1:N2)) {

S = cumsum ( sample(c(-1,1),size=N,replace=TRUE) )

T[i] = max( which (S== 0))

}

hist (T,breaks=seq(-1,201,2),main="")

Stetige Verteilungen

### R examples -2, Stetige Verteilungen

#

# Example 1 - Normal Verteilung

X=rnorm(1000)

hist(X,lwd=2)

lines(seq(-4,4,by=0.1),sqrt(2*pi)*180*dnorm(seq(-4,4,by=0.1)),col=2,lwd=2)

# Dichte und Stammfunktion

x=seq(-4,4,by=0.1)

plot(x,dnorm(x),ylim=c(0,1),type="l",lwd=2)

# Dann die Stammfunktion

lines(x,pnorm(x),col=2,lwd=2)

# Nun die Inverse der Verteilungsfunktion

par(mfrow=c(1,2))

x=seq(-4,4,by=0.1)

plot(x,dnorm(x),ylim=c(0,1),type="l",lwd=2)

lines(x,pnorm(x),col=2,lwd=2)

plot(seq(0,1,by=0.01),qnorm(seq(0,1,by=0.01)),lwd=2,type="l")

# Wir illustrieren die Generierung von Zufallsvariablen durch die Gleichverteilung

# F^{-1}(y) = - \alpha^{-1} \Log( 1-y)

alpha=2

U=runif(100)

X=-1/alpha*log(1-U)

plot(U,X)

# Histogramm - es funktioniert !

hist(-1/alpha*log(runif(500)),main="Histogramm")

lines(seq(0,3,by=0.01),dexp(seq(0,3,by=0.01),rate=alpha)*130,col=2,lwd=2) 

Hypothesentest

x = seq(0,1,by=0.01)

# Hier ist k = 2

k=2

plot(x,pbinom(k,20,x,lower.tail=FALSE),type="l",main="Güte",ylab="",xlab="theta")

lines(c(0.2,0.2),c(0,1),col="red",lty=2)

lines(c(0,1),c(0.05,0.05),col="green",lty=2)

k=4

lines(x,pbinom(k,20,x,lower.tail=FALSE),lty=2)

k=6

lines(x,pbinom(k,20,x,lower.tail=FALSE),lty=3)

k=7

lines(x,pbinom(k,20,x,lower.tail=FALSE),lty=4,col="green",lwd=2)

# Finde p-Wert wenn x=10 beobachtet wurde.

# Hier ist theta_0 = 0.2

1-pbinom (10,20,0.2)