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

  /*          SAS SAMPLE LIBRARY                                  */

  /*                                                              */

  /*    NAME: betafit                                             */

  /*   TITLE: Fits the Parameters of a Simple Beta distribution   */

  /* PRODUCT: STAT                                                */

  /*  SYSTEM: All                                                 */

  /*    KEYS: nonlinear models,                                   */

  /*   PROCS: DATA NLIN MEANS                                     */

  /*    DATA:                                                     */

  /*                                                              */

  /*     REF: Charnes, Frome and Yu, J. Amer. Stat. Soc. V71,     */

  /*          p 169 (1976)                                        */

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





  /* ML ESTIMATES OF BETA PARMS USING NLIN               */

  /* FITTING THE PARAMETERS OF A BETA DISTRIBUTION USING */

  /* MAXIMUM LIKELIHOOD                                  */



  /* Generate Data from the Beta Distribution */



  DATA;

     DO I=1 TO 100; X=BETAINV(RANUNI(1),2,4); OUTPUT; END;

     RUN;



  /* Transform to log(x) and log(1-x)    */

  DATA ONE; SET;

     X1=LOG(X); X2=LOG(1-X);



  PROC MEANS NOPRINT; VAR X X1 X2;

     OUTPUT OUT=ONE SUM=X X1 X2 N=N VAR=V;



  DATA ONE;SET ONE;

   /* METHOD OF MOMENTS INTIAL ESTIMATES */

   /*  OF A AND B.                       */

   /*  mean(x) = a/(a+b), var(x) = ab/( (a+b)**2 (a+b+1) )    */

   /*   hence (a+b) = mu(1-mu)/sigma**2 - 1                   */

   /*         b     = (1-mu)(a+b)                             */



   A_INIT = (X/N)*(1-X/N)/V-1;

   IF ( A_INIT < 0 ) THEN A_INIT=1;

   B_INIT = (1-X / N) * A_INIT;

   A_INIT = A_INIT - B_INIT;

   IF ( A_INIT < 1 ) THEN A_INIT = 1;

   IF ( B_INIT < 1 ) THEN B_INIT = 1;



   /* CREATE TWO OBSERVATIONS, ONE FOR LOG(X) AND ONE FOR LOG(1-X) */

   TYPE=0; X=X1; OUTPUT;

   TYPE=1; X=X2; OUTPUT;

   KEEP X N TYPE A_INIT B_INIT;



  Title "Using NLIN to fit the parameters of the beta distribution";

  /* ------------------------------------------------------------

   *

   *   Equations are:

   *

   *       grad = 0 = du/dp sigma(-1) r

   *

   *         where u are the two means of log(x) and log(1-x)

   *         as functions of a and b, p is the parameter in question,

   *         sigma is the variance covariance matrix of

   *         (log(x), log(1-x) ) and r is the vector of residuals,

   *         ( log(x) - predicted mean, log(1-x) - predicted mean .

   *

   *       hessian = H = du/dp sigma(-1) du/dp.

   *

   *   These equations are realized in the NLIN situation by

   *   factor sigma = L*L` and using

   *     L(-1) * r

   *     L(-1) * (du/dp)

   *   as two independent observations. The effect is to

   *   reproduce the gradient and hessian equation exactly.

   *   Note that the data have been concentrated first.

   *     Sum ( log(x) ) and Sum( log(1-x) ) are sufficient

   *     statistics for this problems.

   *   Note that the resulting variance covariance matrix is

   *     the usual maximum likelihood estimate and that for this

   *     model the expected and observed information matrices are

   *     identical. The actual residual sum of squares is zero in

   *     this implementation, as are the residual degrees of freedom.

   * --------------------------------------------------------------- */



  PROC NLIN SIGSQ=1; /* TO GENERATE THE ML STANDARD ERRORS */

     PARMS A=4  B=4;

     RETAIN Z TA TB TAB WAB WAA WBB LAB LAA LBB UA UB ;



      /* 2 BY 2 CHOLESKY ROOT OF THE VARIANCE MATRIX */

     IF _ITER_ = 0 AND _OBS_=1 THEN DO;

        A=A_INIT;

        B=B_INIT;

     END;

     IF _OBS_=1 THEN DO;

        /* NEED THE TRIGAMMA FUNCTION, A NUMERICAL DERIVATIVE OF */

        /* THE DIGAMMA FUNCTION WILL SUFFICE.                    */

        TA =(DIGAMMA(A+.000001) - DIGAMMA(A-.000001))/.000002 ;

        TB =(DIGAMMA(B+.000001) - DIGAMMA(B-.000001))/.000002 ;

        TAB=(DIGAMMA(A+B+.000001) - DIGAMMA(A+B-.000001))/.000002 ;



        /* ON MVS, CMS 5.18 THERE IS A TRIGAMMA FUNCTION.

        WAA=N*(TRIGAMMA(A)-TRIGAMMA(A+B));

        WAB=N*(           -TRIGAMMA(A+B));

        WBB=N*(TRIGAMMA(B)-TRIGAMMA(A+B));

        */



        WAA=N*(TA-TAB);                      /* VARIANCE MATRIX     */

        WAB=N*(  -TAB);

        WBB=N*(TB-TAB);



        LAA=SQRT(WAA);                       /* CHOLESKY ROOT         */

        LAB=WAB/LAA;

        LBB=SQRT(WBB-LAB*LAB);



        UA = N*(DIGAMMA(A)-DIGAMMA(A+B));    /* EXPECT OF LOG(X)      */

        UB = N*(DIGAMMA(B)-DIGAMMA(A+B));    /* EXPECT OF LOG(1-X)    */

     END;



     IF TYPE=0 THEN DO; /* POOR MANS CHOLESKY ROOT APPLICATION */

        DA=WAA/LAA;     /* DERIVATIVES OF RESIDUALS            */

        DB=WAB/LAA;

        S=(X-UA)/LAA;   /* RESIDUAL: LOG(X)-EXPECTED VALUE OF LOG X */

        Z=S;            /* SAVE FOR OTHER LOG(1-X)                  */

        L=A*X-N*(LGAMMA(A)+LGAMMA(B)-LGAMMA(A+B));

     END;

     ELSE IF TYPE=1 THEN DO;         /* FOR LOG(1-X) COMPONENT */

        DA=(WAB-LAB*WAA/LAA)/LBB;

        DB=(WBB-LAB*WAB/LAA)/LBB;

        S=((X-UB)-LAB*Z/LAA)/LBB;

        L=B*X;

     END;



     /* NLIN FIT */



     MODEL S=0;                 /* S ALREADY SET TO BE RESIDUAL */

     _LOSS_  = -L;              /* LIKELIHOOD FUNCTION FOR LOSS */

     DER.A = DA;                /* DER WITH RESP. TO A AND B    */

     DER.B = DB;


     run;