/* mh1.sas                                      */
 /* first trial of Metropolis-Hastings algorithm */
 /*  using log-series (modified) likelihood      */
 /*  and random walk, uniform trial distribution */
 /*                                              */
  option ls = 80 ;                      * easier to read output ;
  data a ;
 /*                       initial values         */
  told = .5 ;
  pstarold  = 0 ;                       * yes, it's wrong ;
  kaccept = 0 ;
  seed = 5151917 ;
  pstarnew = ranuni(seed) ;
 /*                       main loop              */
  do i = 1 to 16384 ;
  batch = int( (i-1)/512 );              * 32 batches ;
  tnew = told + (ranuni(seed) - .5)/5 ;  * add unif( -.1, +.1 ) ;
  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 ;
run ;
 /*                       some simple statistics */
  proc means data=a n mean var max ;
    var told tnew kaccept ;
    title 'Metropolis-Hastings on log-series posterior' ;
    title2 'demonstration mh1 -- uniform(-.1,+.1) random walk' ;
    title3 'some simple statistics' ;
run ;
 /*                       analyze as AR(1)       */
  proc arima data=a ;
    identify var=told nlag=8 ;
    estimate p=1 ;
    title3 'analyze as AR(1)' ;
run ;
 /*                       spectral density est   */
  proc spectra data=a out=b s adjmean ;
    weights truncat 128 0 ;              * d = 128, constant weight ;
    var told ;
run ;
  data bb ;
    set b ;
    qerf = freq / 6.2831853 ;            * adjust by 2*pi ;
    sf = s_01*6.2831853 ;                *  on both ;
run ;
  proc print data=bb (obs=5) ;
    var qerf sf ;
    title3 'spectral density estimates at small frequencies' ;
run ;
 /*                       analyze batches        */
  proc means mean var data=a noprint nway ;
    class batch ;
    var told ;
    output out=c mean=mean var=var ;     * c holds batch means, vars ; 
run ;
  proc means n mean var data=c ;
    var mean var ;                       * analyze batch means, vars ;
    title3 'statistics on batch statistics' ;
run ;
 /*                 analyze batch means as AR(1) */
  proc arima data=c ;
    identify var=mean nlag=8 ;
    estimate p=1 ;
run ;
 /*                 also analyze batch variances */
  proc arima data=c ;
    identify var=var nlag=8 ;
    estimate p=1 ;
run ;