Ένα τυπικό πρόβλημα που προσπαθούμε να λύσουμε με τη μέθοδο Monte Carlo είναι να εκτιμήσουμε μια μέση τιμή $\mathbb{E}f(X)$ για μια ΤΜ $X$, όπου $f(x)$ είναι μια συνάρτηση την οποία μπορούμε να υπολογίσουμε αν γνωρίζουμε το $x$.
Στο πρόβλημα με το εμβαδό του κύκλου μέσα στο $[0, 1]^2$ που λύσαμε την προηγούμενη φορά η $X$ είναι μια (διανυσματική) ΤΜ ομοιόμορφα κατανεμημένη στο $[0, 1]^2$ και $f:[0, 1]^2 \to \mathbb{R}$ είναι η χαρακτηριστική συνάρτηση του δίσκου με κέντρο το $(\frac12, \frac12)$ και ακτίνα $\frac12$ (παίρνει δηλ. την τιμή 1 μέσα στο δίσκο και 0 εκτός). Η μέση τιμή $\mathbb{E}f(X)$ είναι το εμβαδό που θέλουμε να υπολογίσουμε.
Η μέθοδος, όπως την έχουμε δει μέχρι στιγμής, συνίσταται ουσιαστικά στο να εκτιμούμε τη μέση τιμή $\mathbb{E}f(X)$ από το δειγματικό μέσο $$\frac{1}{N} \sum_{i=1}^N f(X_i)$$ όπου $X_i$, $i=1, 2, \ldots, N$, είναι ανεξάρτητα δείγματα της $X$.
Από το νόμο των μεγάλων αριθμών στη θεωρία πιθανοτήτων είναι σαφές ότι με κάποιες ασθενείς προϋποθέσεις για την ΤΜ $f(X)$ το όριο, για $N \to \infty$, του δειγματικού μέσου είναι, με πιθανότητα 1 (σχεδόν σίγουρα όπως λέμε) η ποσότητα $\mathbb{E}f(X)$. Δεν αρκεί όμως αυτό. Στην πράξη απαιτείται γρήγορη σύγκλιση, να έχουμε δηλ., με μεγάλη πιθανότητα, ότι για μικρό $N$ ο δειγματικός μέσος είναι κοντά στην πραγματική μέση τιμή.
Αυτό δε συμβαίνει αν η συνάρτηση $f$ είναι μεγάλη σε πολύ μικρό μέρος των τιμών της $X$. Φανταστείτε για παράδειγμα ότι θέλουμε να εκτιμήσουμε το εμβαδό ενός χωρίου $\Omega \subseteq [0, 1]^2$, όπου όμως το εμβαδό του $\Omega$ είναι πολύ μικρό, π.χ. της τάξης του $10^{-6}$. Κατά μέσο όρο, τότε, θα χρειαστεί να πάρουμε το $N$ να είναι τουλάχιστον $10^6$ για να δούμε έστω και ένα δείγμα που να πέφτει μέσα στο $\Omega$ (για να είμαστε πιο ακριβείς, η μέση τιμή του $i$ που είναι η πρώτη φορά που θα πέσουμε μέσα στο $\Omega$ είναι ίση με $|\Omega|^{-1} = 10^6$ (η τυχαία αυτή μεταβλητή $i$ ακολουθεί τη γεωμετρική κατανομή με παράμετρο $p=10^{-6}$).
Η μέθοδος του importance sampling προσπαθεί να διορθώσει αυτό το πρόβλημα (πολύ μεγάλα $N$) παίρνοντας δείγματα από μια άλλη κατανομή που είναι μεγάλη εκεί όπου η συνάρτηση $f$ είναι επίσης μεγάλη. Για παράδειγμα, στο πρόβλημα της εκτίμησης του εμβαδού του χωρίου $\Omega$ που αναφέρουμε παραπάνω, αν κάπως έχουμε την πληροφορία ότι το $\Omega$ δεν περιέχεται απλά στο $[0, 1]^2$ αλλά σε ένα πολύ μικρότερο σύνολο, π.χ. στο ορθογώνιο $R = [0, 10^{-3}] \times [0, 1]$, τότε μπορούμε να πάρουμε τα τυχαία μας δείγματα μέσα στο $R$ όπου το $\Omega$ καταλαμβάνει το $10^{-3}$ του χώρου (και όχι το $10^{-6}$ όπως είναι ποσοστό του $\Omega$ που καταλαμβάνει) και έτσι χρειαζόμαστε κατά μέσο όρο μόνο $10^3$ δείγματα για να δούμε ένα δείγμα του $\Omega$ και όχι $10^6$).
Βασικό συστατικό στη μέθοδο importance sampling είναι λοιπόν να ξέρουμε εκ των προτέρων μια άλλη κατανομή, πέραν αυτής του $X$, η οποία να δίνει πιο συχνά μεγάλες τιμές για την $f(X)$.
Μπορείτε να διαβάσετε αναλυτικά για τη μέθοδο αυτή στη σελίδα αυτή, όπου θα δείτε και το παράδειγμα που κάναμε στο μάθημα και το οποίο υλοποιήσαμε στο μάθημα όπως φαίνεται παρακάτω (στη σελίδα θα δείτε μια λίγο διαφορετική υλοποίηση του ίδιου παραδείγματος).
import math
import numpy as np
N = 100000 # ο αριθμός των δειγμάτων και για τις δύο μεθόδους (απλό Monte Carlo και με importance sampling)
def normaldens(t): # Η πυκνότητα της τυπικής κανονικής κατανομής
return 1/np.sqrt(2*math.pi) * np.exp(-t*t/2)
def g(t): # Θέλουμε να υπολογίσουμε τη μέση τιμή Eg(X)
return 10*np.exp(-5*(t-3)**4)
def gg(t): # Στη μέθοδο importance sampling υπολοιγίζουμε τη μέση τιμή Ε gg(Y) όπου Y ακολουθεί
# την πυκνότητα N(3, 1)
return g(t) * normaldens(t) / normaldens(t-3)
x = np.random.normal(0, 1, N) # απλή Monte Carlo
v = g(x)
approx1 = np.mean(v)
var1 = np.var(v) # δειγματική διασπορά των g(X_i)
y = np.random.normal(3, 1, N) # με importance sampling
w = gg(y)
approx2 = np.mean(w)
var2 = np.var(w) # δειγματική διασπορά των gg(Y_i)
print(f"Προσέγγιση απλή: {approx1},\n με importance sampling: {approx2}")
print(f"απλή διασπορά: {var1},\n με importance sampling: {var2}") # όσο μικρότερη τόσο καλύτερα
Αναφερόμαστε στις διαλέξεις Ryan O'Donnell, Probability and Computing, Lecture 14.
Το πρόγραμμα που φαίνεται παρακάτω είναι ο πολύ εύκολος υπολογισμός της 100ής δύναμης του πίνακα μετάβασης $K$, και όπου παρατηρούμε ότι οι γραμμές του είναι πρακτικά ίδιες.
import numpy as np
K = np.array([[0.4, 0.6, 0], [0.1, 0.6, 0.3], [0.5, 0, 0.5]])
t=100
Kt = np.linalg.matrix_power(K, t)
print(Kt)