Algoritmo Monte Carlo

[tetravalenza]

Ciao,

Il libro "Computer Graphics With OpenGL" (Hearn, Baker, Pearson, Pag. 768) dà una breve spiegazione del famoso algoritmo Monte Carlo, in particolare espone la seguente relazione
\[
\int_{a}^{b} f(x)dx\approx h(b-a)\cdot \frac{n_{count}}{n}\\
h =y_{max} - y_{min}\\
x = a + r1(b-a)\\
y = y_{min}+r2\cdot h

\]

Dove h è l'altezza del rettangolo che contiene la curva \(f(x)\), \(n_{count}\) è il numero di punti casuali che cadono tra la \(f(x)\) e l'asse \(x\), n è il numero di intervalli, \(r1, r2\) sono i numeri casuali. Supponiamo io voglia calcolare l'area della curva \(sin(x)\) nell'intervallo da 0 a pi greco. Come verifico che il punto di coordinate (x, y) sta nella regione di piano tra la curva e, per esempio, \(y_{min}\)?

[apatriarca]

Per rispondere alla tua specifica domanda, il tuo punto casuale \((x_i,y_i)\) starà nella regione compresa tra \(y_{min}\) e la curva \(f(x)\) se
\[y_{min} \leq y_i \leq f(x_i).\]

Ho tuttavia alcune osservazioni sul metodo Monte Carlo descritto. L'algoritmo Monte Carlo converge con complessità \(O(\sqrt{n})\) indipendentemente dalla dimensione del dominio di integrazione. Quando la dimensione è piccola (come in questo caso) altri metodi sono MOLTO più efficienti (pure qualcosa di molto semplice come la [url=Regola del Rettangolo]Regola del Rettangolo[/url] converge più velocemente nel tuo esempio), ma diventa utile quando la dimensione del dominio di integrazione cresce. La versione descritta è una modifica dell'algoritmo base che non è particolarmente utile in generale e non fornisce alcun vantaggio rispetto alla sua versione base. Non mi è mai stato quindi chiaro lo scopo di trattarlo.

Quando cerchiamo di calcolare un integrale, lo convertiamo di solito in una somma del tipo
\[ \int_{\Omega} f(\bar{x})\,d\bar{x} \approx \sum_{i=1}^N \omega_i\,f(\bar{x}_i), \]
dove i coefficienti \(\omega_i\) e le posizioni \(\bar{x}_i\) dipendono dall'algoritmo. Nel caso dell'algoritmo descritto sopra dei rettangoli, per un integrale in una singola variable, abbiamo ad esempio \(w_i = \frac{b - a}{N}\) e gli \(x_i\) presi uniformemente nell'intervallo.

Di solito si cerca di scegliere questi punti in modo "ottimale" o che almeno ricoprano in modo uniforme tutto il dominio. Ma con l'aumentare delle dimensioni del dominio il numero dei punti necessari a ricoprire il dominio cresce in modo esponenziale. L'idea del metodo di integrazione Monte Carlo è quella di prendere i punti in modo casuale. L'errore commesso sarà diverso ad ogni esecuzione dell'algoritmo (può anche essere molto grande), ma in media il risultato sarà corretto. L'errore scenderà inoltre lentamente con l'aumentare dei punti.

La formula normalmente usata è quindi:
\[ \int_{\Omega} f(\bar{x})\,d\bar{x} \approx \frac{1}{N} \sum_{i=1}^N \frac{f(\bar{x}_i)}{p(\bar{x}_i)} \]
dove \(p(x_i)\) è la probabilità di scegliere quel punto. Prendendo i punti in modo uniforme nel dominio avrai quindi
\[ \int_{\Omega} f(\bar{x})\,d\bar{x} \approx \frac{V_\Omega}{N} \sum_{i=1}^N f(\bar{x}_i) \]
dove
\[ V_\Omega = \int_\Omega d\bar{x} \]
è il volume/area del dominio di integrazione.

Tornando al tuo esempio, per una singola variabile abbiamo:
\[ \int_a^b f(x)\,dx \approx \frac{b - a}{N} \sum_{i=1}^N f\bigl(a + \xi_i\,(b-a)\bigr) \]
dove gli \(\xi_i\) sono numeri casuali scelti uniformemente in \([0, 1]\). Nel caso particolare di \(f(x) = \sin(x)\) e \([a, b] = [0, \pi],\) è sufficiente calcolare \(N\) valori casuali in \([0, \pi]\), sommare il seno di questi valori e moltiplicare per \(\pi/N\).

Il tuo libro calcola l'integrale in modo diverso. In effetti non sta integrando la tua funzione direttamente, ma sta integrando
\[ H(x, y) = \begin{cases}
0 & y > f(x) \\
1 & y \leq f(x).
\end{cases} \]
nel dominio \( [a, b] \times [y_{min}, y_{max}]. \)

[tetravalenza]

Grazie ancora per la spiegazione.
Ho provato a implementare sia quello del libro che quello dei rettangoli e a confrontare i risultati: sono assai diversi.
Con n = 10000 il Monte Carlo del libro restituisce un'area di 0.8108 mentre il metodo dei rettangoli restituisce 2.0096423977065747 che è corretto.
Con il metodo del libro ho calcolato il rapporto tra il numero dei punti che cadono nell'area alla curva e il numero di test effettuati, ma devo moltiplicarlo per qualche costante?

[apatriarca]

Come hai generato i campioni? Hai moltiplicato per \( (y_{max} - y_{min})\,(b - a) \)?

[tetravalenza]

Come hai generato i campioni? Hai moltiplicato per \( (y_{max} - y_{min})\,(b - a) \)?

No, non avevo moltiplicato il rapporto per i due fattori che hai detto. Ora viene giusto. Grazie mille!