Monte Carlo Methods and Applications

# Monte Carlo Methods: Assignment 03 *Assigned 2023-09-12, due 2023-09-19 on Gradescope* ### Written In general, stratified sampling is guaranteed not to increase the variance of a Monte Carlo estimator—and will typically give a reasonable reduction in variance. A detailed discussion and analysis of stratification can be found in [Owen, §8.4](https://artowen.su.domains/mc/). For this exercise, let's just consider a smooth and bounded function $f: [0,L] \to \mathbb{R}$ (for some interval length $L$), and assume we've chosen $k$ equally-spaced strata of small enough size $h := L/k$ so that \[ V\left[f|_{[(p-1)h,ph]}\right] \leq V[f], \quad p = 1, \ldots, k, \] i.e., such that the variance within each stratum is no greater than the variance of the overall function. (Here $f_{[a,b]}$ denotes the restriction of the function $f$ to the interval $[a,b]$.) Show then that the variance of the stratified estimator \[ \widehat{I}^k_n := \sum_{p=1}^k \frac{h}{n} \sum_{i=1}^n f(X_{p,i}), \quad X_{p,i} \sim \mathcal{U}_{[(p-1)h,ph]} \] (where we take $n$ samples chosen uniformly from each of the strata) is never greater than the variance of the basic Monte Carlo estimator \[ \widehat{I}_N = \frac{L}{N} \sum_{i=1}^N f(X_i), \quad X_i \sim \mathcal{U}_{[0,L]}, \] assuming we take the same total number of samples ($N = kn$). Note that $\mathcal{U}_{[a,b]}$ denotes the uniform distribution on the interval $[a,b]$. *Notes*: The condition on the variance of $f$ in each stratum isn't necessary but simplifies the problem. You're welcome to try a more general proof that doesn't rely on this condition (the statement remains true without it.) ### Coding Download this Jupyter notebook: [download here](https://www.math.cmu.edu/~mcm-cmu/auth/render.ipynb). This week's coding assignment addresses variance reduction techniques in the context of _direct lighting_, a classic problem in physically based rendering -- so you will literally be able to see the results of variance reduction via rendered images. All instructions are contained within text cells in the notebook. When you're done implementing all the specified functions in the notebook, submit your completed notebook on Gradescope **with your images rendered in the notebook**. **Useful clarifications for this assignment** (which have also been added inside the notebook): * It is difficult to see the improvements in the rendered images when they are not viewed side-by-side. I recommend simply downloading the images locally (via "Save Image As") and flicking through the images in quick succession. The newest version of the Python notebook (linked below) also will automatically export rendered images to PNG. * All the sampling methods you implement should return samples in the range $[0,1)^d$, where $d$ is the appropriate dimension ($d=1$ or $d=2$ in this assignment.) For example, the stratified sampler divides $[0,1)^2$ into $\sqrt{N}$ intervals along each axis to yield $N$ total strata. * For all sampling methods in this assignment, use the random number generator stored in the member variable `self.rng` as the RNG (but don't worry about it if you've already completed the assignment by spawning your own.) * We covered Latin Hypercube sampling (and stratfied sampling) in the lecture called "Variance Reduction". For more details about Latin Hypercube sampling (LHS), see [the Wikipedia page](https://en.wikipedia.org/wiki/Latin_hypercube_sampling) or [Chapter 7.2 of PBRT](https://pbr-book.org/3ed-2018/Sampling_and_Reconstruction/Stratified_Sampling), which also have some helpful pictures. In short, the assignment calls for you to generate $N$ 2D samples using LHS, which means that conceptually you are generating random samples in a $N\times N$ grid partition of the unit square such that there is exactly 1 sample in each row and column of the grid. You can do this by first generating a random sample in each diagonal box of the grid, then permuting the rows of this grid, then permuting the columns. * We did not get to covering the [Fisher-Yates shuffle algorithm](https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle#The_modern_algorithm) in lecture, but the algorithm is very simple and I have linked to the relevant Wikipedia section. In particular, the algorithm shuffles in array in-place in linear time, with each possible permutation of the array being equally likely. Here is also some clarification about the multiple importance sampling part of the assignment (also repeated in the Python notebook). In class, we saw the expression for the MIS estimator assuming a _multiple sample model_ , where we generate $N_k$ samples using the $k$th sampling strategy, for each of the $s$ strategies: $$ \widehat{I}_{\rm MIS,\ multi} := \sum_{k=1}^s \frac{1}{n_k}\sum_{i=1}^{n_k} w_k(X_{k,i})\frac{f(X_{k,i})}{p_k(X_{k,i})} $$ where $$ \begin{align*} s &= \text{the number of sampling strategies} \\ n_k &= \text{the number of samples for the $k$th strategy} \\ X_{k,i} &= \text{the $i$th sample drawn using the $k$th strategy} \\ p_k &= \text{the probability density corresponding to the $k$th strategy}. \end{align*} $$ Alternatively, we may assume a _one sample model_ where we generate only one sample, by first choosing one of the $s$ strategies at random and taking a sample from the chosen distribution. This corresponds to the following estimator: $$ \widehat{I}_{\rm MIS,\ single} := \frac{w_I(X_I)f(X_I)}{c_Ip_I(X_I)} $$ where $I\in \{1,\dots,s\}$ is a random variable distributed according to the probabilities $c_k$, and $X_I$ is a sample from the strategy $p_I$. A popular choice for the formula for the weights $w_k$ is the _balance heuristic_: $$ w_k(x) = \frac{p_k(x)}{\sum_{l=1}^s p_l(x)}. $$ In this problem, we will implement the one sample MIS estimator using the balance heuristic, in which case the expression for an estimate becomes $$ \widehat{I}_{\rm MIS,\ single} := \frac{f(X_I)}{c_I\sum_{l=1}^s p_l(X_I)}. $$ In the Python notebook, you will be implementing the above formula. Additional references for the coding portion of the assignment: * [Ch. 7 in _Physically Based Rendering: From Theory to Implementation_ (Pharr, Jakob, Humphreys)](https://pbr-book.org/3ed-2018/contents) (contains details further details about stratified sampling, Halton sequences, and other types of sampling) * Linear congruential generators: [Ch. 3, Vol. 2 in _The Art of Computer Programming_ (Knuth)](https://github.com/djtrack16/thyme/blob/master/computer%20science/Donald.E.Knuth.The.Art.of.Computer.Programming.Volume.2.pdf), [Wikipedia page](https://en.wikipedia.org/wiki/Linear_congruential_generator) * Multiple importance sampling: [Ch. 9 of Eric Veach's PhD thesis](https://graphics.stanford.edu/papers/veach_thesis/thesis.pdf). This thesis is also a great reference for all manner of variance reduction techniques in rendering, including sampling methods. * Generalizations of stratified and Latin Hypercube sampling: [_Orthogonal Array Sampling for Monte Carlo Rendering_ (Jarosz et al.)](https://cs.dartmouth.edu/~wjarosz/publications/jarosz19orthogonal.pdf)