Σήμερα είδαμε τα βασικά της μεθόδου Monte Carlo και της ανάλυσής της στην απλή περίπτωση που επιχειρούμε να εκτιμήσουμε το εμβαδό ενός χωρίου $\Omega \subseteq [0, 1]^2$. Η μέθοδος συνίσταται στο να "πετάμε" τυχαία σημεία (με την ομοιόμορφη κατανομή) στο $[0, 1]^2$ και να μετράμε τι ποσοστό από αυτά πέφτει μέσα στο $\Omega$. Στην περίπτωση π.χ. που το $\Omega$ είναι ένας κυκλικός δίσκος με κέντρο το $(1/2, 1/2)$ και ακτίνα $1/\sqrt{10}$, το σύνολο δηλ. $$\left\{(x,y):\ (x-0.5)^2+(y-0.5)^2 \le \frac{1}{10}\right\}$$ φαίνεται παρακάτω το πρόγραμμα υπολογισμού του εμβαδού του.
import math
import numpy as np
def inout(x, y): # Επιστρέφει True αν και μόνο αν (x, y) είναι στοιχείο του Ω
return (x-0.5)**2+(y-0.5)**2 <= 0.1
area = math.pi/10 # Η πραγματική απάντηση, για λόγους σύγκριση με την προσέγγιση
N = 10000 # πόσα σημεία ρίχνουμε
k = 0 # πόσα από αυτά πέφτουν μέσα
for i in range(N):
x = np.random.uniform() # τυχαία πρώτη συντεταγμένη
y = np.random.uniform() # τυχαία δεύτερη συντεταγμένη
if inout(x, y): # αν πέσει μέσα αυξάνουμε το k κατά 1
k += 1
estimate = k/N # Η εκτίμησή μας για το εμβαδό είναι το ποσοστό των σημείων που έπεσαν μέσα.
print(f"Εκτίμηση Monte Carlo = {estimate}, πραγματική τιμή = {area}")
Ο παραπάνω αλγόριθμος παράγει τυχαία αποτελέσματα. Θέλουμε να μην αποκλίνει πολύ από την πραγματική τιμή, και με όσο γίνεται πιο μεγάλη πιθανότητα. Αυτό φυσικά συμβαίνει όσο πιο μεγάλο πάρουμε το $N$ (πλήθος σημείων) αλλά πόσο μεγάλο πρέπει να είναι αυτό;
Ξεκινάμε πάντα θέτοντας τις προδιαγραφές μας. Αν ονομάσουμε $X$ το αποτέλεσμα του αλγορίθμου τότε η πραγματική τιμή (η σωστή απάντηση) είναι η μέση τιμή $\mathbb{E}(X)$. Αν πάρουμε λοιπόν $\epsilon$ και $\delta$ να είναι δύο μικροί θετικοί αριθμοί ζητάμε πόσο μεγάλο πρέπει να είναι το $N$ ώστε να ισχύει η ανισότητα $$\mathbb{P}\left(|X-\mathbb{E}X| > \epsilon\right) < \delta.\ \ \ \ \ (1)$$ Το $\epsilon$ είναι η παράμετρος ακρίβειας (τι απόκλιση είμαστε διατεθειμένοι να δεχτούμε από τη σωστή τιμή για να μην το θεωρήσουμε αποτυχία) και το $\delta$ είναι η παράμετρος που μας λέει τι πιθανότητα αποτυχίας είμαστε διατεθειμένοι να αποδεχτούμε. Για παράδειγμα, για το πρόβλημα του κύκλου παραπάνω, θα μπορούσαμε να επιλέξουμε $\epsilon=0.001$ και $\delta=0.01$. Αυτό σημαίνει ότι είμαστε διατεθειμένοι να αποδεχτούμε το 1% των φορών ($1-\delta$) να πέφτουμε πάνω από 0.001 μακριά από τη σωστή τιμή.
Αφού επιλέξουμε τις παραμέτρους $\epsilon$ και $\delta$ προσπαθούμε να βρούμε ένα $N_0$, όσο το δυνατόν πιο μικρό, ώστε αν πάρουμε $N \ge N_0$ να έχουμε εξασφαλίσει ότι θα ισχύει η ανισότητα (1) παραπάνω. Αυτό το επιτυγχάνουμε, σε πρώτη φάση, με τη χρήση της ανισότητας Chebyshev.
Για να δείτε πώς γίνεται αυτό πηγαίνετε στις ενότητες 10.7 έως 10.9 (10η εβδομάδα) στο μάθημα Θεωρίας Πιθανοτήτων 2021-22 (θα πρέπει να κάνετε login με τα Πανεπιστημιακά σας στοιχεία και enroll στο μάθημα για να το δείτε).
Στην ανάλυση της μεθόδου Monte Carlo μπορεί κανείς αντί της ανισότητας Chebyshev να χρησιμοποιήσει άλλες καλύτερες ανισότητες που εκμεταλλεύονται το γεγονός ότι η τυχαία μεταβλητή που αναλύουμε είναι άθροισμα πολλών ανεξαρτήτων τυχαίων μεταβλητών. Αυτή η δομή οδηγεί συχνά σε εκθετικά μικρά φράγματα απόκλισης από τη μέση τιμή (αυτό το φαινόμενο ονομάζεται συγκέντρωση). Μια πολύ χρήσιμη από αυτές τις ανισότητες είναι η παρακάτω ανισότητα του Chernoff.
Ανισότητα Chernoff: Αν $X_1, \ldots, X_n$ είναι ανεξάρτητες ΤΜ που παίρνουν μόνο τις τιμές 0 ή 1 και $p_i = \mathbb{P}(X_i=1)$ και $X = \sum_{i=1}^n X_i$, τότε αν $\mu = \mathbb{Ε}(X) = \sum_{i=1}^n p_i$ και $0 < \delta < 1$ ισχύει $$\mathbb{P}(|X-\mu| \ge \delta \mu) \le 2 e^{-\frac{\delta^3}{3}\mu}.$$
Αυτή είναι μια εξαιρετικά ευέλικτη και χρήσιμη ανισότητα (ειδικά για συνδυαστικά προβλήματα) και χρησιμοποιώντας την στο πρόβλημα Monte Carlo παραπάνω μπορεί κανείς να πάρει πολύ καλύτερες (μικρότερες) τιμές για την ποσότητα $N_0$ (ελάχιστος αριθμός σημείων ώστε να είναι εξεασφαλισμένη η επιτυχία με πιθανότητα δεδομένη όπως περιγράφεται παραπάνω).
Η μέθοδος Monte Carlo είναι μια πολύ ευέλικτη μέθοδος και πολλές φορές χρησιμοποιείται για να προσομοιώσει τη λειτουργία διαφόρων συστημάτων.
Στο παρακάτω παράδειγμα μια βιομηχανία παράγει 4 αντικείμενα: τα A, B, C και το κουτί D, μέσα στο οποίο πρέπει να μπαίνουν τα A, B, C, με αυτή τη σειρά, με κενό 1mm ανάμεσά τους και ανάμεσα σε αυτά και τα τοιχάωματα του κουτιού. Αν το κουτί έχει μήκος 100mm, και τα A, B, C έχουν μήκη 30mm, 20mm και 46mm αντίστοιχα, τότε το κουτί D χωράει ακριβώς τα κουτιά A, B, C μαζί με τα διάκενα του 1mm που έχουν προβλεφθεί.
Όμως τίποτα δεν είναι τέλειο και τα A, B, C, D έχουν στην πραγματικότητα μήκη που είναι περίπου τα προβλεφθέντα αλλά όχι ακριβώς. Μπορούμε να το προσομοιώσουμε αυτό λέγοντας ότι κάθε μια από αυτές τις ποσότητες είναι ίση με την ονομαστική της τιμή $\pm$ ένα τυχαίο (και, ελπίζει κανείς, μικρό) μήκος. Ας πούμε λοιπόν ότι τα μήκη αυτά (τα τυχαία) ακολουθούν κανονική κατανομή με μέση τιμή 0 και τυπική απόκλιση $\sigma$, ανεξάρτητα το ένα από το άλλο.
Ποια η πιθανότητα τότε να μη χωράνε τα τρία κουτιά στο κουτί D, ακόμη κι αν θυσιάσουμε τα προβλεφθέντα διάκενα; Ποια η πιθανότητα δηλ. να ισχύει $A+B+C > D$;
Την πιθανότητα αυτή την προσεγγίζουμε εκτελώντας πολλές φορές το αντίστοιχο πείραμα και μετρώντας τι ποσοστό φορών θα ισχύσει αυτό το κακό αποτέλεσμα. Αυτό φαίνεται στο πρόγραμμα παρακάτω. Είναι φανερή η ευελιξία της μεθόδου.
import numpy as np
N = 10000 # πόσες φορές εκτελούμε το πείραμα
sigma = 2 # Η τυπική απόκλιση (σε mm) των τυχαίων μηκών
# που εμφανίζονται ως αποκλίσεις από τις ονομαστικές τιμές.
bad = 0 # Πόσες φορές συμβαίνει το κακό ενδεχόμενο
for i in range(N):
A = 30 + np.random.normal(0, sigma) # Στις επόμενες 4 γραμμές υπολογίζουμε τις A, B, C, D
B = 20 + np.random.normal(0, sigma)
C = 46 + np.random.normal(0, sigma)
D = 100 + np.random.normal(0, sigma)
if A+B+C > D: # Αν συμβεί το κακό ενδεχόμενο αυξάνουμε το μετρητή
bad += 1
estimate = bad/N # Η εκτίμησή μας για την πιθανότητα του κακού ενδεχομένου
print(f"Bad probability estimate: {estimate}")