cmlcount.src

GAUSS软件的CML模块

endp;

/* _cmlc_CLNG - log of the gamma function of x.
**
*/

proc _cmlc_lng(x);
    local x2, x3, c, ys, yp, mask, brkpnt;
    brkpnt = 1;     /* for x below this, use direct computation */
    if not (x > 0);        /* some negative x's */
        errorlog "ERROR: Arguments must all be non-negative.";
        end;
    endif;
    brkpnt = brkpnt + 1;
    if (x > brkpnt);        /* all x's greater than this -- use stirling's
                            :: approx.
                            */
        gosub stirling;
        retp( ys );
    elseif (x <= brkpnt);           /* all x's less than this -- use
                                    :: ln(gamma)
                                    */
        gosub lngamma;
        retp( yp );
    else;           /* some x's of each type -- combine estimates */
        gosub stirling;
        gosub lngamma;
        mask = (x .> brkpnt);
        retp( ys.*mask + yp.*(.not mask) );
    endif;
    /* ----------------- subroutines follow -----------------------  */
stirling:

    /* stirling's approx */
    x2 = x.*x;
    x3 = x.*x2;
    ys = 0.5 * ln(2*pi) + (x - 0.5) .* ln(x) - x + 1./(12.*x) -
        1./(360.*x3) + 1./(1260.*x3.*x2) - 1./(1680.*x3.*x2.*x2) +
        1./(1188.*x3.*x2.*x2.*x2);
    return;
lngamma:

    yp = ln( gamma(x) );
    return;
endp;

/*
** _cmlc_ftosm - matrix field to string
**
**  y = _cmlc_ftosm(sym,n);
**
**  input:  sym = a vector of symbols
**          n = number of characters to use in each symbol
**
**  output: y = a character vector of numbers
*/

proc _cmlc_ftosm(sym,n);
    local i,res;
    res = zeros(rows(sym),1);
    i = 1;
    do while i<=rows(sym);
        res[i] = ftos(sym[i],"*.*lf",n,0);
        i = i+1;
    endo;
    retp(res);
endp;

/*
**   CMLCOUNTPRT -- Print Parameters and Estimates
**
**     { b,vc,logl } = CMLCOUNTPRT(b,vc,logl);
**
** INPUT:
**      b = kx1 coefficient vector
**      vc = kxk variance-covariance matrix
**      logl = log-likelihood
**
** OUTPUT:
**      b = kx1 coefficient vector
**      vc = kxk variance-covariance matrix
**      logl = log-likelihood
**
** GLOBALS:
**  _cmlc_Fix = variable name or column number of var with param
**             constrained to 1.
**  __title = a string with a message to be printed at the top of the output.
**
*/

proc(3) = cmlcountprt(b,vc,logl);
    local x,mask,fmt,hdr,vbls,se,t,rse,iout;
    local dataset,dep,vars;
    print;
    dataset = _cmlc_fn;
    dep = _cmlc_dp;
    vars = _cmlc_vr;

    if type(_cmlc_fn) == 13;
        call header("CMLCOUNT",_cmlc_fn,_cmlc_ver);
    else;
        call header("CMLCOUNT","",_cmlc_ver);
    endif;

    if _cmlc_Inference $== "MAXLIK" or _cmlc_Inference $== "PROFILE";
        rse = real(sqrt(diag(vc)));
        se = real(sqrt(diag(_cml_HessCov)));
        if scalmiss(se);
            errorlog "ERROR:  Covariance matrix of parameters not invertible";
            se = miss(zeros(rows(b),1),0);
        endif;
        if scalmiss(rse);
            errorlog "ERROR:  H-S Covariance matrix of parameters not inver"\
                "tible";
            rse = miss(zeros(rows(b),1),0);
        endif;
        mask = zeros(rows(b),1)~ones(rows(b),3);
        fmt = ("*.*s"~ 10~ 8)| ("*.*lf"~ 10~ _cmlc_Precision)| ("*.*lf" ~
            10~ _cmlc_Precision)| ("*.*lf"~ 10~ _cmlc_Precision);
        if type(dataset)/=13;
            dataset = "Matrix Input";
            if rows(dep)==1;
                dep = "Col."$+ftos(dep,"*.*lf",2,0);
            else;
                dep = "Col."$+ftos(dep[1],"*.*lf",2,0)$+" & " $+ "Col." $+
                    ftos(dep[2],"*.*lf",2,0);
            endif;
            t = ones(rows(vars)-1,1);
            vars[2:rows(vars)] = t.*vars[2:rows(vars)];
        else;
            if rows(dep)==1;
                dep = "" $+ lower(dep);
            else;
                dep = "" $+ lower(dep[1]) $+ " & " $+ lower(dep[2]);
            endif;
        endif;
        vars = lower(vars);
        x = vars~b~se~rse;
        print;
        print "Dependent Variable: "$+dep;
        print;
        print " Parameter  Estimate     S.E.   Het-con S.E.";
        print;
        call printfm(x,mask,fmt);
        print;
        print "log-likelihood = " ftos(logl,"-*.*lf",10,_cmlc_Precision) " ";;
        print " n = " ftos(_cml_NumObs,"-*.*lf",10,0);
        if _cmlc_Fix/=0;
            if round(_cmlc_Fix) == _cmlc_Fix and _cmlc_Fix >= 1 and
                _cmlc_Fix < 131072;
                print ("Coefficient set to 1.0:    " $+ "Col." $+
                    ftos(_cmlc_Fix,"*.*lf ",3,0) $+ _cmlc_Fix);
            else;
                print ("Coefficient set to 1.0:     " $+ _cmlc_Fix);
            endif;
        endif;

    elseif _cmlc_Inference $== "BOOT" or _cmlc_Inference $== "BAYES";

        se = real(sqrt(diag(vc)));
        if scalmiss(se);
            errorlog "ERROR:  Covariance matrix of parameters failed";
            se = miss(zeros(rows(b),1),0);
        endif;
        mask = (zeros(rows(b),1)~ones(rows(b),2));
        fmt = ones(3,1).*("*.*lf"~ 10~ _cmlc_Precision);
        if type(dataset)/=13;
            dataset = "Matrix Input";
            if rows(dep)==1;
                dep = "Col."$+ftos(dep,"*.*lf",2,0);
            else;
                dep = "Col."$+ftos(dep[1],"*.*lf",2,0)$+" & " $+ "Col." $+
                    ftos(dep[2],"*.*lf",2,0);
            endif;
            t = ones(rows(vars)-1,1);
            vars[2:rows(vars)] = t.*vars[2:rows(vars)];
        else;
            if rows(dep)==1;
                dep = "" $+ lower(dep);
            else;
                dep = "" $+ lower(dep[1]) $+ " & " $+ lower(dep[2]);
            endif;
        endif;
        vars = lower(vars);
        x = (vars~b~se);
        print;
        print "Dependent Variable: "$+dep;
        print;
        print "Estimates are means of simulated bootstrapped distributions";
        print;
        print " Parameter  Estimate   Boot S.E.";
        call printfm(x,mask,fmt);
        print;
        hdr = "bootstrap log-likelihood = " $+ ftos(logl,"-*.*lf",10,
            _cmlc_Precision) $+ " ";
        print hdr;
        print " n = " ftos(_cml_NumObs,"-*.*lf",10,0);
        if _cmlc_Fix /= 0;
            if round(_cmlc_Fix) == _cmlc_Fix and _cmlc_Fix >= 1 and
                _cmlc_Fix < 131072;
                print ("Coefficient set to 1.0:    " $+ "Col." $+
                    ftos(_cmlc_Fix,"*.*lf ",3,0) $+ _cmlc_Fix);
            else;
                print ("Coefficient set to 1.0:     " $+ _cmlc_Fix);
            endif;
        endif;
    endif;
    __title = "";

    ndpclex;
    retp(b,vc,logl);
endp;

/*
**  CMLCOUNTCLPrt -- Print Parameters and Estimates
**
**     { b,cl,logl } = CMLCOUNTCLPrt(b,cl,logl);
**
** INPUT:
**      b = kx1 coefficient vector
**      cl = kx2 confidence limits
**      logl = log-likelihood
**
** OUTPUT:
**      b = kx1 coefficient vector
**      cl = kx2 confidence limits
**      logl = log-likelihood
**
** GLOBALS:
**  _cmlc_Fix = variable name or column number of var with param
**             constrained to 1.
**  __title = a string with a message to be printed at the top of the output.
**
*/

proc(3) = cmlcountclprt(b,cl,logl);
    local x,mask,fmt,hdr,se,t,rse;
    local dataset,dep,vars;
    print;
    dataset = _cmlc_fn;
    dep = _cmlc_dp;
    vars = _cmlc_vr;

    if type(_cmlc_fn) == 13;
        call header("CMLCOUNT",_cmlc_fn,_cmlc_ver);
    else;
        call header("CMLCOUNT","",_cmlc_ver);
    endif;

    if rows(cl) /= rows(b);
        if not trapchk(4);
            errorlog "first and second arguments not conformable";
            retp(b,cl,logl);
        endif;
    endif;
    mask = zeros(rows(b),1)~ones(rows(b),3);
    fmt = ("*.*s"~ 10~ 8)| ("*.*lf"~ 12~ _cmlc_Precision)| ("*.*lf"~ 12~
        _cmlc_Precision)| ("*.*lf"~ 12~ _cmlc_Precision);
    if type(dataset)/=13;
        dataset = "Matrix Input";
        if rows(dep)==1;
            dep = "Col."$+ftos(dep,"*.*lf",2,0);
        else;
            dep = "Col."$+ftos(dep[1],"*.*lf",2,0)$+" & " $+ "Col." $+
                ftos(dep[2],"*.*lf",2,0);
        endif;
        t = ones(rows(vars)-1,1);
        vars[2:rows(vars)] = t.*vars[2:rows(vars)];
    else;
        if rows(dep)==1;
            dep = "" $+ lower(dep);
        else;
            dep = "" $+ lower(dep[1]) $+ " & " $+ lower(dep[2]);
        endif;
    endif;
    vars = lower(vars);
    x = vars~b~cl;
    print;
    print "Dependent Variable: "$+dep;
    print;
    print "                           "$+ftos(1-_cml_Alpha,"%*.*lf",3,2)$+""\
        " confidence limits";
    print " Parameter    Estimate   Lower Limit   Upper Limit";
    print;
    call printfm(x,mask,fmt);
    print;
    print "log-likelihood = " ftos(logl,"-*.*lf",10,_cmlc_Precision) " ";;
    print " n = " ftos(_cml_NumObs,"-*.*lf",10,0);
    if _cmlc_Fix/=0;
        if round(_cmlc_Fix) == _cmlc_Fix and _cmlc_Fix >= 1 and _cmlc_Fix
            < 131072;
            print ("Coefficient set to 1.0:    " $+ "Col." $+
                ftos(_cmlc_Fix,"*.*lf ",3,0) $+ _cmlc_Fix);
        else;
            print ("Coefficient set to 1.0:     " $+ _cmlc_Fix);
        endif;
    endif;
    __title = "";

    ndpclex;
    retp(b,cl,logl);

endp;

proc(0) = cmlcountset;
    cmlset;
    _cmlc_Start = 0;
    _cmlc_Precision = 4;
    _cmlc_ZeroTruncate = 1;
    _cmlc_Fix = 0;
    _cmlc_Dispersion = 1;
    _cmlc_Censor = 0;
    _cmlc_Inference = "MAXLIK";
    _cmlc_c1 = 0;
    _cmlc_c2 = 0;
    _cmlc_c3 = 0;
    _cmlc_fn = "";
    _cmlc_dp = "";
    _cmlc_vr = "";
endp;