OC Curve for a p Chart

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: SHWPOC                                              */
/*   TITLE: OC Curve for a p Chart                              */
/* PRODUCT: QC                                                  */
/*  SYSTEM: ALL                                                 */
/*    KEYS: Shewhart Charts, p Charts, OC Curves,               */
/*   PROCS: SHEWHART GPLOT                                      */
/*    DATA:                                                     */
/*                                                              */
/*     REF: SAS/QC Software:  Usage and Reference, Version 6,   */
/*          First Edition, Volume 1 and Volume 2                */
/*                                                              */
/****************************************************************/

data CircOne;
   input Batch Fail @@;
   datalines;
 1   7   2   6   3   6   4   9    5   2
 6  11   7   8   8   8   9   6   10  19
11   7  12   5  13   7  14   5   15   8
16  13  17   7  18  14  19  19   20   5
21   7  22   5  23   7  24   5   25  11
26   4  27   6  28   3  29  11   30   3
;


proc shewhart data=CircOne;
   where Batch ^= 10 and Batch ^= 19;
   pchart Fail*Batch / subgroupn = 500
                       nochart
                       outindex  ='Revised Limits'
                       outlimits = Faillim2;
run;


data ocpchart;
   set Faillim2;
   keep beta fraction _lclp_ _p_ _uclp_;
   nucl=_limitn_*_uclp_;
   nlcl=_limitn_*_lclp_;
   do p=0 to 500;
      fraction=p/1000;
      if nucl=floor(nucl) then
         adjust=probbnml(fraction,_limitn_,nucl) -
                probbnml(fraction,_limitn_,nucl-1);
      else adjust=0;
      if nlcl=0 then
         beta=1 - probbeta(fraction,nucl,_limitn_-nucl+1) + adjust;
      else beta=probbeta(fraction,nlcl,_limitn_-nlcl+1) -
                probbeta(fraction,nucl,_limitn_-nucl+1) +
                adjust;
      if beta >= 0.001 then output;
   end;
   call symput('lcl', put(_lclp_,5.3));
   call symput('mean',put(_p_,   5.3));
   call symput('ucl', put(_uclp_,5.3));
run;

proc sgplot data=ocpchart;
   series x=fraction y=beta / lineattrs=(thickness=2);
   xaxis values=(0 to 0.06 by 0.005) grid;
   yaxis grid;
   label fraction = 'Fraction Nonconforming'
         beta     = 'Beta';
run;