#----------------------------------------# # # # week2a: 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 t Students t-Verteilung binom Binomial-Verteilung pois Poisson-Verteilung ... ... Eine sehr schoene und sehr kompakte, gerade mal 3 Seiten, Uebersicht gibt es in dem 5. Kapitel des Skriptes Grund- lagen der Datenanalyse mit R auf der VL-homepage, http://hsrm-mathematik.de/SS2023/semester4/Stochastik2/R1-UniGiessen-GerritEichner.pdf auf den Seiten 87-89. Weil das wirklich sehr nuetzlich ist, gibt es auch nur diese 3 Seiten nochmal hier, http://hsrm-mathematik.de/SS2023/semester4/Stochastik2/W'keitsverteilungen-in-R.pdf Zu jeder in R vorimplementierten W'keitsverteilung Vert 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 jetzt genauer an: ### Start R-Session ### # Normalverteilte Zufallszahlen: v = rnorm(1000) v head(v) tail(v) plot(v) hist(v) hist(v,breaks=100) hist(v,100) # dasselbe.. ?hist hist(v,100,probability=true) # Fehlermeldung hist(v,100,probability=True) # Fehlermeldung hist(v,100,probability=TRUE) # das geht; ist dasselbe Bild, aber mit anderer Skalierung hist(v,100,prob=TRUE) # geht auch curve(dnorm(x)) curve(dnorm(y)) # Fehlermeldung, der curve-Befehl moechte immer ein x haben.. curve(dnorm(x),from=-5,to=5) # schoene x-Achse ?curve hist(v,breaks=30,prob=TRUE) curve(dnorm(x),add=TRUE) curve(dnorm(x),add=TRUE,col="red") # now 10'000 instead of 1'000 random numbers, dann wird das Bild schoener.. x = rnorm(10000) x plot(x) hist(x) hist(x,100) hist(x,100,prob=TRUE) curve(dnorm(x),add=TRUE,col="red") # zur Uebersichtlichkeit jetzt nur n=10 Stueck: x = rnorm(10) x sortx = sort(x) sortx rankx = rank(x) rankx x # rankx[i] = k bedeuted: x[i] steht in sort(x) an k-ter Stelle, # ist also das k-te kleinste orderx = order(x) orderx # Bedeutung: sortx[i] = x[orderx[i]], also: man moechte das i-te kleinste Element # von x haben, welchen Index muss ich nehmen? # insbesondere: orderx[1] ist der Index fuer das Minimum, orderx[n] ist der Index # fuer das Maximum, testen wir das eben: testx = x[orderx] testx testx == sortx # wieder n=1000 Zufallszahlen: n = 1000 x = rnorm(n) sortx = sort(x) plot(sortx) # q%-Quantil := sortx[ q%*n ], also: q% aller Daten sind kleiner als das q%-Quantil # erstes Quartil := 25%-Quantil # drittes Quartil := 75%-Quantil # Median = zweite Quartil := 50%-Quantil summary(x) quantile( x , 0.25 ) # erste Quartil quantile( x , 0.5 ) # zweite Quartil = Median quantile( x , 0.75 ) # dritte Quartil quantile( x , 0.05 ) sortx[50] sortx[51] sortx[49] mean(x) sd(x) var(x) sd(x)^2 # Jetzt machen wir dasselbe nochmal mit einer diskreten Verteilung, # nehmen wir etwa die Binomial-Verteilung: # # Bedeutung der Binomialverteilung: EIN Zufallsexperiment ist: # wir wuerfeln 100 mal (n=100) und notieren die Anzahl der Einsen oder Zweien =: X (also p=2/6) # dann: # B_{n,p}(k) = (n ueber k) * p^k * (1-p)^{n-k} = Prob[ X = k ] # # dieses Zufallsexperiment machen wir jetzt 1000 mal: X = rbinom(1000, 100, 2/6) # allgemein: x = rVerteilungsname( Anzahl der Zufallszahlen , Verteilungsparameter ) X plot(X) hist(X) hist(X,100) hist(X,1000) # Wir wollen das wieder mit der exakten Formel fuer die Binomialverteilung # vergleichen. Allgemein: 2 Faelle: # # (i) diskrete Verteilung: Prob[ X = k ] : dVert( k , Verteilungsparameter ) # (ii) kontinuierliche Verteilung: Prob[ X in [x,x+dx) ] : dVert( x , Verteilungsparameter ) * dx # ("d" steht also fuer "density") # Testen wir das eben fuer die Binomialverteilung # B_{n,p}(k) = (n ueber k) p^k * (1-p)^(n-k) k = 30 n = 100 p = 2/6 prob = dbinom( k , n , p ) prob # Binomialkoeffizienten mit choose(n,k): choose(4,2) # 4*3/(1*2) choose(6,3) # 6*5*4/(1*2*3) choose(n,k) # (100 ueber 30) = 100*99*...*71/(1*2*...*30) prod(71:100)/prod(1:30) # ok, dasselbe.. checkprob = choose(n,k) * p^k * (1-p)^(n-k) checkprob prob checkprob == prob all.equal(prob,checkprob) # up to numerical noise.. plot( dbinom(1:100, 100, 2/6) ) plot( 1000*dbinom(1:100, 100, 2/6) ,col="red" ) hist(X,1000,add=TRUE) # Jetzt wieder 10'000 Zufallszahlen anstatt von 1'000: X = rbinom(10000, 100, 2/6) plot(X) hist(X) hist(X,100) hist(X,1000) hist(X,10000) plot(dbinom(1:100,100,2/6)) plot(10000*dbinom(1:100,100,2/6),col="red") hist(X,10000,add=TRUE) # Zur Uebersichtlichkeit jetzt wieder n=10: X = rbinom(10, 100, 2/6) X sortX = sort(X) sortX rankX = rank(X) rankX # rankX[i] = k bedeuted: X[i] steht in sort(X) an k-ter Stelle, # ist also das k-te kleinste. # ok, es gibt dann eine Logik wenn das nicht eindeutig ist, # weil ein paar von den X[i]'s gleich sein koennen.. orderX = order(X) orderX # Bedeutung: sortx[i] = x[orderx[i]], also: man moechte das i-te kleinste Element # von x haben, welchen Index muss ich nehmen? # insbesondere: orderx[1] ist der Index fuer das Minimum, orderx[n] ist der Index # fuer das Maximum, testen wir das eben: testX = X[orderX] testX testX == sortX # wieder n=1000: n = 1001 X = rbinom(n, 100, 2/6) sortX = sort(X) plot(sortX) # q%-Quantil := sortx[ q%*n ], also: q% aller Daten sind kleiner als das q%-Quantil # erstes Quartil := 25%-Quantil # drittes Quartil := 75%-Quantil # Median = zweite Quartil := 50%-Quantil summary(X) quantile( X , 0.25 ) quantile( X , 0.5 ) quantile( X , 0.75 ) quantile( X , 0.05 ) sortX[50] sortX[51] sortX[52] mean(X) sd(X) sqrt(100*1/3*2/3) var(X) sd(X)^2 # performance: # 1 Million numbers: x = rnorm(1000000) sx = sort(x) # 10 Million numbers: x = rnorm(10000000) sx = sort(x) # quite fast X = rbinom(10000000,100,2/6)