22S:138 Bayesian Statistics
                           Instructor: Cowles
                                  Lab 4
                               Oct. 7, 2005
             Assessing MCMC Sampler Convergence with WinBUGS

1.  Copy the WinBUGS code called "baseball" from the course web page 
    (Handouts section) into a WinBUGS window.

    We will first fit Model 1 with very weak priors on alpha and beta.

2.  We are preparing to use the Brooks, Gelman and Rubin convergence 
    diagnostic to attempt to determine whether our sampler has converged 
    to the target distribution.  It requires running more than one chain 
    using
     a.  the same model for each chain
     b.  different ("overdispersed") initial values for each chain

3.  Model specification
     a.  Check model
     b.  Load data
     c.  Before you compile, change the number of chains to 3 
         (or however many you want)
     d.  Compile
     e.  Load inits for each chain
     f.  If there were any unitialized nodes, Gen inits once.  This will 
         generate the needed initial values for all the chains.

4.  Running and monitoring the sampler
    We want to observe the sampler output from the first iterations.  Thus 
    we must "set" each of the variables before we do any updating.

    We will then assess how many initial iterations need to be discarded.

     a.  Inference / Samples
         Set   
               p
               alpha
               beta
               betamean

     b.  Model / Update
         Update for 1000 iterations

5.  Processing sampler output
     a.  To request output for all the monitored nodes, go to the Sample 
         Monitor tool and put an asterisk ( *) in the node window.
     b.  To put all output in the same window, use Options / Use log.

6.  Output that is useful in assessing convergence
     a.  History (for all chains superimposed in one graph per 
         parameter, and for separate chains)

         Ideally, we want to see:
             All chains settling into the same range regardless of where
                  they started
             Random "white noise" type pattern within each chain after an
                  initial transient

     b.  Autocorrelation plot (for each separate chain)

         Ideally, we would like to see the autocorrelations drop off very
         rapidly for all parameters.

         If autocorrelations are large,
             the sampler will need many iterations to get close to the
                   target distribution
             once it gets there, we will need to draw a very large number
                   of additional iterations in order to get accurate
                   estimation of features of the posterior
             Monte Carlo sampling error will be high

         Tricks that sometimes may reduce autocorrelation in an MCMC sampler:
             reparameterize the model
             use more precise priors (if this is scientifically justifiable)

     c.  Gelman and Rubin diagnostic (for all chains together)
         Select GR diag.
         Convert plots to numeric output by double-clicking on the plot 
         and then control-left-mouse-click on the window.  

         We want to see:
            convergence of BGR R close to 1  
                rule of thumb:  R no larger than 1.1 or 1.2 for any parameter
            convergence of both the pooled- and within- interval widths to 
                 stability
            
         Even if this happens, it is not a guarantee that the samplers have
         found the true target distribution.  But if it doesn't happen, we
         know they have not.

         Solution:
             run more iterations
             use tricks to reduce autocorrelation if it is high

7.  Also get the Stats output 
    a.  for all chains combined
    b.  for each chain separately
    c.  Notice the Monte Carlo error, and check whether posterior means
        and intervals agree for the different chains. 

8.  Go back to the Update tool and run an additional 4000 iterations.  Then
    do the same checks as in steps 6 and 7 above.

9.  Now we will observe the effect of putting stronger priors on alpha
    and beta.

    Make the following change to the model code:

    alpha ~ dgamma( 0.11, 0.11) ;
    beta ~ dgamma(1.0, 0.33 ) ;
    #alpha ~ dgamma(0.0001, 0.0001) ;
    #beta ~ dgamma(0.0001, 0.0001) ;

10.  Repeat steps 3 - 6 with the new prior.    

11.  Model 2 is a reparameterization of the model.  

     Why would a normal prior be a bad choice for the prior on the p[i]'s?


     What is the effect of the logit transformation?


     I have given you three sets of initial values.  Note that we could have
     omitted the initial values for mu, but it was necessary to provide
     initial values for tausq because it is a precision parameter.