7.4 Αλγόριθμοι Monte Carlo Μαρκοβιανών Αλυσίδων

Έστω μια ακολουθία τ.μ. της μορφής:

$$x^{(0)},x^{(1)},x^{(2)},\mathrm{\dots },x^{(t)}.$$

H παραπάνω σειρά («αλυσίδα») ονομάζεται Μαρκοβιανή αν ισχύει η Μαρκοβιανή ιδιότητα, δηλαδή αν η τιμή της σειράς σε μια χρονική στιγμή, έστω $t+1$, εξαρτάται μόνο από τη τιμή της σειράς στην αμέσως προηγούμενη στιγμή, $t$. Δηλαδή αν η πιθανότητα η τ.μ. να λάβει μια τιμή στην $t+1$ εξαρτάται μόνο από την τιμής της την αμέσως προηγούμενη τιμή:

$$P(x^{(t+1)}=x_{t+1}|x^{(t)}=x_t,\mathrm{\dots },x^{(0)}=x_0)=P(x^{(t+1)}=x_{t+1}|x^{(t)}=x_t)$$

(7.12)

όπου $x_{t+1},x_t,\mathrm{\dots },x_0$ είναι οι τιμές της τ.μ. $X$ στις χρονικές στιγμές $t+1,t,\mathrm{\dots }$ και $0$ αντίστοιχα.

Το παράδειγμα αυτό είναι αρκετά γενικό και μπορεί να περιγράφει ένα πολύ μεγάλο εύρος προβλημάτων. Ως εφαρμογή θεωρήστε ότι $x_n$ είναι τ.μ. που ακολουθούν την κατανομή Bernoulli (π.χ. προέρχονται από ρίψη νομισμάτων) με πιθανές τιμές $x_n=1$ ή $x_n=-1$ και:

$$P(x_n=1)=1-P(x_n=-1)=P.$$

Έστω τώρα η τ.μ. $S$, οι τιμές της οποίας ορίζονται ως:

$$S^{(t)}=x_1+x_2+\mathrm{\cdots }+x_t,S^{(0)}=0$$

δηλαδή η οποία περιγράφει τη διαφορά μεταξύ φορών «γράμματα» και «κορώνας» σε $n$ ρίψεις νομισμάτων. Εύκολα μπορούμε να δείξουμε ότι η $S^{(t)},t=0,1,2,\mathrm{\dots }$, είναι μια Μαρκοβιανή αλυσίδα.

Στο κεφάλαιο αυτό παρουσιάζουμε το πιο γνωστό αλγόριθμο Monte Carlo βασισμένο σε Μαρκοβιανές αλυσίδες, αλλά και γενικότερα το πιο γνωστό αλγόριθμο Monte Carlo. Θεωρούμε το παρακάτω πρόβλημα, που συζητήσαμε και σε προηγούμενο υπο-κεφάλαιο. Έστω ότι θέλουμε να υπολογίσουμε ιδιότητες μέσω προσομοίωσης για ένα στοχαστικό σύστημα το οποίο ακολουθεί την κατανομή Boltzmann:

$$\pi (x)=\frac{\exp \mathrm{}[-h(𝒙)]}{\iiint \exp \mathrm{}[-h(𝒙)]d𝒙}=\frac{1}{Z}\exp \mathrm{}[-h(x)].$$

(7.16)

Στην παραπάνω σχέση η $h(x)$ είναι μια συνάρτηση της ($n$-διάστατης στη γενική μορφή) τυχαίας μεταβλητής $𝒙$ και $Z$ είναι μια σταθερά κανονικοποίησης ($n$-πολλαπλό ολοκλήρωμα).

Η σταθερά $Z$ είναι συνήθως άγνωστη, καθώς ο υπολογισμός της απαιτεί τον υπολογισμό του $n$-διάστατου ολοκληρώματος, πρόβλημα που είναι εξίσου δύσκολο, αν όχι πιο δύσκολο, από το αρχικό πρόβλημα, δηλαδή την προσομοίωση της $\pi (𝒙)$. Το πρόβλημα αυτό είναι ιδιαίτερα σημαντικό καθώς όλα τα μοριακά συστήματα σε ισορροπία ακολουθούν την κατανομή Boltzmann. Σε αυτή την ιδιότητα οφείλεται κυρίως το πολύ μεγάλο εύρος εφαρμογών του αλγορίθμου Metropolis Monte Carlo.

Το παραπάνω πρόβλημα επιλύεται με τη χρήση Μαρκοβιανών αλυσίδων και τον αλγόριθμο Metropolis.

Ο αλγόριθμος Metropolis αναπτύχθηκε από τους N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller and E. Teller (MRRTT) το 1953 και είναι ο πιο γνωστός αλγόριθμος Monte Carlo.

Θεωρείστε μια $n$-διάστατη τυχαία μεταβλητή, στη χρονική στιγμή $t,𝒙(t)=\{x_1(t),x_2(t),\mathrm{\dots },x_n(t)\}$ η οποία ακολουθεί την κατανομή $\pi (𝒙)$, παράδειγμα την κατανομή Boltzmann. Ο αλγόριθμος Metropolis αποτελείται από τα παρακάτω βήματα:

• •

  Βήμα 1: Αρχικά δημιουργούμε μια πρόταση ($n$-διάστατη τυχαία μεταβλητή) $𝒚$ από ένα κανόνα μετάβασης (προτεινόμενη κατανομή) $T(𝒙(t),𝒚)$.

• •

  Βήμα 2: Κατόπιν υπολογίζουμε την ποσότητα $r(𝒙(t),𝒚)$ όπου:

  $$r({𝒙}^{(t)},𝒚)=\min \{1,\frac{\pi (𝒚)}{\pi (𝒙)}\}.$$

  (7.17)

• •

  Βήμα 3: Δημιουργία τυχαίου αριθμού από ομοιόμορφη κατανομή στο διάστημα $(0,1)$, δηλαδή $u\sim U(0,1)$. Η νέα τιμή της Μαρκοβιανής αλυσίδας $𝒙(t+1)$ είναι:

  $${𝒙}^{(t+1)}=\{\begin{array}{cc}𝒚,\text{αν}\hspace{1em}u\le r({𝒙}^{(t)},𝒚)\hfill & \\ {𝒙}^{(t)},\text{αν}\hspace{1em}u>r({𝒙}^{(t)},𝒚)\hfill & \end{array}$$

  (7.18)

Όπως βλέπουμε από τον αλγόριθμο αυτό η νέα τιμή της Τ.Μ., σε χρόνο $t+1$, είναι είτε η προηγούμενη (δηλαδή αυτή σε χρόνο $t$) είτε η νέα, η οποία προτείνεται με βάση κάποιο κανόνα μετάβασης.

Τέλος από την Μαρκοβιανή αλυσίδα της μεταβλητής $𝒙$, μπορούμε να υπολογίσουμε συναρτήσεις της $𝒙$ οι οποίες σχετίζονται με τις ιδιότητες που μας ενδιαφέρουν. Ως παράδειγμα μπορούμε να βρούμε την μέση τιμή συναρτήσεων της μεταβλητής $f(𝒙)$, δηλαδή να υπολογίσουμε $n$-διάστατα ολοκληρώματα της μορφής:

$$E\{f(𝒙)\}=\int f(𝒙)d𝒙.$$

Ο παραπάνω αλγόριθμος υποθέτει ότι ο κανόνας μετάβασης είναι συμμετρικός, δηλαδή $T(𝒙(t),𝒚)=T(𝒚,𝒙(t))$. Λίγα χρόνια αργότερα, το 1970, ο Hastings πρότεινε μια γενίκευση του παραπάνω αλγορίθμου για μη συμμετρικούς κανόνες μετάβασης.

Η μόνη διαφορά το γενικευμένου αλγορίθμου Metropolis – Hastings σε σχέση με τον αρχικό αλγόριθμο Metropolis είναι στον υπολογισμό της ποσότητας $r(𝒙(t),𝒚)$. Πιο συγκεκριμένα, ο αλγόριθμος Metropolis – Hastings (1970) αποτελείται από τα παρακάτω βήματα:

• •

  Βήμα 1: Αρχικά δημιουργούμε μια πρόταση (n-διάστατη τυχαία μεταβλητή) $𝒚$ από ένα κανόνα μετάβασης (προτεινόμενη κατανομή) $T(𝒙(t),𝒚)$.

• •

  Βήμα 2: Κατόπιν υπολογίζουμε την ποσότητα $r(𝒙(t),𝒚)$ όπου:

  $$r({𝒙}^{(t)},𝒚)=\min \{1,\frac{\pi (𝒚)T(𝒚,{𝒙}^{(t)})}{\pi (𝒙)T({𝒙}^{(t)},𝒚)}\}$$

  (7.19)

• •

  Βήμα 3: Δημιουργία τυχαίου αριθμού από ομοιόμορφη κατανομή στο διάστημα $(0,1)$, δηλαδή $u\sim U(0,1)$. Η νέα τιμή της Μαρκοβιανής αλυσίδας ${𝒙}^{(t+1)}$ είναι:

  $${𝒙}^{(t+1)}=\{\begin{array}{cc}𝒚,\text{αν}\hspace{1em}u\le r({𝒙}^{(t)},𝒚)\hfill & \\ {𝒙}^{(t)},\text{αν}\hspace{1em}u>r({𝒙}^{(t)},𝒚)\hfill & \end{array}$$

  (7.20)

Ο αλγόριθμος Metropolis – Hastings (η ο αλγόριθμος Metropolis) χρησιμοποιείται στις μέρες μας κατά κόρον σε ένα τεράστιο εύρος επιστημονικών πεδίων, από τη στατιστική φυσική, τη βιολογία, τη χημεία, και την οικονομία. Στο τέλος του κεφαλαίου υπάρχει πρότυπο ανάπτυξης γενικού αλγορίθμου Metropolis – Hasting σε γλώσσα MATLAB.