Third Problem Set

Due Friday, October 12, 2007

Here's a set of Gibbs sampler problems to be done on the computer.

1. (Casella and George). Consider the following distribution:

f(x,y) = k C(n,x) y^x+a-1(1-y)^n-x+b-1

where k is a fixed normalization constant that you don't need to know for the problem and C(n,x) is the binomial coefficient. For fixed y, this is a binomial distribution on the integer x, and for fixed x it is a beta distribution on the continuous parameter y where y is in the interval [0,1]. You can find out how to use R's binomial and beta functions by executing help(rbinom) and help(rbeta) from within R.

Program a Gibbs sampler for x and y, with n=16, a=2, b=4. Run the Gibbs sampler with several values for the total number of samples. Compute means and 95% HDRs for the distributions on x and y (use R's quantile function for the HDR's). Plot simulation histories and marginal distributions for x and y. WARNING: R's hist function will consolidate the histogram for x=0 and x=1 unless you take special steps. I suggest using

hist(x,0:17-0.5)

to plot the histogram. We can discuss in class how this magic works.

2. (Casella and George)

Consider the following distributions (you don't need to know the normalization constants):

f(x|y) = k exp(-xy) for 0<x<5 and 0 elsewhere
f(y|x) = k exp(-xy) for 0<y<5 and 0 elsewhere

(The reason for the restriction of the range of x is to guarantee the existence of the marginal distributions. You can think of it as due to a prior that is 0 outside this range).

These are truncated exponential distributions. You can simulate draws from a truncated exponential distribution by using the following special function:
rtexp = function(rate=1,cutoff=5) {
    repeat {
        x = rexp(1,rate)
        if(x < cutoff) break
    }
    return(x)
}
where rate is the value of x or y, whichever is fixed, and cutoff=5 (default).

As in problem 1, program a Gibbs sampler for these conditional distributions and run it, computing similar results as in problem 1. You don't need to do anything fancy with the hist function for this, but you may wish to control the number of bins with the second parameter to the function.

I encourage you to work in groups of 2 or 3 on this problem set. Several brains are better than one.