cml.src

GAUSS软件的CML模块

** CMLPrt
**
** FORMAT
**           { x,f,g,vc,ret } = CMLPrt(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
**
**    _cml_ParNames - Kx1 character vector, parameter labels.
*/

proc (5) = CMLprt(x,f,g,h,ret);
    local lbl,mask,fmt,ratio,se,oldfmt;
    oldfmt = formatnv("*.*lf"~8~3);
    print;
    call header("CML",_cml_dsn,_cml_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 "error with initial gradient";
    elseif ret == 9;
        print "error with constraints";
    elseif ret == 10;
        print "secant update failed";
    elseif ret == 11;
        print "maximum time exceeded";
    elseif ret == 12;
        print "error with weights";
    elseif ret == 13;
        print "quadratic program failed";
    elseif ret == 14;
        print "equality constraint Jacobian failed";
    elseif ret == 15;
        print "inequality constraint Jacobian failed";
    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(_cml_NumObs,"%*.*lf",1,0);
    print;

    if _cml_ParNames $== 0 or rows(_cml_ParNames) /= rows(x);
        lbl = 0 $+ "P" $+ ftocv(seqa(1,1,rows(x)),2,0);
    else;
        lbl = _cml_ParNames;
    endif;

    if ret < 3 or ret == 11;
        print "Covariance of the parameters computed by the following method:";
        if _cml_IterData[3] $== "SECANT";
            print "Estimated Hessian from the secant update";
        elseif _cml_IterData[3] $== "HESS";
            print "Inverse of computed Hessian";
        elseif _cml_IterData[3] $== "XPROD";
            print "Cross-product of first derivatives";
        elseif _cml_IterData[3] $== "HETCON";
            print "QML covariance matrix";
        elseif _cml_IterData[3] $== "QML";
            print "QML covariance matrix";
        elseif _cml_IterData[3] $== "BOOT";
            print "Bootstrapped covariance matrix";
        elseif _cml_IterData[3] $== "BAYES";
            print "Bayesian covariance matrix ";
        else;
            print "Not computed or failed to invert";
        endif;
    endif;
    if rows(g) /= rows(x);
         g = miss(zeros(rows(x),1),0);
    endif;
    if not scalmiss(h) and _cml_IterData[3] $/= "NOCOVP" and ret < 3;

        if scalmiss(_cml_A) and scalmiss(_cml_C) and scalmiss(_cml_EqProc)
            and scalmiss(_cml_IneqProc) and scalmiss(_cml_Bounds);

            print;
            print "Parameters    Estimates     Std. err.  Est./s.e.  Prob. "\
                "   Gradient";
            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 "Parameters    Estimates     Std. err.    Gradient";
            print "--------------------------------------------------------"\
                "----------";
            se = real(sqrt(diag(h)));
            mask = 0~1~1~1;
            fmt = { "-*.*s" 9 8, "*.*lf" 14 4, "*.*lf" 14 4, "*.*lf" 12 4 };
            call printfm(lbl~x~se~g,mask,fmt);
        endif;
    else;
        print;
        print "The covariance 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(_cml_IterData[1],"%*.*lf",1,0);
    print "Minutes to convergence  " ftos(_cml_IterData[2],"%*.*lf",10,5);
    call formatnv(oldfmt);

    retp(x,f,g,h,ret);
endp;

/*
** CMLCLPrt
**
** FORMAT
**           { x,f,g,cl,ret } = CMLCLPrt(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
**
**    _cml_Alpha - scalar, confidence level
**
**    _cml_ParNames - Kx1 character vector, parameter labels.
**
**
** REMARKS
**
** CMLPrt as documented in the manual can handle confidence limits as
** well as covariance matrices.  However, when there are two parameters
** CMLPrt cannot distinguish a matrix of confidence limits from a matrix
** of covariance matrices.  In general it might be best to use CMLCLPrt
** for confidence limits.
*/

proc (5) = CMLCLprt(x,f,g,cl,ret);
    local lbl,mask,fmt,oldfmt,str;
    oldfmt = formatnv("*.*lf"~8~3);
    print;
    call header("CML",_cml_dsn,_cml_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 "error with initial gradient";
    elseif ret == 9;
        print "error with constraints";
    elseif ret == 10;
        print "secant update failed";
    elseif ret == 11;
        print "maximum time exceeded";
    elseif ret == 12;
        print "error with weights";
    elseif ret == 13;
        print "quadratic program failed";
    elseif ret == 14;
        print "equality constraint Jacobian failed";
    elseif ret == 15;
        print "inequality constraint Jacobian failed";
    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(_cml_NumObs,"%*.*lf",1,0);
    print;

    if _cml_ParNames $== 0 or rows(_cml_ParNames) /= rows(x);
        lbl = 0 $+ "P" $+ ftocv(seqa(1,1,rows(x)),2,0);
    else;
        lbl = _cml_ParNames;
    endif;

    if rows(lbl) /= rows(x);
        if not trapchk(4);
            errorlog "ERROR:  _cml_Parnames not conformable to coefficient "\
                "vector";
        endif;
        retp(x,f,g,cl,ret);
    endif;

    if rows(g) /= rows(x);
        g = miss(zeros(rows(x),1),0);
    endif;
    if rows(cl) == rows(x);

        print;
        str = "                          " $+
             ftos((1-_cml_Alpha),"%*.*lG",8,6)$+" confidence limits";
        print str;
        print "Parameters    Estimates     Lower Limit   Upper Limit   Grad"\
            "ient";
        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 "The covariance 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(_cml_IterData[1],"%*.*lf",1,0);
    print "Minutes to convergence  " ftos(_cml_IterData[2],"%*.*lf",10,5);
    call formatnv(oldfmt);

    retp(x,f,g,cl,ret);
endp;

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

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