%macro onesidedlimit
(k=, OEL=, nsum=, ytwobar=, tn=, ssybar=, sse=, conflev=, num=100000, check=);
/***************************************************************************** 
 This macro computes one-sided confidence limits for the exceedance 
 probabilities of the mean exposure (denoted by theta in the reference) 
 by implementing the algorithm 1 in the reference. 
 Let xij (j=1,...,ni, i=1,...,k,)be a sample from a lognormal population.

 Input:         k = number of individual worker 
              OEL = occupational exposure limit
             nsum = sum(ni)
              
          ytwobar = overall mean of ln(xij)
               tn = [sum(1/ni)]/k
           ssybar = between sums of squares 
              sse = within sums of squares
          conflev = confidence level of limits
              num = number of generalized variables. Default: num = 100000
            check = choice of the output (1 or 2)

 Output:
            check = 1: upper confidence limits
            check = 2: lower confidence limits

 Reference: Krishnamoorthy, K. Mathew, T.  and Ramachandran, G. (2007). 
 Upper limits for the exceedance probabilities in one-way random effects model.
 Annals of Occupational Hygiene, 51, 397-406.
*******************************************************************************/ 
options ls = 64 ps = 45 nodate nonumber;
title1 'Output of one-sided confidence limits for the  ';
title2 'exceedance probabilities of the mean exposure ';

proc iml;

tn = &tn;
ytwobar = &ytwobar;
sse = &sse;
ssybar = &ssybar;
conflev=&conflev;
check=✓
OEL =&OEL;
k = &k;
nsum = &nsum;
m = #

df1 = k-1;
df2 = nsum-k;
const1 = sqrt(ssybar/k);

seed1 = int(time());
seed2 = int(time()+12345);
seed3 = int(time()+13456);
gv=j(m,1,0); /* create a vector */

do i= 1 to m;
   z = rannor(seed1);
   v1 = 2.0*rangam(seed2,df1*0.5);
   v2 = 2.0*rangam(seed3,df2*0.5);

   gu = ytwobar + z*const1/sqrt(v1);
   gsigtsq = ssybar/v1-tn*sse/v2;
   gsigt = sqrt(max(0.0,gsigtsq));
   gsigesq = sse/v2;
   gg=log(OEL)-gu-0.5*gsigesq;

   if gsigt <= 1.0e-20 then
   	do;
    	if gg > 0.0 then gv[i] = 0.0;
		else gv[i] = 1.0;
    end;
   else
    do;
      g = gg/gsigt;
      gv[i] = 1.0-probnorm(g);
	end; 
end;


L1=int((1.0-conflev)*m);
L2=int(conflev*m);

gv0=gv;
gv[rank(gv)]=gv0; /* sort the gv[i]'s */
gvl1=gv[L1];      /*  (1-alpha) lower limit */
gvl2=gv[L2];      /*  (1-alpha) upper limit */

/* print gv, gv0; */

reset name=noname center=center;

if check=1 then
do;
  print 'm=' m[format=10.0];
  print 'ytwobar=' ytwobar[format=6.3]',' 'ntilde=' tn[format=6.3];
  print 'ssybar=' ssybar[format=6.3] ',' 'sse=' sse[format=6.3] ;
  print 'k=' k[format=4.0] ',' 'N=' nsum[format=4.0] ;
  print conflev[format=4.2] 'Confidence upper limit for theta is',
  gvl2 [format=12.8];
end;

if check=2 then
do;
  print 'm=' m[format=10.0];
  print 'ytwobar=' ytwobar[format=6.3]',' 'ntilde=' tn[format=6.3] ; 
  print 'ssybar=' ssybar[format=6.3] ',' 'sse=' sse[format=6.3] ;
  print 'k=' k[format=4.0] ',' 'N=' nsum[format=4.0] ;
  print conflev[format=4.2] 'Confidence lower limit for theta is',
  gvl2 [format=12.8];
end;

quit;

%mend onesidedlimit;

/* Examples of using the macro onesidedlimit for Table 1 and 
   Table 2 in the reference */

%onesidedlimit(k=23, OEL=1.0, nsum=34, ytwobar=-3.683, tn=0.855, ssybar=16.081, 
		      sse=2.699, conflev=0.99, num=100000, check=1); /*(Table 1)*/

%onesidedlimit(k=20, OEL=1.0, nsum=28, ytwobar=-4.087, tn=0.854, ssybar=19.681, 
		      sse=9.801, conflev=0.99, num=100000, check=1); /*(Table 2)*/ 
run;