Σε αυτή τη σειρά θα λύσουμε ασκήσεις προσομοίωσης για ολοκλήρωση Monte Carlo και άλλα εισαγωγικά παραδείγματα.
Σε όλες τις ασκήσεις η διαδικασία είναι ως εξής:
Πρέπει να εκτιμήσουμε μια μέση τιμή \(\theta = E(X)\) όπου \(X\) τυχαία μεταβλητή από μια δεδομένη κατανομή (γενικά πολυδιάστατη). Για την εκτίμηση η διαδικασία προσομοίωσης ακολουθεί γενικά τα παρακάτω βήματα.
Σχεδιάζουμε μια συνάρτηση που δημιουργεί μια παρατήρηση από την κατανομή της \(X\) (ένα σενάριο ή μια επανάληψη).
Καλούμε τη συνάρτηση αυτή \(N\) φορές και δημιουργούμε ένα ψευδοτυχαίο δείγμα \(N\) από την κατανομή της \(X\). Αποθηκεύουμε αυτό το δείγμα σε ένα διάνυσμα ή πίνακα κατάλληλης διάστασης.
Από το δείγμα υπολογίζουμε τη δειγματική μέση τιμή, τη δειγματική τυπική απόκλιση και με βάση αυτά ένα διάστημα εμπιστοσύνης \((1-\alpha)100\%\) για την άγνωστη παράμετρο \(\theta\).
Πριν προχωρήσουμε θα δημιουργήσουμε μια γενική συνάρτηση που υλοποιεί το βήμα 3, δηλαδή παίρνει ως όρισμα ένα διάνυσμα τιμών από ένα τυχαίο δείγμα και υπολογίζει το αντίστοιχο διάστημα εμπιστοσύνης βασισμένο σε κανονική κατανομή. Η συνάρτηση επιστρέφει ένα διάνυσμα \((L, \hat{\theta}, U)\), όπου \(L,U\) είναι τα άκρα του διαστήματος και \(\hat{\theta}\) η σημειακή εκτίμηση του \(\theta\).
confinterval = function(x, a)
{
n=length(x)
thest=mean(x)
s=sd(x)
halflength=qt(1-a/2,n-1)*s/sqrt(n)
L=thest-halflength
U=thest+halflength
return(c(L,thest,U))
}Για παράδειγμα εφαρμογής της συνάρτησης, δημιουργούμε ένα δείγμα μεγέθους 100 από την ομοιόμορφη κατανομή \(U(0,1)\) και υπολογίζουμε ένα διάστημα εμπιστοσύνης 95% για τη μέση τιμή (που γνωρίζουμε ότι είναι ίση με 1/2):
x=runif(100, 0,1)
cint=confinterval(x, 0.05)
cint## [1] 0.4471484 0.5026930 0.5582376intervallength=cint[3]-cint[1]
intervallength## [1] 0.1110892Φυσικά με δείγμα μεγέθους 10000 έχουμε μεγαλύτερη ακρίβεια, καθώς το μήκος του διαστήματος εμπιστοσύνης μειώνεται στο 1/10 του προηγούμενου.
x=runif(10000, 0,1)
cint=confinterval(x, 0.05)
cint## [1] 0.4948497 0.5005044 0.5061592intervallength=cint[3]-cint[1]
intervallength## [1] 0.01130949Να εκτιμηθεί το ολοκλήρωμα \(\theta = \int_0^1 e^x dx\).
Εκφράζουμε το ολοκλήρωμα ως \(\theta = E(X)\), όπου \(X=e^U\) και \(U \sim U(0,1)\). Για την προσομοίωση μπορούμε να δημιουργήσουμε ένα δείγμα μεγέθους \(n\), \(u=(u_1, \ldots, u_n)\) από την κατανομή \(U(0,1)\) και θέτουμε \(x_j=e^{u_j}, j=1,\ldots,n\). Έπειτα υπολογίζουμε ένα διάστημα εμπιστοσύνης για το \(\theta\).
n=10000
u=runif(n, 0, 1)
x=exp(u)
cint=confinterval(x, 0.05)
cint## [1] 1.713757 1.723420 1.733083Η θεωρητική τιμή του ολοκληρώματος είναι προφανώς \(\theta= e-1\), και παρατηρούμε ότι η προσέγγιση που υπολογίσαμε είναι πολύ ακριβής.
estimate_accuracy=(cint[2]-(exp(1)-1))/(exp(1)-1)
estimate_accuracy## [1] 0.002990072Να εκτιμηθεί το ολοκλήρωμα \(\theta = \int_0^1 \int_0^1 e^{(x+y)^2} dy dx\).
Εκφράζουμε το ολοκλήρωμα ως \(\theta = E(X)\), όπου \(X=e^{(U+V)^2}\) και \(U,V\) ανεξάρτητες και ισόνομες τυχαίες μεταβλητές με κατανομή \(U(0,1)\). Για την προσομοίωση μπορούμε να δημιουργήσουμε δύο δείγματα μεγέθους \(n\), \(u=(u_1, \ldots, u_n)\) και \(v=(v_1, \ldots, v_n)\) από την κατανομή \(U(0,1)\) και θέτουμε \(x_j=e^{(u_j+v_j)^2}, j=1,\ldots,n\). Έπειτα υπολογίζουμε ένα διάστημα εμπιστοσύνης για το \(\theta\).
n=10000
u=runif(n, 0, 1)
v=runif(n, 0, 1)
x=exp((u+v)^2)
cint=confinterval(x, 0.05)
cint## [1] 4.756761 4.873235 4.989710Να εκτιμηθεί το ολοκλήρωμα \(\theta = \int_0^{\infty} \int_0^x e^{-(x+y)} dy dx\).
Eκφράζουμε το ολοκλήρωμα ως \(\theta = \int_0^{\infty} e^{-x} g(x) dx\),όπου \(g(x) = \int_0^x e^{-y} dy\). Για τα παραπάνω ολοκληρώματα παρατηρούμε :
Πρώτα για το \(\int_0^{\infty} e^{-x} g(x) dx\), εφαρμόζοντας το μετασχηματισμό \(z=\frac{1}{1+x}\) έχουμε \(x=\frac{1}{z}-1\), και \(dz = -z^2 dx\), επομένως \[ \theta = \int_0^{1} \frac{1}{z^2} e^{-\left (\frac{1}{z} -1 \right ) } g\left (\frac{1}{z} -1 \right) dz. \] Επιπλέον το ολοκλήρωμα \(g(x) = \int_0^x e^{-y} dy\) μπορεί να εκφραστεί ως \[ g(x) = \int_0^{\infty} 1(y \leq x) e^{-y} dy \] και εφαρμόζοντας πάλι το μετασχηματισμό \(s=\frac{1}{1+y}\), επομένως \(y=\frac{1}{s}-1\) και \(ds=s^2 dy\), παίρνουμε \[ g(x)=\int_0^{1} \frac{1}{s^2} 1\left (\frac{1}{s}-1 \leq x \right ) e^{-\left(\frac{1}{s}-1 \right)} ds. \]
Συνδυάζοντας τα παραπάνω γράφουμε το \(\theta\) ως \[ \theta = \int_0^{1} \int_0^1 \frac{1}{z^2} e^{-\left (\frac{1}{z} -1 \right ) } \frac{1}{s^2} 1\left (\frac{1}{s}-1 \leq \frac{1}{z}-1 \right ) e^{-\left(\frac{1}{s}-1 \right)} ds dz = E\left ( \frac{1}{U_1^2} e^{-\left (\frac{1}{U_1} -1 \right ) } \frac{1}{U_2}^2 e^{-\left (\frac{1}{U_2} -1 \right ) } 1\left(U_1 \leq U_2 \right ) \right ), \] όπου \(U_1, U_2\) α.ι.τ.μ. \(~U(0,1)\).
Για την προσομοίωση μπορούμε να δημιουργήσουμε δύο δείγματα μεγέθους \(n\), \(u=(u_1, \ldots, u_n)\) και \(v=(v_1, \ldots, v_n)\) από την κατανομή \(U(0,1)\) και θέτουμε \[ x_j = \frac{1}{u_j^2} e^{-\left (\frac{1}{u_j} -1 \right ) } \frac{1}{v_j^2} e^{-\left (\frac{1}{v_j} -1 \right ) } 1\left(u_j \leq v_j \right ) e^{-v_j}, j=1,\ldots,n. \] Έπειτα υπολογίζουμε ένα διάστημα εμπιστοσύνης για το \(\theta\).
n=10000
u=runif(n, 0, 1)
v=runif(n, 0, 1)
x=1/u^2 * exp(-(1/u -1)) * 1/v^2 * exp(-(1/v-1))*(u<=v)
cint=confinterval(x, 0.05)
cint## [1] 0.4839929 0.4982433 0.5124937Η τιμή του \(\theta\) μπορεί να υπολογιστεί αναλυτικά: \[ \theta = \int_0^{\infty} e^{-x} \int_0^x e^{-y} dy dx = \int_0^{\infty} e^{-x} (1-e^{-x}) dx = \frac{1}{2}, \] επομένως μπορούμε να ελέγξουμε την ακρίβεια της προσομοίωσης
theta=1/2
estimate_accuracy=(cint[2]-theta)/theta
estimate_accuracy## [1] -0.003513345Να εκτιμηθεί η συνδιακύμανση \(\theta=Cov(U, e^U)\) όπου \(U\sim U(0,1)\).
Η συνδιακύμανση αυτή μπορεί να εκφραστεί ως \(\theta = E(U e^U) - E(U) E(e^U) = E(U e^U) - \frac{1}{2} E(e^U) = E(Χ)\), όπου \(Χ= (U-\frac{1}{2})e^U\), και \(U\sim U(0,1)\). Επομένως μπορούμε να ακολουθήσουμε τη διαδικασία που εφαρμόσαμε και στα προηγούμενα παραδείγματα
n=10000
u=runif(n, 0,1)
x=(u-1/2)*exp(u)
cint=confinterval(x,0.05)
cint## [1] 0.1299220 0.1403054 0.1506888Η τιμή του \(\theta\) μπορεί να υπολογιστεί αναλυτικά \[ \theta = \int_0^1 u e^u du - \frac{1}{2} \int_0^1 e^u du = \frac{3-e}{2}, \] επομένως μπορούμε και εδώ να ελέγξουμε την ακρίβεια της προσομοίωσης
theta=(3-exp(1))/2
estimate_accuracy=(cint[2]-theta)/theta
estimate_accuracy## [1] -0.003930889Εδώ πρέπει να κάνουμε μια σημαντική παρατήρηση. Επειδή θεωρήσαμε την \(Ε(U)\) γνωστή και ίση με 1/2, δεν χρειάστηκε να την εκτιμήσουμε μέσω της προσομοίωσης. Αυτό ήταν μια κρίσιμη απλοποίηση, καθώς επέτρεψε να εκφραστεί το \(\theta\) ως μέση τιμή μιας μεταβλητής \(X\), και επομένως να καταλήξουμε σε αμερόληπτη εκτιμήτρια και να μπορούμε να υπολογίσουμε διαστήματα εμπιστοσύνης. Σε μια πιο γενική περίπτωση όπου ζητείται η συνδιακύμανση \(\theta=Cov(f(U), g(U)) = E(f(U)g(U)) - E(f(U)) Ε(g(U))\) και όλες οι μέσες τιμές πρέπει να εκτιμηθούν, το \(\theta\) δεν μπορεί να εκφραστεί ως \(Ε(Χ)\). Σε αυτή την περίπτωση μπορούμε να πάρουμε μια εκτιμήτρια \[ \hat{\theta} = \frac{1}{n} \sum_{j=1}^n f(u_j) g(u_j) - \frac{1}{n} \sum_{j=1}^n f(u_j) \frac{1}{n} \sum_{j=1}^n g(u_j) \] η οποία όμως είναι συνεπής αλλά όχι αμερόληπτη. Επίσης δεν μπορούμε να υπολογίσουμε ασυμπτωτικό διάστημα εμπιστοσύνης με βάση το κεντρικό οριακό θεώρημα. Για τέτοιες περιπτώσεις θα δούμε αργότερα ότι είναι χρήσιμες μέθοδοι επαναδειγματοληψίας όπως το bootstrap.
Να εκτιμηθεί η μέση τιμή \(\theta=E(N)\), όπου \(N=\min{n: \sum_{i=1}^n U_i >1}\), όπου \(U_1, U_2, \ldots\) α.ι.τ.μ. \(U(0,1)\).
Θα λύσουμε την άσκηση σε πιο γενικό πλαίσιο και θα εκτιμήσουμε την ποσότητα \(\theta(t)=E(N(t)\), όπου \(N(t)=\min{n: \sum_{i=1}^n U_i >t}\). Εδώ για να δημιουργήσουμε τυχαίο δείγμα από την κατανομή της \(N\) πρέπει να προσομοιώσουμε τον τρόπο κατασκευής της \(N\). Επομένως ορίζουμε μια συνάρτηση που δημιουργεί μια παρατήρηση της \(N(t\) προσομοιώνοντας διαδοχικά αθροίσματα \(\sum_{i=1}^n U_i\)
simreplication = function(t)
{
n=0
s=0
while (s<=t)
{
n=n+1
u=runif(1, 0, 1)
s=s+u
}
return(n)
}Τώρα δημιουργούμε ένα τυχαίο δείγμα από την \(N(t)\) π.χ. για \(t=1\).
n=10000
t=1
x=rep(0,n)
for (j in 1:n)
{
x[j]=simreplication(t)
}
cint=confinterval(x,0.05)
cint## [1] 2.696017 2.713200 2.730383Από τη θεωρία των ανανεωτικών στοχαστικών διαδικασιών προκύπτει ότι για \(t \leq 1\) η τιμή του \(\theta(t)\) μπορεί να υπολογιστεί αναλυτικά ως \(\theta(t)=e^t\), επομένως για το προηγούμενο πείραμα η ακρίβεια μπορεί και πάλι να εκτιμηθεί:
theta=exp(t)
estimate_accuracy=(cint[2]-theta)/theta
estimate_accuracy## [1] -0.0018695