Activity: Area Estimation

Monte Carlo Methods Activities

Author

Manuel D. Rossetti, PhD P.E.

We write this program together: the instructor types the code, and you type along. Before each run, write down a prediction; then run and check it.

The question. Let \(f(x) = \sqrt{1 - x^2}\) represent the curve of the quarter circle of radius 1. Let \(\theta\) represent the area under \(f\) from 0 to 1:

\[\theta = \int_0^1 \sqrt{1 - x^2}\,\mathrm{d}x.\]

We estimate \(\theta\) two ways, and ask which way is better.

Setup. Get the starter file AreaEstimation.kt from the code/ folder beside this activity. Put it in src/main/kotlin/student/activities in your KSL project, the package the file declares; create the folder if you do not have it. Run its main. It compiles and runs, but the estimates are NaN: the two functions marked TODO are what we write.

1. The target

  1. What is the exact value of \(\theta\)? (Think about the shape, not the integral.)

  2. The book’s section Simple Monte Carlo Integration lets \(X \sim U(a, b)\) and \(Y = (b - a)\,g(X)\), and shows that \(E[Y]\) is the area under \(g\) from \(a\) to \(b\). Here \(a = 0\), \(b = 1\), and \(g = f\). In words, what will our program average?

2. Crude Monte Carlo

We write crudeMonteCarlo (TODO 1) together.

  1. Before the first run, predict the 95% half-width for \(n = 1000\). You need the standard deviation of \(Y = f(X)\): compute \(\mathrm{Var}[Y] = E[Y^2] - \theta^2\), using \(E[Y^2] = \int_0^1 (1 - x^2)\,\mathrm{d}x\). Then the half-width is about \(1.96\,\sigma/\sqrt{n}\).

  2. Run with \(n = 1000\). Does the 95% confidence interval contain \(\theta\)? How close are the sample standard deviation and the half-width to your predictions?

  3. Predict the half-width for \(n = 4000\). Change n, run again, and check.

3. Hit or miss

A second estimator throws random points \((x, y)\) into the unit square and counts the fraction that lands on or under the curve. We write hitOrMiss (TODO 2) together.

  1. Each point is a hit with probability \(\theta\), so each observation is a Bernoulli random variable. What is its variance? Before running, predict which estimator has the smaller half-width for the same \(n\), and by roughly what factor.

  2. Run with \(n = 1000\), and then with \(n = 4000\). Does each 95% confidence interval contain \(\theta\)? If one does not, is the program wrong?

  3. For the same half-width, how many times as many points does hit or miss need as crude Monte Carlo?

4. How many samples?

  1. The program also prints estimateSampleSize(0.001): the sample size each method needs for a half-width of \(\pm 0.001\), using the formula in the book’s section Determining the Sample Size. Compute the same two numbers from the exact variances. Why do the program’s estimates change between the \(n = 1000\) and the \(n = 4000\) runs?

  2. Which estimator would you use for this integral? Describe a problem where hit or miss is still the natural choice.