/* mh2.sas                                       */
 /* second trial of Metropolis-Hastings algorithm */
 /*  using log-series (modified) likelihood       */
 /*  and random walk, uniform trial distribution  */
 /*                                               */
  options ls = 80 ;                     * easier to read output ;
  data a b ;                            * two datasets ;
 /*                       initial values          */
  told = .5 ;
  pstarold  = 0 ;                       * yes, it's wrong ;
  kaccept = 0 ;
  seed = 5151917 ;
  pstarnew = ranuni(seed) ;
 /*                       main loop               */
  do i = 1 to 16384 ;
  tnew = told + (ranuni(seed)-.5)/2.5 ; * add unif( -.2, +.2 ) ;
  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 ;
  if( i le 2048 ) then output a ;        * first 1/8 of data in a ;
  if( i gt 8192 ) then output b ;        * last 1/2 of data in b ;
  end ; * of loop on i ;
run ;
 /*                       some simple statistics  */
 /*                   on first part (1/8) of data */
  proc means data=a n mean var max ;
    var told tnew kaccept ;
    title 'Metropolis-Hastings on log-series posterior' ;
    title2 'demonstration mh2 -- uniform(-.2,+.2) random walk' ;
    title3 'some simple statistics -- dataset a' ;
run ;
 /*                       spectral density est    */
  proc spectra data=a out=ca s adjmean ;
    weights truncat 45 0 ;               * d = 45, constant weight ;
    var told ;
run ;
  data cca ;
    set ca ;
    qerf = freq / 6.2831853 ;            * adjust by 2*pi ;
    sf = s_01*6.2831853 ;                * on both ;
run ;
  proc print data=cca (obs=5) ;
    var qerf sf ;
    title3 'spectral density estimates at small frequencies -- a' ;
run ;
 /*                       some simple statistics  */
 /*                        on last half of data   */
  proc means data=b n mean var max ;
    var told tnew kaccept ;
    title3 'some simple statistics -- dataset b' ;
run ;
 /*                       spectral density est    */
  proc spectra data=b out=cb s adjmean ;
    weights truncat 90 0 ;               * d = 90, constant weight ;
    var told ;
run ;
  data ccb ;
    set cb ;
    qerf = freq / 6.2831853 ;            * adjust by 2*pi ;
    sf = s_01*6.2831853 ;                * on both ;
run ;
  proc print data=ccb (obs=5) ;
    var qerf sf ;
    title3 'spectral density estimates at small frequencies -- b' ;
run ;