#---------------------------#
#    Loesungen UeBlatt2     #
#---------------------------#

 
#  Aufgabe 1  
#
#  1b)

x = runif(1001)
max(x)
min(x)

which.max(x)
x[which.max(x)]

which.min(x)
x[which.min(x)]

xx = sort(x)
head(xx)
tail(xx)


#  1c)

mean(x)              # Mittelwert
var(x)               # Varianz
sd(x)                # Standardabweichung

quantile(x,0.01)
xx[11]
quantile(x,0.99)
xx[991]
quantile(x,0.25)     # erstes Quartil
xx[251]
quantile(x,0.75)     # drittes Quartil
xx[751]
quantile(x,0.5)      # Median
xx[501]


#  1d)

plot(x)
plot(x, ylim=c(-2,2) )



#  Aufgabe 2  
#
#  2a)

N = 2.53*10^(18)
p = 12/N
N
p


#  2b)

k = 0:25
k
BNpk = dbinom(k,N,p)
plot(k,BNpk)


#  2c)

mu = 12
k = 0:25
Pmuk = dpois(k,mu)
plot(k,Pmuk)

# wir fuegen die Punkte aus (2b) hinzu:
points(k,BNpk, col="red", pch="x")    # pch = plot character


#  2d)

BNpk == Pmuk              # ein paar TRUE, das meiste FALSE
all.equal(BNpk,Pmuk)      # TRUE


#  2e)
#
# wir testen den Wert fuer k=12:

N = 2.53*10^(18)
p = 12/N
prob = dbinom(12,N,p)
prob                                         # 0.1144

testprob = choose(N,12)*p^(12)*(1-p)^(N-12)
testprob                                     # 18613.93, wrong result..

testprob2 = choose(N,12)*p^(12)
testprob2                                    # 18613.93, dasselbe wie testprob..
(1-p)^(N-12)                                 # = 1

# Was klappt hier nicht? p darf nicht zu klein sein: 
# p muss groesser sein als  

.Machine$double.neg.eps                      # ist bei mir 1.110223e-16

# damit 1-p als von 1 verschieden erkannt wird: Etwa fuer
N = 10^(15)
p = 12/N
p                                            # 1.2e-14, ist jetzt groesser
                                             # als .Machine$double.neg.eps
prob = dbinom(12,N,p)                        
prob                                         # 0.1144
testprob = choose(N,12)*p^(12)*(1-p)^(N-12)
testprob                                     # 0.1155, also bei sehr gros * sehr klein 
                                             # muss man etwas vorsichtig sein: 
choose(N,12)                                 # sehr sehr gross         
p^(12)                                       # sehr sehr klein
(1-p)^(N-12)                                 # klein



#  Aufgabe 3  
#
#  3a)

N = 25
k = 0:N
p = 0.1
plot(k, dbinom(k,N,p) )

muN = N*p
sigN = sqrt(N*p*(1-p))

curve( dnorm((x-muN)/sigN)/sigN, add=TRUE, col="red")

# oder, was dasselbe ist:
curve( dnorm(x,muN,sigN), add=TRUE, col="blue")

# es gilt also: 
# dnorm(x, muN, sigN) = dnorm((x-muN)/sigN, 0, 1) / sigN


# Jetzt machen wir 4 Diagramme in einem Plotfenster:

par( mfrow = c(2,2) )       # 2 x 2 Plot-Array, etwa: "make frame by row"

for( N in c(25,50,100,200) )
{ 
  # derselbe code wie oben:
  k = 0:N
  p = 0.1
  plot(k, dbinom(k,N,p) )
  muN = N*p
  sigN = sqrt(N*p*(1-p))
  curve( dnorm((x-muN)/sigN)/sigN, add=TRUE, col="red")
}


#  3b)
#
# Wir machen gleich die 4 Diagramme in einem Plotfenster:

par( mfrow = c(2,2) )

for( N in c(25,50,100,200) )
{ 
  k = 0:N
  p = 0.6
  plot(k, dbinom(k,N,p) )
  muN = N*p
  sigN = sqrt(N*p*(1-p))
  curve( dnorm((x-muN)/sigN)/sigN, add=TRUE, col="red")
}