goptions reset=all;
* filename grafout 'c:grafout\boxes5.ps';

goptions
  targetdevice=pscolor

      lfactor=1 htext=1.5  vsize=6 in;
*********************************************
** Macro for letter with subscript         **
********************************************;

  %macro greek(x,nn);
  function='label';
  position = '5'; style = 'cgreek'; text = &x ; size=2; output;
  position = '9'; style = 'cgreek';x=x+.04; if xsys='5' then x=x+1;
  text = &nn; size=1.6; output;
  position = '5';
  function = 'draw';  size=1;  style='duplex';
  %mend greek;




************************************************
****              CI MACRO                  ****
****   macro to get ellipsoid and slices    ****
************************************************;


   %macro ci(id, ns, color, pct);

    /*-------------------------------------------------------------------
    |   Z is id=1 or 4, X is ID=2 or 5, Y is ID=3 or 6                   |
    |       where PLOT Y*X=Z is used in G3D **;                          |
    |   ns is number of slices you want on ID variable. Color of slice   |
    |   pct controls how far first slice is from end of ellipsoid        |
     -------------------------------------------------------------------*/
     ***  Calculate eigenvalues, vectors, output to STATS dataset **;

    proc iml;

     ****  Linthurst data ********;
    XPXI =   {.4865711 -.0663498 -.0001993,
             -.0663498 0.0146211  -0.0000012,
            -0.0001993 -0.0000012 .00000026} ;
    XPX = inv(XPXI);
                      MSE=160865;
                      C = 3*MSE*finv(0.95,3,42);
        STD = SQRT(VECDIAG(XPXI*MSE)); print STD;

 *** pick one of these variables to slice on ***;
        * The Z variable *;
   %if &id=1 or &id=4 %then %do;  axis = 3;
 A22 = XPX[2:3,2:3] ;
 A11 = XPX[1,1]; A12 = XPX[1,2]; A13 = XPX[1,3];
 STDERR = STD[1,1];
 %end;

        * The X variable *;
   %if &id=2 or &id = 5 %then %do; axis=1;
 A22 = XPX[{1,3}, {1,3}];
 A11 = XPX[2,2]; A12 = XPX[2,1]; A13 = XPX[2,3];
 stderr = std[2,1];
  %end;

        * The Y variable * ;
   %if &id=3 or &id = 6 %then %do; axis=2;
 A22 = XPX[1:2,1:2];
 A11 = XPX[3,3]; A12 = XPX[3,1]; A13 = XPX[3,2];
 STDERR = STD[3,1];
 %end;

   *************************************************
   ***  Your slicing variable is now d1.         ***
   ***  figure out the maximum possible d1 value ***
   ***  output all info to dataset STATS(i)      ***
   *************************************************;

   maxd1 = sqrt(C*inv(a11 - (a12||a13) * inv(A22) * (a12//a13) ));
   L = eigval(A22); T = eigvec(A22);

   data = c || a11 || a12 || a13 || a22[1, ] || a22[2,2]
             || l[1] || l[2] || t[1,1] || t[1,2] || t[2,1] || t[2,2]
             || maxd1 ||axis||stderr;
   name = {'c', 'a11', 'a12', 'a13', 'a22', 'a23', 'a33',
           'L1', 'L2', 'T11', 't12', 't21', 't22' , 'd1max', 'Axis',
            'Stderr' };


   create stats&id from data [colname=name];
   append from data;


      *****************************************************
      ********* Create the ellipses for each slice ********;
      *****************************************************;

   data a; set stats&id;
     keep d1 f1 f2 u1 u2 L1 L2  t11 t12 t21 t22  order bon scheff
        function segment;  segment = "Slices";
             *************************************
             *  First do Bonferroni calculations *
             *************************************;
     file print;
          tcrit = tinv(.99166667,42); bon = stderr*tcrit;
          scheff = sqrt(3*Finv(.95,3,42))*stderr;
          put 'Axis = ' axis '  half width is ' d1max;
            if axis=3 then put @10 'Axis 3 is Z=beta0';
            if axis=1 then put @10 'Axis 1 is X=beta2';
            if axis=2 then put @10 'Axis 2 is Y=beta1';
          put 'Bonferroni half width = ' bon ' = ' stderr ' * ' tcrit;
          put 'Scheffe half width = ' scheff;

             ***   Start slicing  *********;
     d1max=&pct*d1max;  den=&ns-1;
        * special adjustment if just 1 slice *;
        if den > 0 then step1 = 1.9999*d1max/(&ns-1);
        else step1 = d1max*1.000001;
        if den>0 then d1min=-d1max; else d1min=0;

             *** Calculate radius for the slice ***;
        do d1 = d1min to d1max by step1;
                 u1 = t11*2*d1*a12 + t12*2*d1*a13;
                 u2 = t21*2*d1*a12 + t22*2*d1*a13;
        right = c -d1*d1*a11 +u1*u1/4/L1 + u2*u2/4/L2;
        up = .99995*sqrt(right);    step = up/50;

    ****  move to a point on the ellipse ***;
    f1=up; f2 = 0; order = -1; function='move'; output;

    ****  Draw around the ellipse ****;
    do f1 = 0 to up by step; function='draw';
    f2    = .99995*sqrt(right-f1*f1);
                order = up-f1;   output;
    f1 = - f1;  order = up-f1;   output;
    f2 = - f2;  order = 3*up+f1; output;
    f1 = - f1;  order = 3*up+f1; output;
     end;
    *** close it up ****;
     f1=up; f2=0;  order = floor(4*up+2); output;
     end;

         ****************************;
         ***    rescaling         ***;
         ****************************;

   proc sort data=a; by d1 order;
   data a; set a; drop f1 f2 u1 u2;
      Z = d1;
      X = f1/sqrt(L1) - u1/2/L1;
      Y = f2/sqrt(L2) - u2/2/L2;
              output;

         ************************************;
         ***   rotating by eigenvectors   ***;
         ************************************;

   proc sort data = a; by z order ;
   data a; set a; drop t11 t12 t21 t22 d1 d2 d3 ;
     xsys='2'; ysys='2'; zsys='2';
     d2 = x*T11 + y*T21 + .000001*ranuni(1827655);
     d3 = x*T12 + y*T22 + .000001*ranuni(1827655);

     *************************************************
     *  Assigning right names (X, Y, Z)              *
     *************************************************;

     ***** pick corresponding one of these rows   *****;
     *                                                *;
     if &id=1 or &id=4 then do;
            z=d1; x=d2; y=d3;  end;      * Z variable *;
     if &id=2 or &id=5 then do;
            x=d1; y=d3; z=d2;  end;      * X variable *;
     if &id=3 or &id=6 then do;
            y=d1; z=d2; x=d3 ; end;      * Y variable *;
     *                                                *;
     **************************************************;

      *********************************
      *       Create a dataset        *
      *********************************;

  data a&id; length color $ 8;  set a;
     color = "&color";

  %mend ci;


      ********************************
      * Run the CI macro             *
      ********************************;


%ci (4,19,gray,.86);   data a4; set a4; zz=0;
 %ci(5,19,gray,.86);   data a5; set a5; zz=0;
 %ci(6,19,gray,.86);   data a6; set a6; zz=0;
   *** Use 0.6861 for Rawlings pictures.  Use 0.856271 to **
   *** slice on the Bonferroni box.                       **;
  * title h=2 'Bonferroni and Scheffe Boxes';
  * title h=1 'Rawlings slices';
    title ' ';

      *****************************************
      *  Merge ellipsoid and slice datasets   *
      *****************************************;

   data all; length function $8.;  set a4 a5 a6;


      *********************************
      *  Create corner points to plot *
      *********************************;

   proc means noprint data=all;
     output out=corner min = minx miny minz
                       max = maxx maxy maxz;
                       var x y z;
   data corner; set corner;

   ff = 1.75;      *  <--- distance to wall *;
   drop ff;

 x=minx*ff; y=miny*ff; z=minz*ff; output;
 x=maxx*ff; y=maxy*ff; z=maxz*ff; output;




     data all;    set all  end=eof; when = 'B'; size=3;


       *** Shade and resize parts outside Bon box ***;
     if x <-120.935 then do; color='green'; when='A';  size=3; end;
     if y <-   .508 then do; color='green'; when='A';  size=3; end;
     if z > 697.655 then do; color='green'; when='A';  size=3; end;
     if x > 120.95 then do; color='blue'; size=3; end;
     if y >    .51 then do; color='blue'; size=3; end;
     if z <-697.66 then do; color='blue'; size=3; end;

        **********BETAS***BETAS****BETAS*********;
            beta1 = -507;   beta2 = 412;
                    beta3 = -.4871;
        *****************************************;

     z=z + beta1; x=x + beta2; y=y + beta3;  output;


     if eof then do;
         set stats4; bonz=tinv(.9916667,42)*stderr;
          schefz = sqrt(3*Finv(.95,3,42))*stderr;
         set stats5; bonx=tinv(.9916667,42)*stderr;
          schefx = sqrt(3*Finv(.95,3,42))*stderr;
         set stats6; bony=tinv(.9916667,42)*stderr;
          schefy = sqrt(3*Finv(.95,3,42))*stderr;
         color = 'cyan';  
              ************ Bonferroni Box *************;
         z = beta1+bonz; x=beta2+bonx; y=beta3+bony; function='move';
                                                                 output;
         z = z-2*bonz; function = 'draw';                        output;
         x = x - 2*bonx; output; z=z+2*bonz; output; x=x+2*bonx; output;
         y = y -2*bony; output; z=z-2*bonz; output; y = y+2*bony;output;
         z = beta1-bonz; x=beta2-bonx; y=beta3-bony; function='move';
                                                                 output;
         z = z+2*bonz;function = 'draw';  output; x=x+2*bonx;    output;
         z=z-2*bonz; output; x=x-2*bonx;                         output;
         y=y+2*bony; output; z=z+2*bonz; output; y=y-2*bony;     output;
         z=z-2*bonz; output;

              ************ Scheffe box **************;
		 color="red";
         z = beta1+schefz; x=beta2+schefx; y=beta3+schefy; function='move';
                                                                 output;
         z = z-2*schefz; function = 'draw';                        output;
         x = x - 2*schefx; output; z=z+2*schefz; output; x=x+2*schefx; output;
         y = y -2*schefy; output; z=z-2*schefz; output; y = y+2*schefy;output;
         z = beta1-schefz; x=beta2-schefx; y=beta3-schefy; function='move';
                                                                 output;
         z = z+2*schefz;function = 'draw';  output; x=x+2*schefx;    output;
         z=z-2*schefz; output; x=x-2*schefx;                         output;
         y=y+2*schefy; output; z=z+2*schefz; output; y=y-2*schefy;   output;
         z=z-2*schefz; output;
     end;

     data corner;

        **********BETAS***BETAS****BETAS*********;
            beta1 = -507;   beta2 = 412;
                    beta3 = -.4871;
        *****************************************;

     set corner; z=z+beta1; x=x+beta2; y=y+beta3;
         label z="Beta0"; label X="Beta2"; label Y="Beta1";

        **** FINISHED********FINISHED ***********;


   proc sort data=all; by zz;

   data greek;  xsys='5'; ysys='5';
      x=12; y=60;  %greek('b','0');
      x=20; y=14;  %greek('b','3');
      x=65; y= 8;  %greek('b','2');

   data all; set all greek;

                             /*
   proc g3d data = corner annotate = all;
           scatter y*x=z/rotate=3 tilt=85 noneedle shape='point'
            grid xticknum=2 yticknum=2 zticknum=2 ;

   proc g3d data = corner annotate = all;
           scatter y*x=z/rotate=5 tilt=85 noneedle shape='point'
            grid xticknum=2 yticknum=2 zticknum=2 ;
                         */
   proc g3d data = corner annotate = all;
           scatter y*x=z/rotate=20 tilt=85 noneedle shape='point'
            /* grid */ xticknum=2 yticknum=2 zticknum=4 caxis=black;
            label x='.'; label y='.'; label z='.';
			title "Confidence regions"; 
			title2 c=red "Scheffe box outside (tangent)";
			title3 c=cyan "Bonferroni box inside"; 
			title4 "95 pct. Confidence Ellipsoid"; 

                         /*
   proc g3d data = corner annotate = all;
           scatter y*x=z/rotate=65 tilt=85 noneedle shape='point'
            grid xticknum=2 yticknum=2 zticknum=2 ;


   proc plot data=a uniform; plot x*y ='*'; *  by z;
                         */
run;