A Data Scientist working in central Florida with an Applied Math PhD from Northwestern. This space is for code related to my research and recreational projects.
My thesis research has been in rare events in nonlinear optics. This research combines techniques in nonlinear waves, wave propagation, stochastic simulation, and numerical optimization. In a past life I did research in nonlinear water waves, but the most fulfilling parts of my doctoral research have been in stochastic simulation. In particular, I’ve studied a mode-locked laser model and quantified error rates using importance sampling.
Monte Carlo is a numerical simulation technique for calculating integrals and assessing probabilities. It utilizes the fact that the expected value of a random quantity is calculated mathematically via an integral:
\[I=\mathbb{E}_p[f(\mathbf{X})] = \int_{D} f(\mathbf{x})p(\mathbf{x})\,d\mathbf{x}.\]where $f$ is some function and $p$ is the probability distribution that gives rise to the random variable $\mathbf{X}$. Monte Carlo integration works to approximate $I$ by randomly sampling from $p$ and then averaging $f$. For instance, after drawing $N$ samples our approximation to $I$ is
\[\hat{I}_N = \frac{1}{N}\sum_{l=1}^N f(\mathbf{X}_l).\]However, if we’re interested in rare events, it may take an infeasible amount of samples in order to adequately refine the probability. Importance sampling can be used to draw samples from another probability distribution, $p^*$ (called a biasing distribution), which leads to the rare event of interest more frequently. Mathematically, this technique amounts to a simple tweak of the integral:
\[I = \int_{D} f(\mathbf{x})\frac{p(\mathbf{x})}{p^*(\mathbf{x})}p^*(\mathbf{x})\,d\mathbf{x}.\]The Monte Carlo average is then computed via
\[\hat{I}_N = \frac{1}{N}\sum_{l=1}^N f(\mathbf{X}_l)\frac{p(\mathbf{X}_l)}{p^*(\mathbf{X}_l)}.\]where $f$ is now multiplied by the likelihood ratio of the biased noise $L(\mathbf{X}_l)=\frac{p(\mathbf{X}_l)}{p^*(\mathbf{X}_l)}$.
This may sound simple enough, but the difficulty in using importance sampling is the determination of the biasing distribution, which is usually accomplished through some kind of optimization. For instance, in my mode-locked laser research, I studied an escape problem where I was interested in the probability that a pulse of light would escape its bit-slot (or wander outside of the time allotted for its transmission) and cause a communication error. A number of stabilizing physical elements needed to be overcome for this to occur, so the error paths were found to be quite complicated. These paths, found using a combination of asymptotic, analytic, and numerical methods, sometimes involved a large number of oscillations.
This video shows an exit path that goes through two oscillations. The pulse is shown in the top panel, and the drift of the pulse is shown in frequency-position phase space in the bottom panel (the frequency is coupled to the velocity of the pulse so it gives information on the pulse’s movement). For more on importance sampling please view my tutorial repository on Github.