1. Απλή προσομοίωση Monte Carlo

Θα εκτιμήσουμε τη μέση τιμή \(E[sin(X)^2]\), όπου \(X\sim N(0,1)\). Η συγκεκριμένη μέση τιμή δεν υπολογίζεται αναλυτικά και θα την εκτιμήσουμε με τη μέθοδο Monte Carlo, \[z_{N}^{MC}=\sum_{j=1}^Nsin(X_j)^2,\] όπου \(X_j\) ανεξάρτητα και ισόνομα δείγματα \(Ν(0,1)\). Η πιο κάτω συνάρτηση παίρνει ως όρισμα το \(N\) και επιστρέφει τη \(z_N^{MC}\).

my.MCsin<-function(N)
{
  return(mean(sin(rnorm(N))^2))
}

my.MCsin(1000)
## [1] 0.423111

Παίρνουμε 10000 δείγματα της εκτιμήτριας Monte Carlo για \(N=1000\) και \(Ν=10000\) και κάνουμε τα αντίστοιχα ιστογράμματα.

simMC1<-c()
simMC2<-c()

for(i in 1:10000)
{
  simMC1[i]<-my.MCsin(1000)
  simMC2[i]<-my.MCsin(10000)
}
hist(simMC1, breaks= seq(min(simMC1)-0.002, max(simMC1)+0.002, by=0.002), main="samples of MC estimator with N=1000")

hist(simMC2, breaks= seq(min(simMC1)-0.002, max(simMC1)+0.002, by=0.002), main="samples of MC estimator with N=10000")

Βλέπουμε ότι όπως περιμέναμε από τη θεωρία, καθώς το N μεγαλώνει τα δείγματα συγκεντρώνονται γύρω από μια κοινή μέση τιμή και η διακύμανση ελαττώνεται .

var(simMC1)/var(simMC2)
## [1] 10.25663

Μάλιστα η διακύμανση φαίνεται να ελαττώνεται πράγματι με ρυθμό 1/N όπως προβλέπει η Πρόταση 4.2.

Τα δείγματα από την εκτιμήτρια Monte Carlo, φαίνονται να ακολουθούν κανονική κατανομή. Αυτό οφείλεται στο Κεντρικό Οριακό Θεώρημα (παρόμοια με απόδειξη Πρότασης 4.3).

Πράγματι, μπορούμε να κανονικοποιήσουμε το δείγμα από τις δύο εκτιμήτριες (για Ν=1000 και Ν=10000) αντίστοιχα και να επιβεβαιώσουμε ότι ακολουθεί τυπική κανονική κατανομή.

std1<-scale(simMC1)
std2<-scale(simMC2)
ks.test(std1, 'pnorm')
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  std1
## D = 0.0069824, p-value = 0.7142
## alternative hypothesis: two-sided
ks.test(std2, 'pnorm')
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  std2
## D = 0.006105, p-value = 0.8501
## alternative hypothesis: two-sided
qqnorm(std1)
qqline(std1)

qqnorm(std2)
qqline(std2)