/**********************************************************************/
probit(0.1)=-1.2815515655446
/**********************************************************************/
probit(1.e-10)=-6.361340902404
/**********************************************************************/
p=probt(t,df,nc);
if p eq . then do;
nrm=(t*(1-.25/df)-nc)/sqrt(1+.5*t*t/df);
p=probnorm(nrm);
end;
/**********************************************************************/
data one(keep=p df nc tinv approx);
/* given p,df,nc find tinv(p,df,nc) using the */
/* normal approximation for T-distribution */
/* n(x)=p */
/* where x=(t*(1.-.25/df)-nc)/sqrt(1+.5*t*t/df) */
input p df nc;
x=probit(p);
a=(1.-.25/df);
b=-2*nc*a;
c=nc*nc-x*x;
a=a*a-.5*x*x/df;
/* solve a*t*t+b*t+c=0; for t */
/* this gives the normal approximation */
approx=(-b+sqrt(b*b-4*a*c))/(2*a);
tinv=tinv(p,df,nc);
cards;
.90 10 50
.90 100 50
.90 1000 50
.90 100 10
.90 100 50
.90 100 100
.02 100 100
.05 10 100
proc print data=one;
run;
/**********************************************************************/
data;
do nc=1700 to 2000 by 100;
f=60;
df1=30;
df2=10;
p=probf(60,30,10,nc);
put nc= p= e23.15;
end;
run;
NC=1700 P=4.7604885665372000E-01
NC=1800 P=4.2701717860324000E-01
NOTE: Invalid argument to function PROBF .....
NC=1900 P=.
NOTE: Invalid argument to function PROBF .....
NC=2000 P=.
/**********************************************************************/
data two(keep=nc p p1 p2);
do nc=1700 to 2000 by 20;
f=60;
df1=30;
df2=10;
p=probf(f,df1,df2,nc);
/* Print the results if the correct */
/* value P can be computed */
if p<>. then do;
x1=(f-df2/df1*(df1+nc)/(df2-2))/
(df2/df1*sqrt(2/(df2-2)/(df2-4)*((df1+nc)**2
/(df2-2)+df1+2*nc)));
x2=((df1*f/(df1+nc))**(1./3.)*(1-2/9/df2)
-(1-2/9*(df1+2*nc)/(df1+nc)**2))/
sqrt(2/9*(df1+2*nc)/(df1+nc)**2+2/9/df2*
(df1/(df1+nc)*f)**(2/3));
p1=probnorm(x1);
p2=probnorm(x2);
output;
end;
end;
run;
proc print data=two;
run;
/**********************************************************************/