Resources

Model Definition for Gamma Distribution

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

                    SAS Sample Library

        Name: svrtgamm.sas
 Description: Example Program from SAS/ETS User's Guide,
              The SEVERITY Procedure
       Title: Model Definition for Gamma Distribution
     Product: SAS/ETS Software
        Keys: fitting continuous distributions
        PROC: SEVERITY
       Notes: If you run this sample program without any modification, then
              the Sasuser.Svrtdist library contains functions and subroutines
              that are identical to those in the Sashelp.Svrtdist library.

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

proc fcmp outlib=sasuser.svrtdist.models;
    function GAMMA_DESCRIPTION()$32;
        length model $32;
        model = "Gamma Distribution";
        return(model);
    endsub;

    function GAMMA_PARMCOUNT();
        return(2);
    endsub;

    function GAMMA_LOGPDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        return(logpdf("GAMMA", x, Alpha, Theta));
    endsub;
    function GAMMA_PDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        return(pdf("GAMMA", x, Alpha, Theta));
    endsub;

    function GAMMA_CDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */

        eps=3.0e-24;     /* relative accuracy  */
        fpmin=1.0e-37;   /* smallest floating point value */

        y=x/Theta;

        if ((y < 0.0) or (Alpha <= 0.0)) then do;
            return(.);
        end;
        gln = lgamma(Alpha);
        if (y < (Alpha+1.0)) then do;
            /* begin inline of function gser */
            if (y <= 0.0) then do;
                return(0.0);
            end;
            else do;
                ap=Alpha;
                sum=1.0/Alpha;
                del=sum;

                do i=1 to 1000;
                    ap=ap+1;
                    del=del*(y/ap);
                    sum=sum+del;
                    if (abs(del) < abs(sum)*eps) then do;
                        return(sum*exp(-y+Alpha*log(y)-gln));
                    end;
                end;
            end;

            return(.);
            /* end inline of function gser */
        end;
        else do;
            /* begin inline of function gcf */
            b=y+1.0-Alpha;
            c=1.0/fpmin;
            d=1.0/b;
            h=d;

            do i= 1 to 1000;
                an = -i*(i-Alpha);
                b=b+2.0;
                d=an*d+b;

                if (abs(d) < fpmin) then d=fpmin;
                c=b+an/c;
                if (abs(c) < fpmin) then c=fpmin;
                d=1.0/d;
                del=d*c;
                h=h*del;

                if (abs(del-1.0) < eps) then goto MyBreak;
            end;
        MyBreak:
            return(1.0-(exp(-y+Alpha*log(y)-gln)*h));
            /* end inline of function gcf */
        end;
    endsub;
    function GAMMA_LOGCDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        cdf = GAMMA_CDF(x, Theta, Alpha);
        if (not(missing(cdf)) and cdf > 0) then
            return (log(cdf));
        else
            return (.);
    endsub;

    function GAMMA_LOGSDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        cdf = GAMMA_CDF(x, Theta, Alpha);
        if (not(missing(cdf))) then do;
            if (cdf < 1) then
                return (log1px(-cdf));
            else
                return (.);
        end;
        else
            return (.);
    endsub;
    function GAMMA_SDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        cdf = GAMMA_CDF(x, Theta, Alpha);
        if (not(missing(cdf))) then
            return (1-cdf);
        else
            return (.);
    endsub;

    subroutine GAMMA_PARMINIT(dim, x[*], nx[*], F[*], Ftype, Theta, Alpha);
        outargs Theta, Alpha;

        m1 = 0;
        m2 = 0;
        lm = 0;
        nlm = 0;
        n = 0;
        do i=1 to dim;
            t1 = nx[i] * x[i];
            m1 = m1 + t1;
            m2 = m2 + (t1 * x[i]);
            if (x[i] > 0) then do;
                lm = lm + (nx[i] * log(x[i]));
                nlm = nlm + nx[i];
            end;
            n = n + nx[i];
        end;

        if (n = 0) then do;
            Alpha = .;
            Theta = .;
        end;
        else do;
            m1 = m1 / n;
            usemom = 1;
            if (nlm = n) then do;
                /* Compute estimates using approximate maximum likelihood */
                s = log(m1) - lm/n;
                if (s ne 0) then do;
                    s1 = (s-3)**2 + 24*s;
                    if (s1 > 0) then do;
                        Alpha = (3 - s + sqrt(s1))/(12*s);
                        if (Alpha > 0) then do;
                            Theta = m1/Alpha;
                            usemom = 0;
                        end;
                    end;
                end;
            end;

            if (usemom = 1) then do;
                /* Compute estimates using method of moments */
                m2 = m2 / n;
                t1 = m2 - m1**2;
                if (t1 < constant('MACEPS')) then do;
                    Alpha = 1;
                    Theta = m1;
                end;
                else do;
                    Alpha = m1**2 / t1;
                    Theta = t1 / m1;
                end;
            end;
        end;
    endsub;

    function GAMMA_QUANTILE(p, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        return(quantile("GAMMA", p, Alpha, Theta));
    endsub;
quit;