Scientific Computing, Spring 1999

Assignment 8.

Given April 27, due end of finals week.

Objective: Go throught the basics of Monte Carlo

1. Write a program to generate standard normal random variables using the Box Muller algorithm. One call to your routine should make a large number of normals, say ten thousand. Verify the correctness of the program using a histogram. Choose a small but not too small Δx and let the kth bin be the interval

B_{k} = (x_{k} − \frac{Δx}{2}, x_{k} + \frac{Δx}{2}) where x_{k} = kΔx.

Generate a large number, N, of random variables and let N_{k} be the number of them that land in B_{k}. The expected value of N_{k} is

E (N_{k}) = N · Pr (X ∈ B_{k})
= N · ∫_{B_{k}} ρ(x)dx .

For Δx small, the integral is well enough approximated by a panel approximation, with B_{k} being the only panel; try the midpoint rule or 2 point Gauss quadrature. The variance of N_{k} is Np(1 − p) where p = Pr(X ∈ B_{k}). Use this to put error bars on the histogram. Choose a sensible range of k values, and N large enough so that there is relatively little error in the largest bin counts.

2. For a parameter β > 0, we want to sample the one dimensional density

ρ_{β}(x) = \frac{1}{Z(β)}e^{−βh(x)} , where h(x) = \frac{1}{2} (x^{2} + x^{4}) .

Write a program to do this by rejection from a gaussian with appropriate β:

ρ_{0β}(x) = \frac{1}{Z_{0}(β)}e^{−βh_{0}(x)} , where h_{0}(x) = \frac{1}{2}x^{2} .

Verify that the generator is correct using the histogram method from part 1. For this verification, and only for this, you need to find Z(β). Comment on how the efficiency depends on β.

3. Now define

S_{n} = Δx (X_{1} + ⋯ + X_{n}) , β(s) = \frac{1}{2} + s^{2} , β_{k} = β(S_{k}) ,

and choose X_{k+1} ∼ ρ_{β_{k}}, with the convention that β_{0} = 1/2. This corresponds to a rondom walk in a variable temperature environment. Estimate

E (S_{n}^{2}) with n = 100, Δx = .1.

Error bars are required.

4. (Extra credit, do as time permits.) We want an accurate estimate of the very small

Pr (S_{n} ≥ r)

with r = 4, Δx = .1, and n = 100. Do this by importance sampling, replacing h(x) by

h_{μ}(x) = \frac{1}{2} (x^{2} + x^{4} − μx) .

Experiment with μ to find a value that gives the largest variance reduction.

5. (Time permitting but with no class credit.) Relax, you’re done.