cmlcount.src

GAUSS软件的CML模块

/*
** cmlcount.src - Event Count and Duration Regression
**
**
** (C) Copyright 1994-1996  Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**        CONTENTS                                           LINE
**        --------                                           ----
**        CLOLS         compute logged OLS                   137
**        CMLCountPrt   print estimates                      392
**        CMLCountClPrt print confidence limits              537
**        CMLCountSet   reset globals to defaults            626
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
** >CLOLS
**
**  Purpose:    To compute LOGGED least squares regression. Ln(y+.5)/x
**
**  Format:     b = CLOLS(dataset,depvar,indvars);
**
**  Input:     dataset      string, name of data set.
**
**                          If this is a null string, the procedure assumes
**                          that the actual data has been passed in the
**                          next two arguments.
**
**             depvar       dependent variable.
**
**                          If dataset contains the name of a dataset, this is
**                          interpreted as:
**
**                           string, name of dependent variable
**                               or
**                           scalar, index of dependent variable. If scalar 0,
**                           the last column of the data set will be used.
**
**                          If dataset is a null string or 0, this is
**                          interpreted as:
**
**                           Nx1 vector, the dependent variable
**
**             indvars      independent variables.
**
**                          If dataset contains the name of a dataset, this is
**                          interpreted as:
**
**                            Kx1 character vector, names of independent
**                            variables
**                                 or
**                            Kx1 numeric vector, indices of independent
**                            variables
**
**                            These can be any size subset of the variables in
**                            the data set, and can be in any order.  If a
**                            scalar 0 is passed, all columns of the data set
**                            will be used except for the one used for the
**                            dependent variable.
**
**                          If dataset is a null string or 0, this is
**                          interpreted as:
**
**                            Nx1 vector, the dependent variable
**
**
**  Output:         b       Dx1 vector, the least squares estimates of
**                          parameters
**
**
**  Globals:     __con      global scalar, default 1.
**
**                             1   a constant term will be added, D = K+1.
**
**                             0   no constant term will be added, D = K.
**
**                          A constant term is always used in
**                          constructing the moment matrix m.
**
**
**
**             __miss       global scalar, default 1.
**
**                           0   there are no missing values (fastest).
**
**                           1   listwise deletion, (default) drop any
**                               cases in which missings occur.
**
**                           2   pairwise deletion, this is equivalent to
**                               setting missings to 0 when calculating m.
**                               The number of cases computed is equal to
**                               the total number of cases in the data set.
**
**              __row       global scalar, the number of rows to read per
**                          iteration of the read loop. Default 0.
**
**                          If 0, the number of rows will be calculated
**                          internally.  If you get an "Insufficient memory"
**                          error while executing OLS you can supply a value
**                          for __row that works on your system.
**
**                          The answers may vary slightly due to rounding
**                          error differences when a different number of
**                          rows is read per iteration.  You can use __row
**                          to control this if you want to get exactly the
**                          same rounding effects between several runs.
**
*/

/* common external declarations */

#include cmlcount.ext
#include cml.ext

#ifDLLCALL
external proc indices2;
#else
external proc indices2,indexcat;
#endif

proc clols(dataset,depvar,indvars);

    local idat,resid,fin,tobs,depindx,indindx,nvar,nr,dta,y0,mn,nc,cy, i,
        constflg,constvlu,vardx,std,vnames,cxx,cxxi,cxy,df, sse,nobs,mobs,
        be,b,vc,ms, prcn,nvar1,cvec,old,u2,m,constant, stdb,fout,u,str,tv,
        oldtrp,u0,cov,stdest,dd;

    clear constflg,mobs,constant;
    fin = -1;
    constvlu = 1;
    if type(dataset) == 13;
        dataset = "" $+ dataset;
        idat = 0;
    /* open file using name in variable DATASET */
        open fin = ^dataset;
        if fin == -1;
#ifUNIX
            errorlog "ERROR:  Data set " $+ dataset $+ " not found";
#else
            errorlog "ERROR:  Data set " $+ upper(dataset) $+ " not found";
#endif
            retp(zeros(indvars+1,1));
        endif;
        tobs = rowsf(fin);
        { depvar,depindx,indvars,indindx } = indices2(dataset,depvar,indvars);
        vardx = indindx|depindx;
        nobs = tobs;
        nvar = rows(indindx);
        nvar1 = nvar+1;

    /* Computation of max number of rows to read at one time */
        nr = getnr(3.5,colsf(fin));
    else;
        dta = dataset[.,indvars]~dataset[.,depvar];
        idat = 1;
        tobs = rows(dta);
        nobs = tobs;
        nvar1 = cols(dta);
        nvar = nvar1-1;
        vardx = seqa(1,1,nvar1);
    endif;

    if __miss == 2;
        old = ndpcntrl(0,0);
        call ndpcntrl(1,1);
        clear mn,nc,m,i,nobs;
        constflg = ones(1,nvar1);
        do until i == tobs;
            i = i+1;
            if idat;
                y0 = dta[i,.];
            else;
                y0 = readr(fin,1);
                y0 = y0[vardx];
            endif;
            cy = (y0 .> 0 .or y0 .< 1);
            ndpclex;
            y0 = missrv(y0,0);
            y0[.,nvar1] = ln(y0[.,nvar1] + 0.5);
            m = m+y0'*y0;
            mn = mn + y0'*cy;
            nc = nc+cy'*cy;
            nobs = nobs+(cy /= 0);
        endo;
        call ndpcntrl(old,0xffff);
        if nc == 0;
            errorlog "ERROR:  No cases left.";
            retp(zeros(nvar1,1));
        endif;
        mobs = tobs-nobs;
        mn = mn./nc;
        m = m./nc;
    elseif idat;
        if __miss == 0 and ismiss(dta);
            errorlog "ERROR:  Missing data found - __miss set to 1 (listwise)";
            __miss = 1;
        endif;
        if __miss == 1;
            dta = packr(dta);
            if scalmiss(dta);
                errorlog "ERROR:  No cases left.";
                retp(zeros(nvar1,1));
            endif;
            nobs = rows(dta);
            mobs = tobs-nobs;
        endif;
        dta[.,nvar1] = ln(dta[.,nvar1] + 0.5);
        mn = meanc(dta);
        m = moment(dta,0)/nobs;
    else;
        clear mn,m,nc;
        do until eof(fin);
            y0 = readr(fin,nr);
            y0 = y0[.,vardx];
            if __miss == 1;
                y0 = packr(y0);
                nc = nc+rows(y0);
            elseif ismiss(y0);
                errorlog "ERROR:  Missing data found.  __miss set to 1 (lis"\
                    "twise deletion)";
                __miss = 1;
                y0 = packr(y0);
                nc = nc+rows(y0);
            endif;
            if not scalmiss(y0);
                y0[.,nvar1] = ln(y0[.,nvar1] + 0.5);
                m = m+moment(y0,0);
                mn = mn + sumc(y0);
            endif;
        endo;
        if __miss == 1;
            if nc == 0;
                errorlog "ERROR:  No cases left.";
                retp(zeros(nvar1,1));
            endif;
            nobs = nc;
            mobs = tobs-nobs;
        endif;
        mn = mn/nobs;
        m = m/nobs;
    endif;

    if __miss == 2;
        constflg = indexcat(dotfeq(diag(m),diag(mn)^2),1);
    else;
        constflg = indexcat(dotfeq(diag(m),mn^2),1);
    endif;

    if scalmiss(constflg);
        constflg = 0;
    elseif rows(constflg) > 1;
        errorlog "ERROR:  Too many constants.";
        retp(zeros(nvar1,1));
    endif;

    if constflg;
        cvec = packr(miss(seqa(1,1,rows(mn)),constflg));
        if __miss == 2;
            constvlu = mn[constflg,constflg];
            mn = mn[cvec,cvec];
        else;
            constvlu = mn[constflg];
            mn = mn[cvec];
        endif;
        m = m[cvec,cvec];
        nvar1 = rows(cvec);
        nvar = nvar1 - 1;
    endif;

    if __miss == 2;
        mn = diag(mn);
    endif;

    if __con == 1 and constflg;
        errorlog "ERROR:  A column of constant value has been included amon"\
            "g the in dependent variables.";
        __con = 0;
    endif;

    if __con or constflg;
        cov = m - mn*mn';
    else;
        cov = m;
    endif;

    cxy = cov[1:nvar,nvar1];
    cxx = cov[1:nvar,1:nvar];

    oldtrp = trapchk(1);
    trap 1,1;
    cxxi = invpd(cxx);
    trap oldtrp,1;
    if scalmiss(cxxi);
        errorlog "ERROR:  Covariance matrix of independent variables failed"\
            " to invert; regression coefficients set to zero";
        retp(zeros(nvar1,1));
    endif;
    b = cxxi*cxy;
    if __con or constflg;
        constant = (mn[nvar1]-mn[1:nvar]'*b)/constvlu;
        b = constant|b;
    endif;
    if fin > 0;
        fin = close(fin);
    endif;
    retp(b);