/**********************************************************************/


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; 


/**********************************************************************/