/*program to compare F variance ratio test size with robust rank test when
random variables are Cauchy*/

n1=25;n2=25;                            @ partitions of sample@
fpnt=2.2693;                            @ F, 2.5%, 24,24@
sig=ln(fpnt)|1.96;                      @ sig is 2 by 1@
n=n1+n2;                                @ total sample size@
nrep=1000;                              @ no. of replications@

r=seqa(1,1,n)/(n+1);                    @scaled ranks 1..n,  /(n+1)@
psi=(tan(pi*(r-0.5)))^2;
psi=2*psi./(1+psi)-1;psim=psi-meanc(psi); @transformed scaled ranks@
z=ones(n1,1)|zeros(n2,1); zm=z-meanc(z);  @ 1-0 dummy @
sd=sqrt(((psim'psim)*(zm'zm))/(n-1));     @ denominator of test, fixed@

tim=hsec;                               @ store time in 1/100-ths sec@
irep=1;icount=zeros(2,1);               @ initialise loop, & count@
do until irep > nrep;                   @ start loop@
    e=rndn(n,1)./rndn(n,1);             @ N(0,1)/N(0,1) is Cauchy@
    e1=e[1:n1]; e1=e1-meanc(e1);        @ split sample@
    e2=e[n1+1:n]; e2=e2-meanc(e2);
    frat=ln(e1'e1/(e2'e2));             @find log of F for 2-tailed test@

    /*optimal rank test correlates z and rank mapped into psi*/

    wk=sortc(zm~e,2);                   @ sort on e@
    eta=(wk[.,1]'psim)/sd;              @ correlate sorted z in wk with
                                          transformed rank psi@
    tst=frat|eta;                       @ stack in tests 2 by 1 col@

    icount=icount + (abs(tst).>sig);    @ () expression gives Boolean
                                           count increment@
    irep=irep+1;                        @ increment loop counter@
endo;
"time" 0.01*(hsec-tim);; "secs    nrep" nrep;          @ output results@
"rejection % f, rank" 100*(icount')/nrep;
"Monte Carlo precision for 5% size  +/-"  196*sqrt(0.05*0.95/nrep);