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)

Έστω τώρα ότι θέλουμε να εκτιμήσουμε την πιθανότητα \(P(X\leq x)\) για \(X\sim N(0,1)\). Πιο κάτω γράφουμε συνάρτηση, η οποία παίρνει ως όρισμα το \(x\) και ένα \(\epsilon>0\) και επιστρέφει εκτιμήτρια Monte Carlo της πιο πάνω πιθανότητας, η οποία έχει Μέσο Τετραγωνικό Σφάλμα μικρότερο από \(\epsilon^2\). Η συνάρτηση βρίσκει πόσο μεγάλο πρέπει να είναι το πλήθος των δειγμάτων που χρησιμοποιούμε στην εκτιμήτρια, με χρήση του φράγματος \(Ν\geq \frac1{4\epsilon^2}\) που είδαμε στο Παράδειγμα Α της θεωρίας.

MC.normcdf<-function(x,eps)
{
  N=ceiling(1/(4*(eps^2)))
  return(mean(rnorm(N)<=x))
}

Δοκιμάζουμε τη συνάρτησή μας για διάφορες τιμές του \(x\) και βλέπουμε ότι τα αποτελέσματα συμφωνούν με τη συνάρτηση pnorm της R.

p0=MC.normcdf(0,0.001)
p1=MC.normcdf(-1,0.001) 
pnorm(-1)
## [1] 0.1586553
abs(p1-pnorm(-1))
## [1] 0.001608746
p2=MC.normcdf(0.5, 0.001)
pnorm(0.5)
## [1] 0.6914625
abs(p2-pnorm(0.5))
## [1] 0.0007655387

Ένας πιο σωστός τρόπος να μετρήσουμε το σφάλμα είναι με χρήση του σχετικού σφάλματος.

abs(p1-pnorm(-1))/pnorm(-1)
## [1] 0.01013989
p3=MC.normcdf(-3,0.001)
abs(p3-pnorm(-3))
## [1] 9.389803e-05
abs(p3-pnorm(-3))/pnorm(-3)
## [1] 0.06955935
p4=MC.normcdf(-5,0.001)
abs(p4-pnorm(-5))
## [1] 2.866516e-07
abs(p4-pnorm(-5))/pnorm(-5)
## [1] 1

Βλέπουμε ότι για \(x<-3\) παρόλο που το απόλυτο σφάλμα παραμένει μικρό, το σχετικό σφάλμα μεγαλώνει και τείνει προς το 1. Δηλαδή έχουμε πολύ κακή απόδοση. Ο λόγος είναι ότι για \(x<-3\), η πιθανότητα \(P(X\leq x)\) είναι πολύ μικρή, άρα ακόμα και ένα ΜΤΣ της τάξης του 0.001^2 μπορεί να είναι μεγάλο σε σχέση με την πραγματική τιμή της πιθανότητας. Το πρόβλημα με την εκτιμήτρια Monte Carlo, είναι ότι πχ για \(x=-5\), η πιθανότητα \(P(X\leq -5)=2.87\times 10^{-7}\) και άρα στα \(Ν=250000\) δείγματα που παίρνουμε για να πετύχουμε μέσο τετραγωνικό σφάλμα που είναι μικρότερο από 0.001, με μεγάλη πιθανότητα δεν υπάρχει κανένα που είναι μικρότερο του \(-5\). Θα δούμε ότι σε τέτοιες περιπτώσεις αποδίδει πολύ καλύτερα η μέθοδος Importance Sampling με κατάλληλη proposal κατανομή.

2. Importance Sampling

Θα εκτιμήσουμε την πιθανότητα \(P(X\leq x)\), όπου \(X\sim N(0,1)\) για \(x\leq-3\), με χρήση Importance Sampling. Είχαμε δει ότι η απλή μέθοδος Monte Carlo έχει κακή απόδοση σε αυτή την περίπτωση. Σύμφωνα με το Παράδειγμα Β της θεωρίας, για να βελτιώσουμε την απόδοση (να ελαττώσουμε τη διακύμανση και άρα και το μέσο τετραγωνικό σφάλμα) πρέπει να επιλέξουμε proposal κατανομή g, τέτοια ώστε για \(x<-3\) να ισχύει \(f_X(x)<g(x)\). Επιλέγουμε proposal κατανομή μια αρνητική εκθετική \(Y=-Z\) όπου \(Ζ\sim Exp(1)\), δηλαδή \(g(y)=e^y, y<0\) και \(g(y)=0, y\geq0\). Τότε \[h(y)\frac{f_X(y)}{g(y)}=1_{(-\infty,x]}(y)\frac{e^{-\frac{y^2}2-y} }{\sqrt{2\pi}}.\]

IS.normcdf<-function(x,N)
{
  if(x>=0){
    stop("x needs to be negative")
  }
  y<-log(runif(N)) #arnitiki ekthetiki me parametro 1
  return(mean((y < x)*exp(-(y^2)/2-y)/sqrt(2*pi)))
}

Εφαρμόζουμε τη συνάρτηση για να εκτιμήσουμε την \(P(X<-5)\) με \(N=250000\) που δίνει ΜΤΣ \(0.001^2\) για την απλή μέθοδο Monte Carlo.

p5<-IS.normcdf(-5,250000)
abs(p5-pnorm(-5))
## [1] 2.957133e-09
abs(p5-pnorm(-5))/pnorm(-5)
## [1] 0.01031612

Βλέπουμε ότι τόσο το απόλυτο όσο και το σχετικό σφάλμα έχουν βελτιωθεί σε τεράστιο βαθμό.

simIS<-c()
simMC<-c()

for(i in 1:1000)
{
  simIS[i]<-IS.normcdf(-5,250000)
  simMC[i]<-MC.normcdf(-5,0.001)
}
var(simMC)/var(simIS) 
## [1] 10693.22

Βλέπουμε ότι η διακύμανση στο Importance Sampling είναι 4 τάξεις μεγέθους μικρότερη από αυτή στο απλό Monte Carlo (για τον ίδιο αριθμό δειγμάτων).

3. Εφαρμογές Monte Carlo

Εφαρμογή Α, εκτίμηση μεροληψίας

Έστω \(\rho\in[-1,1]\) και \(Χ, \eta \sim N(0,1)\) ανεξάρτητες. Ορίζουμε \[Y=\rho X+\sqrt{1-\rho^2}\eta.\] Είδαμε ότι \(Corr(X,Y)=\rho\). Θα μελετήσουμε τη μεροληψία στην εκτίμηση του \(\rho\) από το δειγματικό συντελεστή συσχέτισης του Pearson \[r(X,Y)=\hat{\rho}(X,Y)=\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sqrt{\sum_{i=1}^n(X_i-\bar{X})^2\sum_{i=1}^n(Y_i-\bar{Y})^2}},\] όπου \((Χ_i, Y_i)\) ανεξάρτητα δείγματα από την \((X,Y)\).

Για ένα δοθέν \(\rho\), η πιο κάτω συνάρτηση παράγει \(N\) ανεξάρτητα δείγματα από την \(\hat{\rho}(X^{(j)},Y^{(j)})\) και εκτιμά με χρήση Monte Carlo τη μέση τιμή \(E[\hat{\rho}]\) και συνεπώς τη μεροληψία \(E[\hat{\rho}]-\rho\). Κάθε εκτιμήτρια \(\hat{\rho}\), υπολογίζεται με χρήση δείγματος μεγέθους \(n\) από τη \((X,Y)\).

MC.cor.bias <- function(N,n,r)
{
  rhat <- c()
  for(j in 1:N)
  {
    x <- rnorm(n)
    y <- r*x + sqrt(1 - r^2) * rnorm(n)
    rhat[j] <- cor(x,y) #by default einai Pearson correlation coefficient
  }
  return(mean(rhat) - r)
}

Χρησιμοποιούμε τη συνάρτηση αυτή για διάφορες τιμές του \(\rho\in[-1,1]\) για να εκτιμήσουμε τη μεροληψία της \(\hat{\rho}\) για \(n=10\), με χρήση Monte Carlo με \(N=10000\).

rvector=seq(-1,1,by=0.1)
K=length(rvector)
bias=c()
for(k in 1:K)
{
  bias[k]=MC.cor.bias(10000,10,rvector[k])
}
plot(rvector, bias, xlab="rho", ylab="estimated bias")

Κοιτάζοντας το γράφημα μπορούμε να συμπεράνουμε ότι η μεροληψία κυμαίνεται μεταξύ -0.025 και 0.025 και φαίνεται να είναι μέγιστη όταν \(|\rho|=0.5\), ενώ μηδενίζεται όταν έχουμε τέλεια γραμμική συσχέτιση \(Y=\pm Χ\) (\(\rho=\pm1\)) ή όταν \(Y\) και \(X\) ανεξάρτητα, \(Y=\eta\) (\(\rho=0\)). Eπίσης έχουμε θετική μεροληψία για αρνητικά \(\rho\) και αρνητική για θετικά \(\rho\).

Αξιολόγηση \(1-\alpha\) προσεγγιστικού διαστήματος εμπιστοσύνης για τη μέση τιμή δείγματος από την κατανομή Poisson

Θεωρούμε ανεξάρτητο δείγμα \(X_1,\dots,X_n\) από την κατανομή \(Poisson(\lambda)\) για άγνωστο \(\lambda=E[Poisson(\lambda)]\). Σε αυτό το παράδειγμα θα εκτιμήσουμε το επίπεδο εμπιστοσύνης του \(95\%\) προσεγγιστικού διαστήματος εμπιστοσύνης

\[I(X)=[\bar{X}-\frac{p_{n,0.05}\hat{\sigma}}{{\sqrt{n}}}, \bar{X}+\frac{p_{n,0.05}\hat{\sigma}}{{\sqrt{n}}}],\] όπου \(p_{n,\alpha}\) είναι το \((1-\alpha/2)\)-ποσοστημόριο της κατανομής Student με \(n-1\) βαθμούς ελευθερίας.

Για δοθέν \(\lambda\), η πιο κάτω συνάρτηση χρησιμοποιεί \(Ν\) δείγματα του \(Ι(Χ^{(j)})\), για να εκτιμήσει την πιθανότητα\[P_\lambda(\lambda\in I(X))\] με χρήση Monte Carlo. Κάθε διάστημα \(I(X^{(j)})\), υπολογίζεται με χρήση δείγματος μεγέθους \(n\) από την \(Poisson(\lambda)\).

MC.CI.lambda <- function(N,n,lambda){
  p=0
  for(j in 1:N)
  {
    x <- rpois(n,lambda)
    m <- mean(x)
    s <- sd(x)
    U <- m - qt(p = 0.975, df = n-1)*s/sqrt(n)
    V <- m + qt(p = 0.975, df = n-1)*s/sqrt(n)
    p <- p + (U<lambda)*(V>lambda)
  }
  return(p/N)
}

Χρησιμοποιούμε τη συνάρτηση για διάφορες τιμές του \(\lambda\in(0,1)\) για να εκτιμήσουμε το επίπεδο εμπιστοσύνης του \(I(X)\) για \(n=10\), με χρήση Monte Carlo με \(Ν=10000\).

lambda.vector <- seq(0.02, 1, by=0.02)
L <- length(lambda.vector)
CI.vector <- c()

for(l in 1:L )
{
  CI.vector[l] <- MC.CI.lambda(10000,10, lambda.vector[l])
}  

plot(lambda.vector, CI.vector, xlab="lambda", ylab="Confidence coefficient")

Παρατηρούμε ότι για μικρές τιμές του \(\lambda\) το επίπεδο εμπιστοσύνης είναι αρκετά μικρότερο από \(95\%\). Eπαναλαμβάνουμε για \(n=100\) και βλέπουμε ότι το επίπεδο εμπιστοσύνης είναι πιο κοντά στο επιθυμητό \(95\%\).

lambda.vector <- seq(0.02, 1, by=0.02)
L <- length(lambda.vector)
CI.vector <- c()

for(l in 1:L )
{
  CI.vector[l] <- MC.CI.lambda(10000,100, lambda.vector[l])
}  

plot(lambda.vector, CI.vector, xlab="lambda", ylab="Confidence coefficient")