/* mh3.sas                                      */
 /* third trial of Metropolis-Hastings algorithm */
 /*  using log-series (modified) likelihood      */
 /*  and uniform trial distribution              */
 /*                                              */
 /* Gelman-Rubin replications                    */
  options ls=80 ;                       * easier to read output ;
  data a ;
 /*                       initial values         */
  seed = 5151917 ;
  pstarnew = ranuni(seed) ;             * waste one to match .f95 ;
  kaccept = 0 ;
  do rep = 1 to 16 ;
 /*                       rep initialization     */
  told = ranuni(seed) ;                 * uniform random starts ;
  pstarold  = 0 ;                       * yes, it's wrong ;
 /*                       main loop              */
  do i = 1 to 1024 ;
  tnew = ranuni(seed) ;                * acceptance / rejection ;
  pstarnew = 0 ;
  if( (tnew > 0) and (tnew < 1) ) then do ;
             lp = 16*log(tnew) + log(1-tnew) - 10*log(-log(1-tnew)) ;
             if( lp > -35 ) then pstarnew = exp(lp) ;
                                       end ;
  alpha = 1 ;                         * get acceptance prob ;
  if( pstarnew < pstarold ) then alpha = pstarnew / pstarold ;
 /*                                              */
 /*                       acceptance trial       */
  if( ranuni(seed) < alpha ) then do ;
 /*                       acceptance             */
                          kaccept = kaccept + 1 ;
                          told = tnew ;
                          pstarold = pstarnew ;
                          end ;
  output ;
  end ; * of loop on i ;
  end ; * of loop on rep ;
run ;
 /*                       ANOVA computations     */
  proc means data=a mean var ;
    class rep ;
    var told ;
    title 'Demonstration mh3 -- Metropolis-Hastings' ;
    title2 'on log-series posterior, independence chain' ;
    title3 'Gelman-Rubin replications -- means' ;
run ;
  proc glm data=a ;
    class rep ;
    model told = rep ;
    title3 'Gelman-Rubin replications -- ANOVA' ;
run ;
 /*                                              */
 /*                  spectral variance estimates */
  proc spectra data=a out=b s adjmean ;
    weights truncat 32 0 ;
    var told ;
    by rep ;                                     * for each rep ;
run ;
  data bb ;
    set b ;
    if( freq ne 0 ) then delete ;
    sf = s_01*6.2831853 ;
run ;
  proc print data=bb  ;
    *var rep sf ;
    title3 'Gelman-Rubin replications -- spectral var est';
run ;
  proc means data=bb ;
  var sf ;
run ;
proc means max data=a ;
 var kaccept ;
run ;