Resources

Model Definition for Weibull Distribution

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

                    SAS Sample Library

        Name: svrtwbul.sas
 Description: Example Program from SAS/ETS User's Guide,
              The SEVERITY Procedure
       Title: Model Definition for Weibull 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.
              Further, if you run this sample without any modification to the
              names of functions and subroutines, then PROC FCMP writes
              warnings of the following nature to the SAS log:
              WARNING: Function  was defined in a previous package.
                       Function  as defined in the current program will
                       be used as default when the package is not specified.
              You can ignore such warnings; they appear because the functions
              defined in this sample are already defined in the input library
              Sashelp.Svrtdist.

--------------------------------------------------------------*/
proc fcmp library=sashelp.svrtdist outlib=sasuser.svrtdist.models;
    function WEIBULL_DESCRIPTION() $32;
        length model $32;
        model = "Weibull Distribution";
        return(model);
    endsub;

    function WEIBULL_PARMCOUNT();
        return(2);
    endsub;

    function WEIBULL_LOGPDF(x, Theta, Tau);
        /* Theta : Scale */
        /* Tau   : Shape */
        if (x >= constant('MACEPS')) then do;
            z = (x/Theta)**Tau;
            return (log(Tau * z / x) - z);
        end;
        else
            return (.);
    endsub;
    function WEIBULL_PDF(x, Theta, Tau);
        /* Theta : Scale */
        /* Tau   : Shape */
        if (x >= constant('MACEPS')) then do;
            z = (x/Theta)**Tau;
            return ((Tau / x) * z * exp(-z));
        end;
        else
            return (0);
    endsub;

    function WEIBULL_LOGCDF(x, Theta, Tau);
        /* Theta : Scale */
        /* Tau   : Shape */
        if (x >= constant('MACEPS')) then do;
            z = (x/Theta)**Tau;
            return (log1px(-exp(-z)));
        end;
        else
            return (.);
    endsub;
    function WEIBULL_CDF(x, Theta, Tau);
        /* Theta : Scale */
        /* Tau   : Shape */
        if (x >= constant('MACEPS')) then do;
            z = (x/Theta)**Tau;
            return (1 - exp(-z));
        end;
        else
            return (0);
    endsub;

    function WEIBULL_LOGSDF(x, Theta, Tau);
        /* Theta : Scale */
        /* Tau   : Shape */
        if (x >= constant('MACEPS')) then do;
            z = (x/Theta)**Tau;
            return (-z);
        end;
        else
            return (0);
    endsub;
    function WEIBULL_SDF(x, Theta, Tau);
        /* Theta : Scale */
        /* Tau   : Shape */
        if (x >= constant('MACEPS')) then do;
            z = (x/Theta)**Tau;
            return (exp(-z));
        end;
        else
            return (1);
    endsub;

    subroutine WEIBULL_PARMINIT(dim, x[*], nx[*], F[*], Ftype, Theta, Tau);
        outargs Theta, Tau;

        /* Compute estimates using percentile matching method */
        q1 = svrtutil_percentile(0.25, dim, x, F, Ftype); /* First quartile */
        q3 = svrtutil_percentile(0.75, dim, x, F, Ftype); /* Third quartile */
        if (missing(q1) or missing(q3)) then do;
            Theta = .;
            Tau   = .;
        end;
        else do;
            ll4 = log(log(4));
            lq3 = log(q3);
            lratio = ll4/log(log(4/3));
            ltheta = (lratio*log(q1) - lq3)/(lratio-1);

            if (ltheta < constant('LOGBIG')) then do;
                Theta = exp(ltheta);

                lq3 = lq3 - ltheta;
                if (lq3 > constant('MACEPS')) then
                    Tau = ll4/lq3;
                else
                    Tau = .;
            end;
            else do;
                Theta = .;
                Tau   = .;
            end;
        end;
    endsub;

    function WEIBULL_QUANTILE(p, Theta, Tau);
        /* Theta : Scale */
        /* Tau   : Shape */
        if (p >= 1-constant('MACEPS')) then
            return (.);
        return (Theta*(-log1px(-p))**(1/Tau));
    endsub;
quit;