PROC IML;


    start KL_alpha(C,G,T);
       alphaC=1-probchi(T,1);
       mu=(C+2*1*G*T)*alphaC;
       alphaF=1-exp(-mu);
       return(alphaF);
    finish KL_alpha;

    start kl_cv(C,G,alpha);
       Tmin=0;
       al_min=1;
       Tmax=1;
       al_max=kl_alpha(C,G,Tmax);
       do while(al_max > alpha);
           tmin=tmax;
           al_min=al_max;
           tmax=2*tmax;
           al_max=kl_alpha(C,G,Tmax);
           end;
       T=Tmax;
       do while(al_min-al_max >1.e-4);
           T=(Tmax+tmin)/2;
           al=kl_alpha(C,G,T);
           if al < alpha then do;
               tmax=T;
               al_max=al;
               end;
           else do;
               tmin=T;
               al_min=al;
               end;
           end;
       return(t);
    finish kl_cv;
    store module=(kl_cv kl_alpha);
quit;


Proc IML;
     load _all_;
    alpha=.05;
    df=1;
    QTLvar=.1;
    lam_mult=QTLvar/2;
    cv=cinv(1-alpha,df);
    c=4;
    g=600;
    cv2=kl_cv(c,g,alpha);
    cv3=kl_cv(1,150,alpha);
    create power var {n lambda power power2 power3 cv};
    nmin=0;
    nmin2=0;
	nmin3=0;
    do n=10 to 1000 by 10;
        lambda=n*lam_mult;
        power=1-probchi(cv,df,lambda);
        power2=1-probchi(cv2,df,lambda);
        power3=1-probchi(cv3,df,lambda);
        if power > .9 & nmin=0 then do;
            nmin=n;
            print n power;
            end;
        if power2 > .9 & nmin2=0 then do;
            nmin2=n;
            print n power2;
            end;
        if power3 > .9 & nmin3=0 then do;
            nmin3=n;
            print n power3;
            end;
        power=100*power;
        power2=100*power2;
        power3=100*power3;
        append;
        end;
    quit;
    run;

symbol1 i=join line=1 value=none color=black ;
symbol2 i=join line=3 value=none color=black;
symbol3 i=join line=5 value=none color=black;
axis1 label=(angle=90 'Power' ) order=(0 to 100 by 10);
axis2 label=('Sampls Size');
goption device=psepsf  hsize=7in vsize=4in gsfmode=replace;

legend1 label=("Error Rate") Value=("Comp." "Chrom." "Genome");
filename grafout 'bc_power.eps';
goption gsfname=grafout;
proc gplot data=power;
    plot power*n power3*n power2*n/overlay haxis=axis2 vaxis=axis1 vref=(90) legend=legend1;
run;