Given Feb 25; Due March 4


1. Gibbs Sampler to study fertility changes:

    Fertility fluctuations F_1, F_2, ..., F_m over a sequence of fields
    are modeled with the joint density

    p(f_1,...,f_m) = constant exp[ -0.5 psi sum_(i~j) ( f_i - f_j )^2 ]

    Here the sum is over all adjacent fields, and the fields are in
    a linear order, so, for example, field 3 is adjacent to fields 2 and 4.
    Also, the precision psi is known to be psi=1/2.  Crop yields
    Y_1, ..., Y_m are modeled as conditionally indpendent given fertilities,
    with   p(y_j|f_j) = Normal[ mean=(mu+f_j), variance=sigma^2 ]
    and further both the overall mean mu and variance sigma^2 are known
    to equal 10 and 1 respectively.

    a. Show that the full conditional distribution of F_j is normal,
       and compute its mean and variance.

    b. There are m=13 fields and we have observed yields

      10.83,7.48,7.34,7.12,8.70,7.57,6.61,8.25,6.91,5.73,6.60,5.30,4.87
 
       in some standardized units.
       Consider the Gibbs sampler for simulating the conditional
       distribution p(f_1,...,f_m|y_1,...,y_m).  Note you may do this
       by imagining the Gibbs sampler simulating the joint distribution,
       but in which you fail to update any of the y_i's, leaving them at
       their observed values.

       i  Write code to run a systematic scan Gibbs sampler over the 
          fertilities.
       ii Implement the sampler from the starting position (0,...,0)
          and run it for at least B=1000 scans.  
          Use the output to approximate the conditional distributions
          p(f_j|y_1,...,y_m) for each j.  You may calculate means and
          variances, and/or report histograms in a compact way.
      iii Focus on one variable, say F_1.  Use the sampler above again
          to obtain 1000 realized f_1's.  Calculate the sample autocorrelation
          of this time series.  
       iv  Run the sampler over the 13 fertilities as above, but until B=5 only.
           Repeat this 100 times. Compare the 5 resulting distributions of F_1,
           for example using parallel boxplots.  Comment.
        v  Repeat (iv) using the starting position seq(5,-5,length=13).