/*
maxit=1;
stat=zeros(maxit,5);
seed=100;
it=1;
do while it<=maxit;
it;
e=rndns(50,1,seed);
e=(e-meanc(e))./stdc(e);
stat[it,.]=loop(e);
it=it+1;
endo;
print "Rejection Rate" sumc(stat .> 2.2)'/maxit;
stat;
*/

proc(1)=loop(e);
local n,v,f,g,df,iqr,a,unit,inc,z,u,uu,nu,nuu,nwt,ndwt,wt,dwt,lim;
local nf,ndf,ng,n2,n1,u1,u2,nu1,nu2,s1,s2,w,d1,d2,dv,g1,g2,ccinv1,ccinv2;
local cc,gg,correct1a,correct2a,correct2,test1,test,test2,stat,m,correct1,i;
local m1,m2,vv,pi2,h,xx,sm,m0;
stat=zeros(1,5);
n=rows(e);
v=sortc(e,1);
m=(.05/.9)*(1-1/sqrt(n));
m0=m;
f=zeros(n,1);
nf=zeros(n,1);
g=zeros(n,1);
df=zeros(n,1);
ndf=zeros(n,1);
s1=zeros(n,1);
s2=zeros(n,1);
iqr=v[round(3*n/4),1]-v[round(n/4),1];
a=minc(stdc(e)|iqr/1.34);
h=1.06*a*n^(-0.2);
pi2=sqrt(2*pi);
u=zeros(n,1);
nu=zeros(n,1);
inc=(v[rows(v)]-v[1]+4*h)/(n-1);
z=seqa(v[1]-2*h,inc,n);
i=1;
do while i <= n;
u=(v[i]-v)/h;
nu=(-v[i]-v)/h;
f[i,1]=meanc(exp(-0.5*u^2))/(pi2*h);
nf[i,1]=meanc(exp(-0.5*nu^2))/(pi2*h);
df[i,1]=meanc((-u/h).*exp(-0.5*u^2))/(pi2*h);
ndf[i,1]=meanc((-nu/h).*exp(-0.5*nu^2))/(pi2*h);
s1[i]=sumc(e .<= v[i]);
s2[i]=sumc(-e .<= v[i]);
i=i+1;
endo;

f=f.*(f.> m)+ m * (f .<= m);
nf=nf.*(nf.> m)+ m*(nf .<= m);
g=df./f;
ng=ndf./nf;
/*
xy(v,f);
density(v,1);
*/

/*
m1=maxc(sumc(v .< -2*stdc(v))|sumc(v .> 2*stdc(v)));
m2=m1;
*/
m1=m0*n;
m2=m1;






/* n2 is the number of -ve v*/
n2=sumc(v .< 0);
n1=sumc(v.>= 0);




w=(s1-s2)/sqrt(n);       /* the untransformed test */
d1=zeros(n-1,1);
d2=zeros(n-1,1);

/* the g dW */

i=1;
do while i<=n2;
d1[i]=sumc(g.*(v .<= v[i])-ng.* (-v .<= v[i]));
i=i+1;
endo;

i=n1;
do while i<=n-1;
d2[i]=sumc(g.* (v .> v[i])-ng .* ( -v .> v[i]) );
i=i+1;
endo;

dv=(v[2:n]-v[1:n-1]);
g1=(df.*df./f);
g2=g1[2:n].*dv;

/* the inv(C) */
ccinv1=zeros(n-1,1);
ccinv2=zeros(n-1,1);
cc=zeros(n-1,1);
correct1a=zeros(n-1,1);
correct1=zeros(n-1,1);
correct2=zeros(n-1,1);
correct2a=zeros(n-1,1);
gg=(df[2:n]+df[1:n-1])/2;



i=1;
do while i<=n2;
cc[i]=sumc(g2[1:i]);
ccinv1[i]=(1/cc[i]);
correct1a[i]=gg[i].*ccinv1[i].*d1[i].*dv[i].*( f[i].> m);
i=i+1;
endo;
i=n1;
do while i<=n-1;
cc[i]=sumc(g2[i:n-1]);
ccinv2[i]=(1/cc[i]);
correct2a[i]=gg[i].*ccinv2[i].*d2[i].*dv[i].*(f[i].> m);
i=i+1;
endo;



i=1;
do while i<=n2;
correct1[i]=sumc(correct1a[i:n2] .* (f[i].> m))/sqrt(n);
i=i+1;
endo;

i=n1;
do while i<=n-1;
correct2[i]=sumc(correct2a[n1:i] .* (f[i].> m))/sqrt(n);
i=i+1;
endo;




test1=(w[m1+1:n2]+correct1[m1+1:n2]-w[n2]);
test2=(w[n1:n-1-m2]-correct2[n1:n-1-m2]-w[n2]) ;

test1=maxc(abs(test1));
test2=maxc(abs(test2));
test=maxc(test1|test2);

stat[1,1]=test;
stat[1,2]=maxc(abs(w));
stat[1,3]=test1;
stat[1,4]=test2;
stat[1,5]=(test1+test2)/2;
retp(stat);
endp;