maxlik.src

来自「没有说明」· SRC 代码 · 共 962 行 · 第 1/3 页

    _max_key = 1;
    _max_MaxTime = 1e5;   /* maximum time for bootstrapping */
    _max_Active = 1;
    _max_GradStep = {.};      /* size of increment for computing gradient */
    _max_GradCheckTol = 0;
    _max_Diagnostic = {.};
    _max_Alpha = .05;
    _max_Switch = 0;

    _max_BayesAlpha = 1.4;
    _max_PriorProc = {.};
    _max_Increment = 0;     /* if nonzero, histogram increments */
    _max_Center = 0;        /* if nonzero, center points for histogram */
    _max_Width = 2;         /* width of histogram = _max_Width * sd's */
    _max_NumSample = 50;   /* bootstrap sample size */
    _max_NumCat = 16;       /* # of cat's for bootstrapped histogram */
    _max_Select = {.};


endp;



/*
** PROC MAXPrt
**
** FORMAT
**           { x,f,g,vc,ret } = MAXPrt(x,f,g,vc,ret);
**
** INPUT
**           x - Kx1 vector, estimated parameters
**           f - scalar, function at minimum
**           g - Kx1 vector, gradient evaluated at x
**          vc - KxK matrix, covariance matrix of parameters
**     retcode - scalar, return code:
**
** OUTPUT
**         same as input
**
**
** GLOBALS
**
**    _max_ParNames - Kx1 character vector, parameter labels.
*/

proc (5) = maxprt(x,f,g,h,ret);
    local lbl,mask,fmt,se,ratio,oldfmt,str;
    oldfmt = formatnv("*.*lf"~8~3);
    print;
    call header("MAXLIK",_max_dsn,_ml_ver);
    print;
    print "return code = " ftos(ret,"%*.*lf",4,0);
    if ret == 0;
       print "normal convergence";
    elseif ret == 1;
       print "forced termination";
    elseif ret == 2;
       print "maximum number of iterations exceeded";
    elseif ret == 3;
       print "function calculation failed";
    elseif ret == 4;
       print "gradient calculation failed";
    elseif ret == 5;
       print "Hessian calculation failed";
    elseif ret == 6;
       print "step length calculation failed";
    elseif ret == 7;
       print "function cannot be evaluated at initial parameter values";
    elseif ret == 8;
       print "number of elements in the gradient vector inconsistent";
       print "with number of starting values";
    elseif ret == 9;
       print "gradient function returned a column vector rather than";
       print "the required row vector";
    elseif ret == 10;
       print "secant update failed";
    elseif ret == 11;
       print "maximum time exceeded";
    elseif ret == 99;
       print "termination condition unknown";
    elseif ret == 20;
       print "Hessian failed to invert";
    elseif ret == 34;
       print "data set could not be opened";
       retp(x,f,g,h,ret);
    endif;
    print;
    print "Mean log-likelihood " ftos(f,"%#*.*lG",15,6);
    print "Number of cases     " ftos(_max_NumObs,"%*.*lf",1,0);
    print;
    if ret < 3 or ret == 1;
        print "Covariance matrix of the parameters computed by the following"\
         " method:";
        if _max_IterData[3] $== "SECANT";
            print "Estimated Hessian from the secant update";
        elseif _max_IterData[3] $== "HESS";
            print "Inverse of computed Hessian";
        elseif _max_IterData[3] $== "XPROD";
            print "Cross-product of first derivatives";
        elseif _max_IterData[3] $== "QML";
            print "QML covariance matrix";
        elseif _max_IterData[3] $== "BOOT";
            print "Bootstrapped covariance matrix";
        elseif _max_IterData[3] $== "BAYES";
            print "Bayesian covariance matrix";
        else;
            print "Not computed or failed to invert";
        endif;
    endif;
    if _max_ParNames $== 0 or rows(_max_ParNames) /= rows(x);
        lbl = 0 $+ "P" $+ ftocv(seqa(1,1,rows(x)),2,0);
    else;
        lbl = _max_ParNames;
    endif;
    if rows(g) /= rows(x);
        g = miss(zeros(rows(x),1),0);
    endif;

    if not scalmiss(h) and _max_IterData[3] $/= "NOCOVP" and ret < 3;
        print;
        print "Parameters    Estimates     Std. err.  Est./s.e.  Prob.    G"\
            "radient";
        print "------------------------------------------------------------"\
            "------";
        se = real(sqrt(diag(h)));
        mask = 0~1~1~1~1~1;
        fmt = { "-*.*s"  9 8,
                "*.*lf" 14 4,
                "*.*lf" 14 4,
                "*.*lf"  9 3,
                "*.*lf"  9 4,
                "*.*lf" 12 4 };
        ratio = x./se;
        call printfm(lbl~x~se~ratio~2*cdfnc(abs(ratio))~g,mask,fmt);

        print;
        print "Correlation matrix of the parameters";
        call printfmt(h./se./se',1);
    else;
        print;
        print "The covariance matrix of the parameters failed to invert";
        print;
        print "Parameters     Estimates      Gradient";
        print "----------------------------------------------";
        mask = 0~1~1;
        let fmt[3,3] = "-*.*s" 9 8 "*.*lf" 14 6 "*.*lf" 14 6;
        call printfm(lbl~x~g,mask,fmt);
    endif;
    print;
    print "Number of iterations    " ftos(_max_IterData[1],"%*.*lf",1,0);
    print "Minutes to convergence  " ftos(_max_IterData[2],"%*.*lf",10,5);
    call formatnv(oldfmt);
    retp(x,f,g,h,ret);
endp;

/*
** PROC MAXCLPrt
**
** FORMAT
**           { x,f,g,cl,ret } = MAXCLPrt(x,f,g,cl,ret);
**
** INPUT
**           x - Kx1 vector, estimated parameters
**           f - scalar, function at minimum
**           g - Kx1 vector, gradient evaluated at x
**          cl - Kx2 matrix, confidence limits
**     retcode - scalar, return code:
**
** OUTPUT
**         same as input
**
**
** GLOBALS
**
**    _max_Alpha - scalar, (1-_max_Alpha)% confidence limits are computed
**
**    _max_ParNames - Kx1 character vector, parameter labels.
**
**
** REMARKS
**
** MAXPrt as documented in the manual can handle confidence limits as
** well as covariance matrices.  However, when there are two parameters
** MAXPrt cannot distinguish a matrix of confidence limits from a matrix
** of covariance matrices.  In general it might be best to use MAXCLPrt
** for confidence limits.
*/

proc (5) = maxclprt(x,f,g,cl,ret);
    local lbl,mask,fmt,oldfmt,str;
    oldfmt = formatnv("*.*lf"~8~3);
    print;
    call header("MAXLIK",_max_dsn,_ml_ver);
    print;
    print "return code = " ftos(ret,"%*.*lf",4,0);
    if ret == 0;
       print "normal convergence";
    elseif ret == 1;
       print "forced termination";
    elseif ret == 2;
       print "maximum number of iterations exceeded";
    elseif ret == 3;
       print "function calculation failed";
    elseif ret == 4;
       print "gradient calculation failed";
    elseif ret == 5;
       print "Hessian calculation failed";
    elseif ret == 6;
       print "step length calculation failed";
    elseif ret == 7;
       print "function cannot be evaluated at initial parameter values";
    elseif ret == 8;
       print "number of elements in the gradient vector inconsistent";
       print "with number of starting values";
    elseif ret == 9;
       print "gradient function returned a column vector rather than";
       print "the required row vector";
    elseif ret == 10;
       print "secant update failed";
    elseif ret == 11;
       print "maximum time exceeded";
    elseif ret == 99;
       print "termination condition unknown";
    elseif ret == 20;
       print "Hessian failed to invert";
    elseif ret == 34;
       print "data set could not be opened";
       retp(x,f,g,cl,ret);
    endif;
    print;
    print "Mean log-likelihood " ftos(f,"%#*.*lG",15,6);
    print "Number of cases     " ftos(_max_NumObs,"%*.*lf",1,0);
    print;
    if _max_ParNames $== 0 or rows(_max_ParNames) /= rows(x);
        lbl = 0 $+ "P" $+ ftocv(seqa(1,1,rows(x)),2,0);
    else;
        lbl = _max_ParNames;
    endif;
    if rows(g) /= rows(x);
        g = miss(zeros(rows(x),1),0);
    endif;

    if rows(cl) == rows(x);
            print;
            str = "                          " $+
                     ftos((1-_max_Alpha),"%*.*lG",8,6)$+" confidence limits";
            print str;
            print "Parameters    Estimates     Lower Limit   Upper Limit"\
                "   Gradient";
            print "---------------------------------------------------------"\
                "---------";
            mask = 0~1~1~1~1;
            fmt = { "-*.*s"  9 8,
                    "*.*lf" 14 4,
                    "*.*lf" 14 4,
                    "*.*lf" 14 4,
                    "*.*lf" 12 4 };
            call printfm(lbl~x~cl~g,mask,fmt);
    else;
        print;
        print "Confidence limits are not available";
        print;
        print "Parameters     Estimates      Gradient";
        print "----------------------------------------------";
        mask = 0~1~1;
        let fmt[3,3] = "-*.*s" 9 8 "*.*lf" 14 6 "*.*lf" 14 6;
        call printfm(lbl~x~g,mask,fmt);
    endif;
    print;
    print "Number of iterations    " ftos(_max_IterData[1],"%*.*lf",1,0);
    print "Minutes to convergence  " ftos(_max_IterData[2],"%*.*lf",10,5);
    call formatnv(oldfmt);
    retp(x,f,g,cl,ret);
endp;





/*
** PROC MAXTlimits
**
** FORMAT
**          cl = MAXTlimits(b,cov);
**
** INPUT
**          b - Kx1 vector, parameter estimates
**
**        cov - KxK matrix, covariance matrix of estimates
**
** OUTPUT
**         cl - Kx2 matrix, lower (first column) and upper (second column)
**                          (1-_max_Alpha)% confidence limits
**
** GLOBALS
**         _max_Alpha - scalar, (1-_max_Alpha)% confidence limits are computed
**
**         _max_NumObs - scalar, number of observations
**
** REMARKS
** _max_NumObs must be set.
**
** Computes (1-_max_Alpha)% confidence limits given parameter estimates and
** covariance matrix of the parameters.
*/


proc maxtlimits(b,cov);
     local dv;
     if scalmiss(_max_NumObs) or _max_NumObs == 0 or
         (_max_NumObs - rows(b)) <= 0;
         if not trapchk(4);
             errorlog "_max_NumObs not set or negative DF";
         endif;
         retp(error(0));
     endif;
     dv = cdftci(0.5*_max_Alpha,_max_NumObs-rows(b))*real(sqrt(diag(cov)));
     retp((b-dv)~(b+dv));
endp;