count.src

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

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

proc _cn_lng(x);
    local x2, x3, 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;

/*
** _cn_ftosm - matrix field to string
**
**  y = _cn_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 _cn_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;

/* COUNTPRT -- Print Parameters and Estimates
**
**     { b,vc,logl } = COUNTPRT(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:
**  _cn_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) = countprt(b,vc,logl);
    local x,mask,fmt,hdr,se,t,rse;
    local dataset,dep,vars;
    print;
    dataset = _cn_fn;
    dep = _cn_dp;
    vars = _cn_vr;

    if type(_cn_fn) == 13;
        call header("COUNT",_cn_fn,_cn_ver);
    else;
        call header("COUNT","",_cn_ver);
    endif;

    if upper(_cn_Inference) $== "MAXLIK" or upper(_cn_Inference) $== "PROFILE";
        rse = real(sqrt(diag(vc)));
        se = real(sqrt(diag(_max_HessCov)));
        if ismiss(se);
            errorlog "ERROR:  Covariance matrix of parameters not invertible";
            se = miss(zeros(rows(b),1),0);
        endif;
        if ismiss(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~ _cn_Precision)| ("*.*lf"~ 10~
            _cn_Precision)| ("*.*lf"~ 10~ _cn_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,_cn_Precision) " ";;
        print " n = " ftos(_max_NumObs,"-*.*lf",10,0);
        if _cn_Fix/=0;
            if round(_cn_Fix) == _cn_Fix and _cn_Fix >= 1 and _cn_Fix < 131072;
                print ("Coefficient set to 1.0:    " $+ "Col." $+
                    ftos(_cn_Fix,"*.*lf ",3,0) $+ _cn_Fix);
            else;
                print ("Coefficient set to 1.0:     " $+ _cn_Fix);
            endif;
        endif;

    elseif upper(_cn_Inference) $== "BOOT";

        se = real(sqrt(diag(vc)));
        if ismiss(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~ _cn_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 " Parameter  Estimate   Boot S.E.";
        call printfm(x,mask,fmt);
        print;
        hdr = "bootstrap log-likelihood = " $+ ftos(logl,"-*.*lf",10,
            _cn_Precision) $+ " ";
        print hdr;
        print " n = " ftos(_max_NumObs,"-*.*lf",10,0);
        if _cn_Fix /= 0;
            if round(_cn_Fix) == _cn_Fix and _cn_Fix >= 1 and _cn_Fix < 131072;
                print ("Coefficient set to 1.0:    " $+ "Col." $+
                    ftos(_cn_Fix,"*.*lf ",3,0) $+ _cn_Fix);
            else;
                print ("Coefficient set to 1.0:     " $+ _cn_Fix);
            endif;
        endif;
    endif;
    __title = "";

    ndpclex;
    retp(b,vc,logl);

endp;

/* COUNTCLPrt -- Print Parameters and Estimates
**
**     { b,cl,logl } = COUNTCLPrt(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:
**  _cn_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) = countclprt(b,cl,logl);
    local x,mask,fmt,hdr,se,t,rse;
    local dataset,dep,vars;
    print;
    dataset = _cn_fn;
    dep = _cn_dp;
    vars = _cn_vr;

    if type(_cn_fn) == 13;
        call header("COUNT",_cn_fn,_cn_ver);
    else;
        call header("COUNT","",_cn_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~ _cn_Precision)| ("*.*lf"~ 12~
        _cn_Precision)| ("*.*lf"~ 12~ _cn_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-_max_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,_cn_Precision) " ";;
    print " n = " ftos(_max_NumObs,"-*.*lf",10,0);
    if _cn_Fix/=0;
        if round(_cn_Fix) == _cn_Fix and _cn_Fix >= 1 and _cn_Fix < 131072;
            print ("Coefficient set to 1.0:    " $+ "Col." $+ ftos(_cn_Fix,
                "*.*lf ",3,0) $+ _cn_Fix);
        else;
            print ("Coefficient set to 1.0:     " $+ _cn_Fix);
        endif;
    endif;
    __title = "";

    ndpclex;
    retp(b,cl,logl);

endp;

proc(0) = countset;
    maxset;
    _cn_Start = 0;
    _cn_Precision = 4;
    _cn_ZeroTruncate = 1;
    _cn_Fix = 0;
    _cn_Dispersion = 3;
    _cn_Censor = 0;
    _cn_Inference = "MAXLIK";
    _cn_c1 = 0;
    _cn_c2 = 0;
    _cn_c3 = 0;
    _cn_fn = "";
    _cn_dp = "";
    _cn_vr = "";
endp;