Fourth Problem Set

Due Monday, October 22, 2007

Here's a Metropolis-Hastings sampler problem to be done in R. I encourage you to work in groups. Email me the code, and turn in hard copy including all code, plots and comments.

Suppose you have a data set that is supposed to be fit by a straight line y_k=a+bx_k, where the error in y_k is Cauchy rather than normal. We presume that the error in y_k is drawn from a Cauchy distribution with an unknown scale parameter s and location parameter a+bx_k (i.e., in R construct the log likelihood from sum(dcauchy(y,a+b*x,s,log=TRUE)) where x and y are given below). I have generated some data you can use at the end of the assignment. You can copy them out of the web page and paste them into your program.

Program a Metropolis-Hastings sampler for this problem, to obtain a sample from the joint posterior distribution of the unknown parameters a, b and s. I suggest that you use a symmetric sampler such as we discussed in class, where you propose a move to a distance +/- da, db or ds from your current location, drawing the proposal from a uniform distribution centered on your current location. Choose appropriate priors for a, b and s and discuss the reasons for your choices. Run your sampler for a number of choices of (a) starting points, (b) proposal step size (may be different for different parameters), seeking to find optimal values for each that gives a sampler with good mixing properties. Also, plot summary information such as: histograms of posterior marginals; sampling history plots; plots of the samples of a vs. b, a vs. s and b vs. s. Compute appropriate credible intervals for a, b and s, as well as posterior means, medians, and variances.

Hint: Since the scale parameter s>0, you may generate a proposal for s that violates this constraint. You should check whether the constraint is violated before attempting to evaluate the log likelihood, or you may generate an exception! If the constraint is violated, just reject the proposed step. In effect, you would be setting the prior to 0 whenever s is 0 or less; so if a proposal took you into this forbidden territory, the M-H factor would be zero and you wouldn't accept the step.

Hint: A straightforward way to do this is to modify the code I posted for the linear problem. Just plug in the correct likelihood, and delete the code that samples on the X variables, since this is not an errors-in-variables problem.

Comment: There are several ways to design the sampler. You can sample a, b and s in sequence, or you can sample them simultaneously, or even sample some parameters simultaneously but in sequence with others (e.g., sample a and b simultaneously, then sample s). It's your choice, but be sure to discuss your reasons for your choice. You will probably find it challenging to get good "mixing" with this problem, just do the best you can, but be sure to stay out of regions where the sampler "sticks" or on the other hand makes steps that are much too small. Since this problem doesn't mix well, you may need a larger sample to get adequate estimates.

Discuss your results intelligently. Don't forget to include your code!

DATA SET

x <- c(1.808,1.799,1.179,0.574,3.727,0.946,3.719,1.566,3.596,3.253)

y <- c(1.816,1.281,-1.382,0.573,3.793,0.935,1.775,1.474,3.679,3.889)