/* It will take some time. Wait... */
new; cls;

p=pureseqa(0,0.01,101);
upperbound=1;
lowerbound=0;
alpha=1/2;
beta=4;
fn f(x)=cdfbeta(x,alpha,beta);    /* Suppose cdf is Beta(exactly bounded between 0 and 1). */

library pgraph;
graphset;
xy(searchcdfi(p,lowerbound,upperbound,&f),p);


/*
**  searchcdfi.txt - Calculation of Inverse CDF from any CDF function by iterated line search
**                   with a direction and bounds.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**
**  Purpose:    Computes inverse CDF of any distribution between lower bound and upper bound.
**
**  Format:     x=searchcdfi(p,lowerbound,upperbound,&f)
**
**  Input:      p           n x 1 vector, probabilities ( 0 <= p <= 1 )
**
**              lowerbound  scalar, lower bound when p=0
**
**              upperbound  scalar, upper bound when p=1
**
**              &f          a procedure, cdf function pre-setting as fn f(x)=cdf(x,para1,para2,...); 
**
**  Output:     x           n x 1 vector, values
**
**
**  Notice:     Set 'f(x)=cdfname(x,para1,para2,....);' before using this procedure by regarding
**              parameter para1, para2... in original cdf function procedure as global variables.
**
**              Parameter lowerbound and upperbound should be set properly not to blow up the 
**              original cdf function procedure.
**
**              It may take a long time to run this procedure, but the result should be accurate
**              if you set both lower and upper bounds properly beacuse CDF is monotone.
*/

proc searchcdfi(p,lowerbound,upperbound,&f);
   local Max_Iter,Tolerance,Points,n,j,i,x,xi,g,step,f:proc;
   Max_Iter=10;
   Tolerance=1e-8;
   Points=10000;                    /* Maximum for Light version of GAUSS */
   n=rows(p);
   x=zeros(n,1);
   j=1;
   do while j<=n;
      if p[j]==0;
         x[j]=lowerbound;
      elseif p[j]==1;
         x[j]=upperbound;
      elseif p[j]<0 or p[j]>1;
         x[j]=miss(0,0);
      else;
         step=(upperbound-lowerbound)/(Points-1);
         xi=seqa(lowerbound,step,Points); xi[Points]=upperbound;
         i=1;
         do while i<=Max_Iter;
            g=abs(f(xi)-p[j]);
            x[j]=xi[minindc(g)];
            if abs(f(x[j])-p[j])<Tolerance;
               break;
            elseif f(x[j])-p[j]>0;
               step=step/100;
               xi=seqa(x[j],-step,Points);
               if xi[Points]<lowerbound;
                  step=(x[j]-lowerbound)/(Points-1);
                  xi=seqa(x[j],-step,Points); xi[Points]=lowerbound;
               endif;
            else;
               step=step/100;
               xi=seqa(x[j],step,Points); 
               if xi[Points]>upperbound;
                  step=(upperbound-x[j])/(Points-1);
                  xi=seqa(x[j],step,Points); xi[Points]=upperbound;
               endif;
            endif;
            i=i+1;
         endo;
      endif;
      j=j+1;
   endo;
   retp(x);
endp;

proc phi(x);
   local y,a1,a2,a3,a4,a5,n,xx,i;
   n=rows(x);
   xx=zeros(n,1);
   i=1;
   do while i<=n;
      if x[i]>=0;
         y=1./(1+0.2316419*x[i]);
         a1=0.319381530;
         a2=-0.356563782;
         a3=1.781477937;
         a4=-1.821255978;
         a5=1.330274429;
         xx[i]=1-1/sqrt(2*pi).*exp(-x[i]^2/2).*(a1*y+a2*y^2+a3*y^3+a4*y^4+a5*y^5);
      else;
         x[i]=-x[i];
         y=1./(1+0.2316419*x[i]);
         a1=0.319381530;
         a2=-0.356563782;
         a3=1.781477937;
         a4=-1.821255978;
         a5=1.330274429;
         xx[i]=1-(1-1/sqrt(2*pi).*exp(-x[i]^2/2).*(a1*y+a2*y^2+a3*y^3+a4*y^4+a5*y^5));
      endif;
      i=i+1;
   endo;
   retp(xx);
endp;

proc pureseqa(ini,step,n);
   local x,i;
   x=zeros(n,1);
   i=1;
   do while i<=n;
      x[i]=ini+step*(i-1);
      i=i+1;
   endo;
   retp(x);
endp;