/*********************************************************/
/*   TITLE: Probability Plotting                         */
/* PRODUCT: GRAPH                                        */
/*  SYSTEM: ALL                                          */
/*    KEYS: NORMAL PROBABILITY PLOT                      */
/*   PROCS: GPLOT                                        */
/*                                                       */
/*    MISC:                                              */
/*                                                       */
/*********************************************************/

goptions reset=all;
goptions nocharacters nocell notext82;

* set up desired range for the percentile axis ;
data dummy;
  do p = .01, .05, .10, .50, .75, .90, .95, .99;
     prob=probit(p);
     label=p*100;
     output;
  end;
run;

*create annotate data set;
data axis;
  set dummy;
  length function color text $8;
  retain ysys '2' xsys '1' color 'blue';
  y=prob;

  * draw the reference line;
  function='move'; x=0; output;
  function='move'; x=100; size=1.5; output;

  *create tickmark labels;
  text=left(trim(label));

  * select desired tickmark labels;
  if text='1'   |  text='5'   |
     text='10'  |  text='50'  |  text='75' |
     text='90'  |  text='95'  |  text='99'
  then do;
     function='label';
     text=right(put(text,2.))||'-';
     x=0;
     position='4';
     size=1;
     output;
  end;
run;

* generating a test data set, RODS;
 data rods;
    input diameter @@;
    label diameter='Diameter in mm';
    cards;
 2.091  3.612  0.001  2.123  3.694
 1.623  1.561  0.563  1.452  3.562
 4.544  2.345  4.561  0.452  2.822
 5.501  4.251  5.404  7.366  5.445
 6.138  7.126  2.782  5.127  3.436
 0.576  1.607  4.200  3.977  5.177
 4.332  1.399  0.661  4.512  1.252
 1.404  3.739  3.525  2.160  0.410
 2.823  2.376  1.202  1.470  3.410
 3.394  0.146  4.244  0.309  1.480
 1.388  2.399  1.360  3.368  2.394
 3.248  5.409  5.304  2.239  1.781
 ;

* sort your data;
proc sort data=rods;
  by diameter;
run;

* compute the normal quantiles;
data rods;
  set rods nobs=n;
  y=(_n_-(3/8)) / (n+(1/4));
  prob=probit(y);
run;

* combine data sets;
data all;
  set dummy rods;
run;

* define and enhance axes;
axis1 value=none major=none minor=none
      label=(justify=c h=2 "     "
             justify=c h=1 c=black a=-90 r=90 "Normal Percentiles");
axis2 minor=none;

* generate probability plot;
proc gplot data=all anno=axis;
  plot prob*diameter / vaxis=axis1
                       haxis=axis2
  ;
run;

title;

/* Here is the same plot in SAS/QC, if available */
/*
proc capability data=rods graphics noprint;
  probplot diameter / normal
                      rotate
  ;
run;
*/