/* chex133.sas -- Example 13.3                 */
 /*  pump example -- Metropolis within Gibbs    */
 /*                                             */
  options ls=80 ;                 * easier to read output ;
  data a ;
    array x(10) ;
    array y(10) ;
    array f(10) ;                 * rate parameters ;
    input x1-x10 ;                * operating times ;
    input y1-y10 ;                * failures ;
 /*                initial values */
    seed = 5151917 ;
    a = .7 ;
    b = .9 ;
    kaccept = 0 ;
 /*                main loop(s)   */
    do i = 1 to 1000 ;            * 1000 replications ;
 /*                first phi's    */
    slf = 0. ;
    sumf = 0. ;
    do j = 1 to 10 ;
    f(j) = rangam(seed, y(j)+a ) / (x(j)+b)  ;
    slf = slf + log(f(j)) ;
    sumf = sumf + f(j) ;
    end ;  * loop on j ;
 /*                then beta      */
    b = rangam(seed, 10*a+.1) ;   * c = .1, J = 10 ;
    b = b / ( 1 + sumf ) ;
 /*                M-H for alpha  */
 /* pstarold changes each time    */
    pstarold = exp( a*J*log(b) + a*slf - a - 10*lgamma(a) ) ;
 /*                candidate anew */
    anew = a + (ranuni(seed) - .5 ) / 1 ;
    pstarnew = 0. ;
    if( anew > 0 ) then do ;
       pstarnew
         = exp( anew*J*log(b) + anew*slf - anew - 10*lgamma(anew) ) ;
                        end ;
    alpha = 1 ;
    if( pstarnew < pstarold ) then alpha = pstarnew / pstarold ;
 /*                acceptance     */
    if( ranuni(seed) < alpha ) then do ;
 /*                accept         */
                   kaccept = kaccept + 1 ;
                   a = anew ;
                                    end ;
    output ;
    end ;  * of loop on i ;
  cards ;
  94.3  15.7  62.9  126  5.24  31.4  1.05  1.05  2.1  10.5
   5     1     5     14  3     19    1     1     4    22
 ;
run ;
/*        some analysis and plots */
  proc means data=a ;
    var a b f1-f10 kaccept ;
    title 'chex133 -- pump example' ;
run ;
  goptions rotate=landscape ;
  proc gplot data=a ;
    plot b*a ;
run ;
  proc gplot data=a ;
    plot a*i ;
    title2 'trace plot of alpha' ;
run ;
  proc gplot data=a ;
    plot b*i ;
    title2 'trace plot of beta' ;
run ;
  proc gplot data=a (obs=50) ;
    plot b*a ;
    symbol1 c=black i=join ;
    title2 'trace plot of first 50 obs' ;
run ;