/* The %BayesIntervals Macro */

/* This macro computes Bayesian simultaneous confidence intervals. */


/*---------------------------------------------------------------*/
/* Name:      BayesIntervals                                     */
/* Title:     Bayesian simultaneous intervals based on           */
/*            percentiles                                        */
/* Author:    Russell D. Wolfinger, sasrdw@sas.com               */
/* Release:   Version 6 or later                                 */
/*---------------------------------------------------------------*/
/* Inputs:                                                       */
/*                                                               */
/*    data=   data set representing a sample from the            */
/*            posterior distribution (required)                  */
/*                                                               */
/*    vars=   the variables in the data set for which to compute */
/*            simultaneous intervals (required)                  */
/*                                                               */
/*    alpha=  the joint significance level, default=0.05         */
/*                                                               */
/*    tail=   the tail type of the intervals, must be L, U, or   */
/*            2, default=2                                       */
/*                                                               */
/*    maxit=  maximum number of iterations in the search,        */
/*            default=50                                         */
/*                                                               */
/*    tol=    convergence tolerance, default=0.001               */
/*                                                               */
/* Output:                                                       */
/*                                                               */
/*    The macro begins by computing the joint coverage of the    */
/*    naive unadjusted alpha-level intervals, computed as        */
/*    percentiles across the sample.  It then decreases alpha    */
/*    using a bisection search until the joint coverage comes    */
/*    within the tol= value of alpha.  A history of this search  */
/*    is printed and then the simultaneous intervals.            */  
/*                                                               */
/*---------------------------------------------------------------*/

%macro BayesIntervals(data=&syslast, out=BayesIntervals, 
    var=, vars=&var, alpha=.05, tail=2, maxit=50, tol=0.001);

 %let data=%unquote(&data);
    
 options nonotes nomprint;

 /*---do bisection search to find adjusted alpha---*/
 %let lowera = 0;
 data _null_;
    uppera = &alpha * 2;
    call symput('uppera',uppera);
 run;
 %let iter = 0;

 %put %str(        )The BayesIntervals Macro;
 %put;
 %put Iteration    Alpha          Coverage;

 %do %while(&iter < &maxit);

    /*---compute quantiles---*/
    data _null_;
       alf = (&lowera + &uppera)/2;
       %if (&tail = 2) %then %do;
          lowerp = 100 * alf/2;
          upperp = 100*(1 - alf/2);
       %end;
       %else %if (&tail = L) %then %do;
          lowerp = 100 * alf;
          upperp = 100;
       %end;
       %else %do;
          lowerp = 0;
          upperp = 100*(1 - alf);
       %end;
       call symput('alf',left(alf));
       call symput('lowerp',left(lowerp));
       call symput('upperp',left(upperp));
    run;
  
    proc univariate data=&data pctldef=1 noprint;
       var &vars;
       output pctlpts=&lowerp,&upperp pctlpre=&vars pctlname=l u out=p;
    run;
  
    proc transpose data=p out=pt;
    run;
  
    /*---load limits and variable names into macro variables---*/
    data _null_;
       set pt nobs=count end=last;
       retain i 0;
       if (mod(_n_,2)=1) then do;
          i = i + 1;
          mname = "v" || left(put(i,8.));      
          len = length(_NAME_);
          vname = substr(_NAME_,1,len-1);
          call symput(mname,left(vname));
          mname = "lv" || left(put(i,8.));      
          call symput(mname,left(put(COL1,best8.)));
       end;
       else do;
          mname = "uv" || left(put(i,8.));      
          call symput(mname,left(put(COL1,best8.)));
       end;
       if last then do;
          call symput('nvar',left(put(count/2,8.)));
       end;
    run;
  
    /*---pass through data and determine simultaneous coverage---*/
    data _null_;
       set &data nobs=no end=last;
       retain count 0;
       bad = 0;
       %do i = 1 %to &nvar;
          if (&&v&i < &&lv&i) or (&&v&i > &&uv&i) then bad = 1;
       %end;
       if (bad = 0) then count = count + 1;
       if last then do;
          coverage = count / no;
          target = 1 - &alpha;
          if (abs(coverage - target) < &tol) then conv = 1;
          else do;
             conv = 0;
             alf = &alf;
             if (coverage < target) then do;
                call symput('uppera',left(alf));
             end;
             else do;
                call symput('lowera',left(alf));
             end;
          end;
          call symput('conv',left(conv));
          call symput('coverage',left(coverage));
       end;
    run;

    %let iter = %eval(&iter+1);
    %put  %str(   ) &iter %str(      ) &alf %str(  ) &coverage;
    %if (&conv=1) %then %let iter=&maxit;

 %end;

 options notes;
 %if (&conv=1) %then %do;
    data &out;
       %do i = 1 %to &nvar;
          _NAME_ = "&&v&i";
          Lower = &&lv&i;
          Upper = &&uv&i;
          output;
       %end;
    run;
    
    proc print data=&out;
    run;
 %end;
 %else %do;
    %put Did not converge;
 %end;

 options notes mprint;

%mend BayesIntervals;