#-------------------------------------#
#                                     #
#    Week8: Monte Carlo Simulation    #
#                                     #
#-------------------------------------#


Schauen Sie sich zunaechst die Sachen in dem 

                week8.pdf

an, dort gibt es den theoretischen Hintergrund 
zur Monte Carlo Methode und es werden die drei 
Beispiele spezifiziert, die wir hier jetzt be-
trachten wollen:



###   Start R-Session   ###



# Beispiel 1)


N = 1000000               # die Bilder sind dann etwas schoener
                          # als nur fuer N = 50000
x = rnorm(N)

lambda = 1
exactres = exp(lambda^2/2)

# Integral mit MC:

terms = exp(lambda*x)
mcsum = cumsum(terms)/(1:N)

# plot(mcsum)             # this takes longer for N = 1'000'000

# only 1000,2000,...,1'000'000 for plot:

plotindices = seq(from=1000, to=N, by=1000)

plot(mcsum[plotindices])
abline(a=exactres,b=0,col="red")


# Monte Carlo Fehler

mcerror      = exactres - mcsum
sqrtNmcerror = sqrt(1:N) * mcerror
Nmcerror     =     (1:N) * mcerror

par(mfrow=c(2,2))         # 2 times 2 plot array

plot(mcsum[plotindices])
abline(a=exactres,b=0,col="red")

plot(mcerror[plotindices])
abline(a=0,b=0,col="red")

plot(sqrtNmcerror[plotindices])
abline(a=0,b=0,col="red")

plot(Nmcerror[plotindices])
abline(a=0,b=0,col="red")



# Beispiel 2)


N = 1000000
x = rexp(N,rate=1)        # exponential-verteilte 
                          # Zufallszahlen
n = 4
exactres = factorial(n)

# Integral mit MC:

terms = x^n
mcsum = cumsum(terms)/(1:N)

# plot(mcsum)             # this takes longer for N = 1'000'000

# only 1000,2000,...,N for plot:

plotindices = seq(from=1000, to=N, by=1000)

plot(mcsum[plotindices])
abline(a=exactres,b=0,col="red")

# Monte Carlo Fehler

mcerror      = exactres - mcsum
sqrtNmcerror = sqrt(1:N) * mcerror
Nmcerror     =     (1:N) * mcerror

par(mfrow=c(2,2))
plot(mcsum[plotindices])
abline(a=exactres,b=0,col="red")
plot(mcerror[plotindices])
abline(a=0,b=0,col="red")
plot(sqrtNmcerror[plotindices])
abline(a=0,b=0,col="red")
plot(Nmcerror[plotindices])
abline(a=0,b=0,col="red")



# Beispiel 3)


N = 1000000

x = runif(N)
y = runif(N)
z = runif(N)

exactres = 1/6

# Integral mit MC:

terms = ifelse( x+y+z<=1 , 1 , 0 )         # recht nuetzlicher Befehl
mcsum = cumsum(terms)/(1:N)

# plot(mcsum)             # this takes longer for N = 1'000'000

# only 1000,2000,...,N for plot:

plotindices = seq(from=1000, to=N, by=1000)

plot(mcsum[plotindices])
abline(a=exactres,b=0,col="red")

# Monte Carlo Fehler

mcerror      = exactres - mcsum
sqrtNmcerror = sqrt(1:N) * mcerror
Nmcerror     =     (1:N) * mcerror

par(mfrow=c(2,2))
plot(mcsum[plotindices])
abline(a=exactres,b=0,col="red")
plot(mcerror[plotindices])
abline(a=0,b=0,col="red")
plot(sqrtNmcerror[plotindices])
abline(a=0,b=0,col="red")
plot(Nmcerror[plotindices])
abline(a=0,b=0,col="red")