diff --git a/semester3/numcs/numcs-summary.pdf b/semester3/numcs/numcs-summary.pdf index ec20943..077ae66 100644 Binary files a/semester3/numcs/numcs-summary.pdf and b/semester3/numcs/numcs-summary.pdf differ diff --git a/semester3/numcs/numcs-summary.tex b/semester3/numcs/numcs-summary.tex index 7a65b33..3c7be3b 100644 --- a/semester3/numcs/numcs-summary.tex +++ b/semester3/numcs/numcs-summary.tex @@ -74,6 +74,11 @@ The numbering should match the script's numbering exactly (apart from the cases Many of the figures in this summary were taken directly from the Script or Lecture notes created by Professor Vasile Gradinaru. +We have also taken some explanations and code examples from the slides of our TA, Nils Müller, whose slides can be found +\color{MidnightBlue}\fbox{\href{https://n.ethz.ch/~muellerni/courses/numcs25.php}{here}}\color{black} +% TODO: Update this when n.ethz is taken offline completely + + To add to the one quote regarding Python and type annotation: This is objectively wrong and a really hot take. Yes, this applies for small projects, but libraries \textit{DO} need type annotation, as you can't possibly read the entire library's code to use it. The reason this quote was even included here is that his coding style is really awful (yes, there were semicolons in his Python-code sometimes) diff --git a/semester3/numcs/parts/02_quadrature/04_in-rd.tex b/semester3/numcs/parts/02_quadrature/04_in-rd.tex index fbac474..4810a2c 100644 --- a/semester3/numcs/parts/02_quadrature/04_in-rd.tex +++ b/semester3/numcs/parts/02_quadrature/04_in-rd.tex @@ -1,4 +1,4 @@ -\newsection +\newsectionNoPB \subsection{Quadratur in $\R^d$ und dünne Gitter} Eine einfache Option wäre natürlich, zwei eindimensionale Quadraturformeln aneinander zu hängen. Für zweidimensionale Funktionen sieht dies so aus: diff --git a/semester3/numcs/parts/02_quadrature/05_monte-carlo.tex b/semester3/numcs/parts/02_quadrature/05_monte-carlo.tex index 57a696f..e4700af 100644 --- a/semester3/numcs/parts/02_quadrature/05_monte-carlo.tex +++ b/semester3/numcs/parts/02_quadrature/05_monte-carlo.tex @@ -3,10 +3,23 @@ \subsection{Monte-Carlo Quadratur} Bei der Monte-Carlo Quadratur wird, wie bei anderen Monte-Carlo-Algorithmen der Zufall genutzt +\shade{Orange}{Grundrezept} Wir nehmen $N$ zahlen $t_i$ die zufällig aus der uniformen Verteilung auf $[0, 1]$ gewählt werden. +\rmvspace +\begin{align*} + I = \int_{0}^{1} z(t) \dx t \approx \frac{1}{N} \sum_{i = 1}^{N} z(t_i) +\end{align*} + +\drmvspace Auf einem anderen Intervall $[a, b]$ haben wir dann für $s_i = a + (b - a) \cdot t_i$ +\rmvspace +\begin{align*} + I = \int_{a}^{b} z(s) \dx s \approx |b - a| \frac{1}{N} \sum_{i = 1}^{z(s_i)} =: I_N +\end{align*} + \inlineremark Die Konvergenz ist sehr langsam ($\sqrt{N}$), aber nicht abhängig von der Dimension oder Glattheit. Zudem kann das Ergebnis falsch sein, da es probabilistisch ist. -Jede Monte-Carlo-Methode benötigt folgendes mit $X = [I_N - \tsigma_N, I_N + \tsigma_N]$: +Jede Monte-Carlo-Methode benötigt folgendes mit $X = [I_N - \tsigma_N, I_N + \tsigma_N]$ und $X$ enthält den wahren Wert $\int_{[0, 1]^d} z(t) \dx t$ in ungefähr $68.3\%$ der Fälle +für $N$ $t_i$ uniform Verteilt in $[0, 1]^d$: \rmvspace \begin{multicols}{2} \begin{itemize}[noitemsep] @@ -16,28 +29,33 @@ Jede Monte-Carlo-Methode benötigt folgendes mit $X = [I_N - \tsigma_N, I_N + \t \item Darstellung des Ergebnis, hier $\Pr[I \in X] = 0.683$ \end{itemize} \end{multicols} -Mit $\displaystyle I_N = \int_{0}^{1} z(t) \dx t = \frac{1}{N} \sum_{i = 1}^{N} z(t_i)$, wobei $t_i$ Zufallszahlen sind und -\rmvspace + +\drmvspace\drmvspace \begin{align*} - \tsigma_N = \sqrt{\frac{\frac{1}{N} \sum_{i = 1}^{N} z(t_i)^2 - \left( \frac{1}{N} \sum_{i = 1}^{N} z(t_i) \right)^2}{N - 1}} = \frac{\sigma_N}{\sqrt{N}} + \text{Sei }\tsigma_N = \sqrt{\frac{\frac{1}{N} \sum_{i = 1}^{N} z(t_i)^2 - \left( \frac{1}{N} \sum_{i = 1}^{N} z(t_i) \right)^2}{N - 1}} = \frac{\sigma_N}{\sqrt{N}} \end{align*} + +% ──────────────────────────────────────────────────────────────────── + % TODO: Consider adding some of the theory on random variables here (especially consider normal distribution) +% ──────────────────────────────────────────────────────────────────── Das Monte-Carlo-Verfahren beruht auf folgendem: \rmvspace \begin{align*} - \int_{[0, 1]^d} z(x) \dx x = \E z(\mathcal{X}) \text{ mit } \mathcal{X} \sim \mathcal{U}([0, 1]^d) + \int_{[0, 1]^d} z(x) \dx x = \E z(\cX) \text{ mit } \cX \sim \cU([0, 1]^d) \end{align*} -% TODO: Clarify what \mathcal{U} is \drmvspace +wobei $\cU([0, 1]^d)$ die uniforme Verteilung der Zufallsvariable auf dem $d$-dimensionalen Intervall $[0, 1]^d$ ist. + Das Ziel der Monte-Carlo-Methode ist es, den Erwartungswert durch den Mittelwert der Funktionswerte der simulierten Zufallsvariable mit einem Schätzer $m_N(z(\mathcal{X}))$, bzw. einer Schätzung $m_N(z(x))$ zu approximieren: \rmvspace \begin{align*} m_N(z(\cX)) & := \frac{1}{N} \sum_{i = 1}^{N} z(\cX_i) & - m_N(z(x)) & := \frac{1}{N} \sum_{i = 1}^{N} z(x_i) + m_N(z(x)) & := \frac{1}{N} \sum_{i = 1}^{N} z(x_i) \end{align*} \setLabelNumber{all}{16} @@ -48,8 +66,82 @@ Das Ziel der Monte-Carlo-Methode ist es, den Erwartungswert durch den Mittelwert \end{align*} \drmvspace -Die Approximation ist besser, je kleiner die Varianz ist: +Die Approximation von $\E z(\cX)$ durch $m_N(x)$ ist besser, je kleiner die Varianz ist: \rmvspace \begin{align*} \V m_N(z(\cX)) = \V \left( \frac{1}{N} \sum_{i = 1}^{N} z(\cX_i) \right) = \frac{1}{N^2} N \V(z(\cX)) = \frac{1}{N} \V(z(\cX)) \rightarrow 0 \end{align*} + +\drmvspace +Der Zentrale Grenzwertsatz (Central Limit Theorem) besagt, dass für grosse $N$ +\rmvspace +\begin{align*} + \frac{m_N(z(\cX)) - \E z(\cX)}{\sqrt{\frac{1}{N} \V(z(\cX))}} +\end{align*} + +\drmvspace +sich fast wie eine normalverteilte Zufallsvariable $\cY \sim \cN(0, 1)$ verhält. +Daraus folgt mit $\sigma(z(\cX)) = \sqrt{\V(z(\cX))}$: +\rmvspace +\begin{align*} + \Pr\left[ |m_N(z(\cX)) - \E z(\cX)| \leq \frac{\lambda \sigma}{\sqrt{N}} \right] + = \frac{1}{\sqrt{2\pi}} \int_{-\lambda}^{\lambda} e^{-\frac{x^2}{2}} \dx x + \tco{\frac{1}{\sqrt{N}}} + = p(\lambda) + \tco{\frac{1}{N}} +\end{align*} + +\drmvspace +wobei $2 \lambda$ die Länge des Integrationsintervalls ist, mit $\lambda$ definiert als die Länge des untersuchten Intervalls. +Wir haben also eine langsame Konvergenz mit Fehler in $\tco{\frac{1}{\sqrt{N}}}$. + +Oben (Bemerkung \ref{all:5-8-1}) wurde bereits erwähnt, dass wir $68.3\%$ als die Wahrscheinlichkeit haben, dass wir den exakten Wert in unserem Intervall haben. +Woher kommt aber dieser Wert? +Wenn wir $\lambda = 1$ wählen (wie es der Fall ist für das gewählte Intervall), dann erhalten wir $p(1) \approx 0.683$. +Dies trifft allerdings nur zu, wenn $N$ genügen gross ($N \geq 30$) ist. + +Um die Präzision abzuschätzen, benötigen wir einen Schätzer für $\V(z(\cX)) = \E(z(\cX)^2) - (\E(z(\cX)))^2$: +\rmvspace +\begin{align*} + d^*_2(z(\cX)) := \frac{1}{N - 1} \sum_{j = 1}^{N} (z(\cX_j) - m_N(z(\cX)))^2 = \frac{N}{N - 1} \left( m_N(z(\cX))^2 - (m_N(z(\cX)))^2 \right) +\end{align*} + +\drmvspace +Aus dem kann die Schätzung für $\tsigma_N$ von oben hergeleitet werden: +\rmvspace +\begin{align*} + \sqrt{\frac{N}{N - 1} \left( m_N(z(x)^2)) - (m_N(z(x)))^2 \right)} +\end{align*} + +\drmvspace für $y = \sqrt{\V(z(\cX))}$ (oder ohne das $N$ im Zähler für $\frac{y}{\sqrt{N}}$) und $\V(z(\cX)) = \E d_2^*(z(\cX))$ + +\begin{theorem}[]{Vertrauensintervall} + Das Intervall mit $y = \lambda d^*_2(z(\cX))$ + \rmvspace + \begin{align*} + \left[ m_N(z(\cX)) - \frac{y}{\sqrt{N}}, m_N(z(\cX)) + \frac{y}{\sqrt{N}} \right] + \end{align*} + + \rmvspace + enthält den wahren Wert $I = \E(z(\cX))$ mit Wahrscheinlichkeit $p(\lambda)$ bis auf einen Fehler in $\tco{\frac{1}{\sqrt{N}}}$ und das Vertrauensintervall ist: + \drmvspace + \begin{align*} + [m_N(z(\cX)) - \lambda \tsigma_N, m_N(z(\cX)) + \lambda \tsigma_N] + \end{align*} +\end{theorem} +Wiederholen wir das Experiment $M$ mal, mit $M$ gross, so enthält das Vertrauensintervall den Wert $I$ in ungefähr $100 p(\lambda)$ Prozent der $M$ Fälle. + +\setLabelNumber{all}{21} +\inlineremark Ein kurzes Vertrauensintervall entspricht einer hohen Wahrscheinlichkeit, dass das Vertrauensintervall den wahren Wert enthält. +Im Grundrezept wird $\lambda = 1$ mit $p(1) = 0.683$ verwendet. +Um ein kürzeres Vertrauensintervall zu erzielen können wir $N$ erhöhen, oder die Berechnungen so reorganisieren, dass die Varianz kleiner wird. +Wir brauchen also für einen Fehler $\varepsilon$ ungefähr $\frac{1}{\varepsilon^2}$ Evaluierungen. + +Seiten 171 bis 176 im Skript enthalten eine Implementation. Unten eine simple Implementation ohne Plotting: +\begin{code}{python} + import numpy as np + def monte_carlo_integral(func, a, b, N): + t = np.random.uniform(a, b, N) + fx = func(t) + I = np.mean(fx) * (b - a) + var = np.std(fx, ddof=1) / np.sqrt(N) + return I, var +\end{code} diff --git a/semester3/numcs/parts/02_quadrature/06_reduction-of-variance.tex b/semester3/numcs/parts/02_quadrature/06_reduction-of-variance.tex index a11f95d..ccba90f 100644 --- a/semester3/numcs/parts/02_quadrature/06_reduction-of-variance.tex +++ b/semester3/numcs/parts/02_quadrature/06_reduction-of-variance.tex @@ -1,2 +1,79 @@ \newsection \subsection{Methoden zur Reduktion der Varianz} +% NOTE: Mostly from TA slides, as the script is quite convoluted there +Bei höheren Dimensionen $d$ ist die Monte-Carlo-Methode oft die einzige praktikable Option. +Deshalb ist es wichtig, Methoden zu haben, um die Varianz zu verringern. + +\subsubsection{Control Variates} +Die Idee ist hier, bekannte Integrale zu verwenden, um die Varianz zu reduzieren. +Wir schreiben unser Integral unter Verwendung eines bekannten, exakten Integrals $\varphi(x)$ neu: +\rmvspace +\begin{align*} + \iiint f(x) = \iiint (z(x) - \varphi(x)) = \iiint z(x) + \iiint \varphi(x) +\end{align*} + +\drmvspace +Das Ganze funktioniert natürlich für jedes $d \in \N$. +Oft wird die Taylor-Entwicklung von $f(x)$ gewählt, da diese einfach analytisch integrierbar ist. + +Die Varianz wird dadurch reduziert, dass wir nur noch für $z(x)$ einen Fehler haben. + +\innumpy können wir dies folgendermassen implementieren: +\begin{code}{python} + def control_variate_mc(func, phi, analytic_int_phi, a, b, N): + t = np.random.uniform(a, b, N) + val = func(t) - phi(t) + I = np.mean(val) * (b - a) + analytic_int_phi + return I +\end{code} + + +\subsubsection{Importance Sampling} +\begin{intuition}[]{Importance Sampling} + \begin{itemize} + \item Nicht alle Punkte, die während der Monte-Carlo Integration gezogen werden sind, sind gleich wichtig + \item Importance Sampling optimiert die Verteilung der Punkte + \item Man gewichtet die Punkte mit einer Dichtefunktion $g(x)$, die wichtige Bereiche betont + \item Der Erwartungswert wird als gewichteter Mittelwert berechnet, sodass keine Verzerrung auftritt + \end{itemize} +\end{intuition} +Wir schreiben unser zu berechnendes Integral mit $D = [0, 1]^d$ ein Intervall +\rmvspace +\begin{align*} + I = \int_D f(x) \dx x +\end{align*} + +\drmvspace +mit der Hilfsdichte $g(x)$ (für welche gilt $\int_{D} \dx = 1$) +\rmvspace +\begin{align*} + I = \int_D \frac{z(x)}{g(x)} g(x) \dx x = \E_g \left( \frac{z(\cX)}{g(\cX)} \right) +\end{align*} + +\drmvspace +Der entsprechende Monte-Carlo-Schätzer mit $N$ Stichproben $\cX_i \sim g$ ist +\rmvspace +\begin{align*} + \hat{I}_N = \frac{1}{N} \sum_{i = 1}^{N} \frac{z(\cX_i)}{g(\cX_i)} +\end{align*} + +\drmvspace +und dessen Varianz ist +\rmvspace +\begin{align*} + \V_g\left( \frac{z(\cX)}{g(\cX)} \right) = \int_{D} \frac{z^2(x)}{g(x)} \dx x - I^2 +\end{align*} + +\drmvspace +Ideal ist $g(x) \propto |f(x)|$, also proportional zum Betrag von $f(x)$ + + + +\subsubsection{Quasi-Monte-Carlo} +Oft ist ein deterministischer Fehler nützlich, weshalb man bei der Quasi-Monte-Carlo-Methode die Zufallszahlen durch quasi-zufällige Folgen ersetzt. +Diese Folgen decken den Integrationsbereich systematisch ab. + +Dies führt dazu, dass unser Fehler mit $\tco{N^{-1} \cdot (\log(N))^d}$ abnimmt. +Für kleine $d$ haben wir ungefähr eine Abnahme in $\tco{N^1}$, aber bei grossen $d$ ist die Verbesserung kaum mehr sichtbar. + +In der Realität sind diese Methoden (besonders die Sobol-Sequenzen) trotzdem effektiv, da viele Integrale ``effektiv niedrigdimensional'' sind.