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 */
        minVal = 2.220446E-16;      /* alternatives:
                                       MACEPS     = 2.220446E-16
                                       sqrt(SMALL)= 0.1491668147e-153 */
        if (x < minVal) then do;
            x1 = minVal;
            /* assume x1/Theta~0, because it is too small */
            lp = (Alpha-1)*log(x1) - lgamma(Alpha) - Alpha*log(Theta);
        end;
        else
            lp = logpdf("GAMMA", x, Alpha, Theta);
        return(lp);
    endsub;
    function GAMMA_PDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        minVal = 2.220446E-16;      /* alternatives:
                                       MACEPS     = 2.220446E-16
                                       sqrt(SMALL)= 0.1491668147e-153 */
        if (x < minVal) then do;
            x1 = minVal;
            /* assume exp(-x1/Theta)~1, because x1/Theta is too small */
            p = x1**(Alpha-1) / (gamma(Alpha) * (Theta**Alpha));
        end;
        else
            p = pdf("GAMMA", x, Alpha, Theta);
        return(p);
    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 (missing(cdf)) then
            return (.);

        if (cdf > 0) then
            return (log(cdf));
        else
            return (-180.218266); /* constant('LOGSMALL') */
    endsub;

    function GAMMA_LOGSDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        cdf = GAMMA_CDF(x, Theta, Alpha);
        if (missing(cdf)) then
            return (.);

        if (cdf < 1) then
            return (log1px(-cdf));
        else
            return (-180.218266); /* constant('LOGSMALL') */
    endsub;
    function GAMMA_SDF(x, Theta, Alpha);
        /* Theta : Scale */
        /* Alpha : Shape */
        cdf = GAMMA_CDF(x, Theta, Alpha);
        if (missing(cdf)) then
            return (.);

        return (1-cdf);
    endsub;

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

        m1 = 0;
        m2 = 0;
        lm = 0;
        n = 0;
        minVal = 2.220446E-16;      /* alternatives:
                                       MACEPS     = 2.220446E-16
                                       sqrt(SMALL)= 0.1491668147e-153 */
        do i=1 to dim;
            x1 = max(x[i], minVal);
            t1 = nx[i] * x1;
            m1 = m1 + t1;
            m2 = m2 + (t1 * x1);
            lm = lm + (nx[i] * log(x1));
            n = n + nx[i];
        end;

        if (n = 0) then do;
            Alpha = .;
            Theta = .;
        end;
        else do;
            m1 = m1 / n;
            usemom = 1;
            /* 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;

            if (usemom = 1) then do;
                /* Compute estimates using method of moments */
                m2 = m2 / n;
                t1 = m2 - m1**2;
                if (t1 < 2.220446E-16) then do; /* constant('MACEPS') */
                    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;

    function GAMMA_MEAN(x, Theta, Alpha);
        return (Theta*Alpha);
    endsub;
quit;