%macro sison(data=_last_,
             alpha=0.05,
             decimal=12.4);

proc iml;
 reset nocenter;
 reset noname;
 use &data;
 setin &data;
 read all into x;
 data={&data};
 alpha={&alpha};
 decimal={&decimal};
 n=sum(x);
 k=nrow(x);
 p=x/n;
print 'Data set used:         ' data;
print 'Alpha =         ' alpha;
print 'Number of cells=' k;
print 'N =             ' n;
print 'Observed cell counts ',x;

start moments;
   a=lambda+c;
   b=lambda-c;
   if b<0 then b=0;
   if b>0 then den=poisson(lambda,a)-poisson(lambda,b-1);
   if b=0 then den=poisson(lambda,a);
   mu=j(4,1,0); mom=j(5,1,0);
   do r = 1 to 4;
      poisA=0;
      poisB=0;
      if a-r>=0 then poisA=poisson(lambda,a)-poisson(lambda,a-r);
      if a-r< 0 then poisA=poisson(lambda,a);
      if b-r-1>=0 then poisB=poisson(lambda,b-1)-poisson(lambda,b-r-1);
      if b-r-1<0 && b-1>=0 then poisB=poisson(lambda,b-1);
      if b-r-1<0 && b-1<0 then poisB=0;
      mu[r]=(lambda**r)*(1-(poisA-poisB)/den);
   end;
   mom[1]=mu[1];
   mom[2]=mu[2]+mu[1]-mu[1]**2;
   mom[3]=mu[3]+mu[2]*(3-3*mu[1])+(mu[1]-3*mu[1]**2+2*mu[1]**3);
   mom[4]=mu[4]+mu[3]*(6-4*mu[1])+mu[2]*(7-12*mu[1]+6*mu[1]**2)
        +mu[1]-4*mu[1]**2+6*mu[1]**3-3*mu[1]**4;
   mom[5]=den;

finish;

start truncpoi;
   m=j(k,5,0);
   do i=1 to k;
    lambda=x[i];
    run moments;
    do j=1 to 5;
     m[i,j]=mom[j];
    end;
   end;
   s1=m[+,1];
   s2=m[+,2];
   s3=m[+,3];
   do i=1 to k;
    m[i,4]=m[i,4]-3*m[i,2]**2;
   end;
   s4=m[+,4];
   probn=1/(poisson(n,n)-poisson(n,n-1));
   z=(n-s1)/sqrt(s2);
   g1=s3/(s2**(3/2));
   g2=s4/(s2**2);
   poly=1+g1*(z**3-3*z)/6+g2*(z**4-6*z**2+3)/24
         +g1**2*(z**6-15*z**4+45*z**2-15)/72;
   f=poly*exp(-z**2/2)/(sqrt(2)*gamma(0.5));
   probx=1;
   do i=1 to k;
    probx=probx*m[i,5];
   end;
   p=probn*probx*f/sqrt(s2);
finish;

start;
  pold=0;
  do c=1 to n;
   run truncpoi;
   if p > 1-alpha && pold < 1-alpha then goto done;
   pold=p;
  end;

 done:
  delta=(1-alpha-pold)/(p-pold);
  out=j(k,5,0);
  num=j(k,1,0);
  c=c-1;
  vol1=1;
  vol2=1;
  do i=1 to k;
   num[i,1]=i;
   obsp=x[i]/n;
   out[i,1]=obsp;
   out[i,2]=obsp-c/n;
   out[i,3]=obsp+c/n+2*delta/n;
   if out[i,2]<0 then out[i,2]=0;
   if out[i,3]>1 then out[i,3]=1;
   out[i,4]=obsp-c/n-1/n;
   out[i,5]=obsp+c/n+1/n;
   if out[i,2]<0 then out[i,2]=0;
   if out[i,3]>1 then out[i,3]=1;
   vol1=vol1*(out[i,3]-out[i,2]);
   vol2=vol2*(out[i,5]-out[i,4]);
  end;
  c1={'PROPORTION', 'LOWER(SG)','UPPER(SG)','LOWER(C+1)','UPPER(C+1)'};
  cov=100*(1-alpha);
  sg=(x+delta)/n;
  c2={'SG-midpoint'};
  print '-------------------------------------------------------------';
  print cov'% SIMULTANEOUS CONFIDENCE INTERVALS';
  print  '       BASED ON THE METHODS OF SISON AND GLAZ';
  print '-------------------------------------------------------------';
  print 'C = ' c;
  print 'P(c+1) = ' p(|format=5.4|);
  print 'P(c)   = ' pold(|format=5.4|);
  print 'delta =  ' delta(|format=5.4|);
  print 'Volume(SG) = ' vol1;
  print 'Volume(C+1)= ' vol2;
  print num(|format=3.0|) out(|colname=c1 format=&decimal|);
  print num(|format=3.0|) sg(|colname=c2 format=&decimal|);
finish;
run;
quit;

%mend;

quit;
run;
data one;
 input x @@;
cards;
56 72 73 59 62 87 58

run;
%sison(data=one, alpha=0.05);
run;
quit;