#----------------------------------------#
#                                        #
#     week5: W'keitsverteilungen und     #
#          Zufallszahlen in R            #
#                                        #
#----------------------------------------#


In R sind eine Reihe von ueblichen Wahrscheinlichkeits-
verteilungen vorimplementiert. Etwa:

R-Name        mathematische Bezeichnung

norm          Normalverteilung
unif          Gleichverteilung oder Uniforme Verteilung
binom         Binomial-Verteilung
pois          Poisson-Verteilung 
t             t-Verteilung
...
...

Eine sehr schoene und sehr kompakte Uebersicht gibt es auf 
den folgenden 3 Seiten:

https://hsrm-mathematik.de/SS2025/Wirtschaftsmathematik3/W'keitsverteilungen-in-R.pdf 


Zu jeder in R vorimplementierten W'keitsverteilung

 Vert in { norm , unif , binom , pois , t , ... }
 
gibt es die 4 Funktionen 

  rVert()    -> erzeugt Zufallszahlen 
  dVert()    -> Dichte der Verteilung
  pVert()    -> Integral oder Summe ueber die Dichte, kummulierte Verteilungsfunktion
  qVert()    -> das Inverse von pVert(), die Quantil-Funktion

Schauen wir uns das genauer an:



### Start R-Session ###


# Beispiel 1: Gleichverteilung
#
# Plotten wir die Dichte einer Gleichverteilung auf dem Intervall, sagen wir, 
# [2,6]:  dunif(x,2,6) 
# Anstatt des plot()-Befehls nehmen wir den curve()-Befehl, wenn man damit 
# Funktionen von mehreren Variablen plotten moechte, das dunif() hat 3 Argumente, 
# muss man ein x in das Argument schreiben, welches auf der waagerechten  Achse, 
# der x-Achse, abgetragen werden soll:


curve( dunif(x,2,6), from=0, to=10 )

curve( dunif(y,2,6), from=0, to=10 )        # liefert Fehlermeldung, der curve()-Befehl 
                                            # moechte immer ein x


# Integral ueber dunif = punif:

curve( punif(x,2,6), from=0, to=10 ) 

# die Dichte dunif diesem Bild hinzufuegen, etwa in rot:

curve( dunif(x,2,6), col="red", add=TRUE )



# Quantil-Funktion, die ist immer nur auf dem Intervall [0,1] definiert:
# (schauen Sie im week5.pdf nach fuer mehr Erklaerungen zu dieser Funktion)

curve( qunif(x,2,6), from=0, to=1 )

curve( 4*x+2, col="red", add=TRUE )         # das passt



# 1000 auf dem Intervall [2,6] gleichverteilte Zufallszahlen:

z = runif(1000, min=2, max=6 )
z
plot(z)
plot(z, ylim=c(0,10), main="alle Zahlen sind zwischen 2 und 6" )


# Histogramm, zeigt im Wesentlichen (bis auf Skalierung) das dunif:

hist(z)                                     # ok, gleichverteilt..


# wie kommt das erste Rechteck zustande? Das ist die Anzahl der Zufallszahlen 
# zwischen 2 und 2.5, ueberpruefen wir das:

test = ifelse( z>2 & z<2.5 , 1 , 0 )
test
sum(test)                                   # ok, passt in etwa..


# wenn wir eine W'keitsverteilung sehen wollen, muss die Flaeche aller 
# Rechtecke gleich 1 sein, geht mit prob = TRUE:

hist(z, prob=TRUE )
hist(z, xlim=c(0,10), prob=TRUE )

# Vergleich mit theoretischer Dichte: dunif(x, min=2, max=6 )

curve( dunif(x, min=2, max=6 ), add=TRUE, col="red" )       # das passt.. 





# Beispiel 2: Normalverteilung


# W'keitsdichte:  dnorm(x, mean=mu, sd=sigma )

curve( dnorm(x, mean=0, sd=1), from=-10, to=10 )

# plotten wir das fuer verschiedene Mittelwerte mu, jedes mu 
# bekommt eine eigene Farbe:

mu = seq(-5, 5, by=0.5 )
n = length(mu)

mycol = rainbow(n)      # n Farben, Regenbogen-maessig geordnet

for( i in 1:n )
{
  curve( dnorm(x, mean=mu[i], sd=1), col=mycol[i], add=TRUE )
}


# jetzt variieren wir das sigma bei festem mu=0:

sigma = seq(0.5, 5, by=0.1 )
n = length(sigma)
mycol = rainbow(n)

curve( dnorm(x, mean=0, sd=0.5), from=-10, to=10 )
for( i in 1:n )
{
  curve( dnorm(x, mean=0, sd=sigma[i] ), col=mycol[i], add=TRUE )
}

# ok, man koennte die Aufloesung etwas erhoehen, da ist nicht wirklich eine Spitze..



# Integral ueber die Dichte: pnorm(x, mean=mu, sd=sigma )

mu = seq(-5, 5, by=0.5 )
n = length(mu)
mycol = rainbow(n)    

curve( pnorm(x, mean=0, sd=1), from=-10, to=10 )
for( i in 1:n )
{
  curve( pnorm(x, mean=mu[i], sd=1), col=mycol[i], add=TRUE )
}

# packen wir noch die Dichten mit dazu:
for( i in 1:n )
{
  curve( dnorm(x, mean=mu[i], sd=1), col=mycol[i], add=TRUE )
}


# fuer variables sigma:

sigma = seq(0.5, 5, by=0.1 )
n = length(sigma)
mycol = rainbow(n)

curve( pnorm(x, mean=0, sd=0.5), from=-10, to=10 )
for( i in 1:n )
{
  curve( pnorm(x, mean=0, sd=sigma[i] ), col=mycol[i], add=TRUE )
}



# n = 10000 (mu=5,sigma=1)-normalverteilte Zufallszahlen:

z = rnorm(10000, mean=5, sd=1 )
z
plot(z)
plot(z, ylim=c(0,10) )


# Histogramm, zeigt die Verteilung:

hist(z)

# wenn wir eine W'keitsverteilung sehen wollen, muss die Flaeche aller 
# Rechtecke gleich 1 sein, geht mit prob = TRUE:

hist(z, prob=TRUE )
hist(z, breaks=40, prob=TRUE )


# Vergleich mit theoretischer Dichte: dnorm(x, mean=5, sd=1 )

curve( dnorm(x, mean=5, sd=1 ), add=TRUE, col="red" )



# Uebungsaufgabe: 

# Wie gross ist die W'keit, dass eine (mu=5,sigma=1)-normalverteilte Zufallszahl 
# zwischen 3 und 4 liegt?

pnorm(4,mean=5,sd=1) - pnorm(3,mean=5,sd=1)          # = 0.1359

# bei 10000 Zufallszahlen muessten dann in etwa 

10000*0.1359               # = 1359 

# also in etwa 1360 Zahlen zwischen 3 und 4 liegen, checken wir das mal: 

test = ifelse(z>3 & z<4 , 1 , 0 ) 
test
sum(test)                  # ok, so mehr oder weniger, sind ja Zufallszahlen..




# Schauen wir uns noch die Quantil-Funktion an:  qnorm(x, mean=5, sd=1 ) 
# Quantil-Funktionen sind immer auf dem Intervall [0,1] definiert:  
 
curve( qnorm(x, mean=5, sd=1 ), from=0, to=1 )

# etwas hoehere Aufloesung (wird jetzt an n=10000 Stellen berechnet, default ist 100):

curve( qnorm(x, mean=5, sd=1 ), from=0, to=1, n=10000 )


sz = sort(z)               # ordnet die Zufallszahlen..
sz
plot(sz)                   # sieht sehr aehnlich aus..

# versuchen wir, das exakt zu matchen:  

curve( qnorm(x/10000, mean=5, sd=1), add=TRUE, col="red", lwd=3 )         # das passt!

# lwd = line width, default ist 1, ist etwas duenn..





# Beispiel 3: Diskrete Verteilung
#
# Versuchen wir noch, fuer eine diskrete Verteilung das Histogramm 
# von Zufallszahlen mit der theoretischen Dichte zu matchen: 

N = 10000
z = rbinom(N, size=100, prob=1/6 )
plot(z)
hist(z)
hist(z, prob=TRUE )

curve(dbinom(x,size=100,prob=1/6), add=TRUE, col="red" )       

# liefert Warnungen und sieht etwas komisch aus..  ..curve ist fuer das Plotten 
# von kontinuierlichen Funktionen, dbinom ist diskret..


hist(z, prob=TRUE )
kk = min(z):max(z)
kk
points(dbinom(kk,size=100, prob=1/6), col="red" )       # ..geht in die richtige Richtung..


hist(z, prob=TRUE )
points(kk, dbinom(kk,size=100,prob=1/6), col="red" )    # ..schon ganz gut..


nbreaks = max(z) - min(z)
nbreaks
hist(z, prob=TRUE, breaks = nbreaks )
points(kk, dbinom(kk,size=100,prob=1/6), col="red" )

# ok, ist fast perfekt, aber die Rechtecke sind noch nicht 'in der Mitte'.. 
# koennen wir so fixen:

xbreak = seq(from=min(z)-0.5, to=max(z)+0.5, by=1 )
xbreak
hist(z, prob=TRUE, breaks = xbreak )
points(kk, dbinom(kk,size=100,prob=1/6), col="red" )    # ok, das passt!