- CA5 due Monday, June 4 (Dropbox)
- Final projects due Friday June 8 (Dropbox)
Add this to the end of CA5:
8) Simulate a 5-pixel data set with flux per pixel described by F(z) = a*(z-1)^2 + b*(z-1) + c, with a=3, b=0.2, c=4. Then use an MCMC routine to produce the two-dimensional PDF for {a,b,c}, as well as the individual 1-dimensional PDFs for a, b and c.
The basic algorithm for MCMC is as follows:
- Evaluate the likelihood of the data at an initial, guessed set of parameters {a,b,c}. Call this L0
- Evaluate a new likelihood L1 at a new set of parameters, say {a+delta,b,c}, where delta is a small perturbation drawn from a "jump distribution." This jump is usually just a call to a normally-distributed random number generator centered on the current value. The parameter to perturb is also drawn at random (there are other ways to do this, but this is the easiest in my experience).
- If L1/L0 > 1, then add the new set of parameters to the chain (array of parameters). If L1/L0 < 1, then add the new parameters to the chain if a uniform random number is less than L1/L0. Else, just add the previous set of parameters to the chain.
- Goto 1 until the desired number of chain links is calculated (usually 1e5 to 1e7 links).
For three parameters, the array of chains will be a 3xN array, where N is the number of chains run. Each chain turns out to be the marginalized, posterior pdf of that parameter, which is pretty cool.
For more info, check out Wikipedia, or Appendix A of this classic astronomy article:
http://arxiv.org/pdf/astro-ph/0310723v2.pdf
I will be available during office hours Wednesday (10am-11am) to go over MCMC in greater detail.
No comments:
Post a Comment