probit.src

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

/*
** probit.src - Probit Analysis
**
**
** (C) Copyright 1988-1995  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.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**> probit
**
**
**  Purpose:  To estimate the binomial probit model using a GAUSS data set.
**            The details of the model are presented in the remarks section.
**
**  Format:   { vnames,b,vc,ndtran,pctall,meanx,sdx,fit,df,tol }
**              = probit(dataset,depvar,indvars);
**
**  Input:    dataset    string, name of data file
**
**            depvar     string, name of dependent variable
**                          - or -
**                       scalar, index of dependent variable.
**
**                       The value of depvar will be truncated before analysis.
**                       Thus, 1.4 is treated as category 1.
**
**            indvars    Kx1 character vector, names of independent variables.
**                          - or -
**                       Kx1 numeric vector, indices of independent variables.
**
**                       The program adds one variable for the constant term.
**
**            Defaults are provided for the following global input
**            variables.  They can be ignored unless you need control
**            over the other options provided by this procedure.
**
**            WARNING:  If you change the defaults in a command file, the new
**            values will apply in the next program you run using PROBIT unless
**            you change them back.  This can be done by running the files
**            gauss.dec and quantal.dec that contains the declare statements
**            for the globals.  If you want to make permanent changes to the
**            globals, they should be changed in gauss.dec and quantal.dec.
**
**            __altnam    provides alternative names for the variables.
**
**                        if 0 (default), the original names of the variables
**                        are used.
**
**                        if a ((1+NIVAR)x1) character vector, the first name in
**                        this vector will be used to label the dependent
**                        variable and the remaining NIVAR names will be used to
**                        label the independent variables.
**
**            __miss      global scalar, default 1.
**
**                        if 0, there are no missing values (fastest).
**
**                        if 1, do listwise deletion, drop an observation
**                        if there are any missing values among the independent
**                        and dependent variables.
**
**            __output    global scalar, default 1.
**
**                        if 1, sends results to the output device (including
**                        the screen).
**
**                        if 0, no information is sent to output.
**
**            __range     2*1 vector.  The range of record in data
**                        set used for analysis.  The first element is the
**                        starting row index, the second element is the
**                        endding row index.
**
**                        Default is the whole dataset.
**
**            __row       global scalar, default 0.
**
**                        if 0, the number of rows to read per iteration of
**                        the read loop is calculated by the program.
**
**                        if not 0, the specified number of rows will be read.
**
**            __tol       global scalar controlling the iterations.  __tol
**                        indicates the maximum difference between estimates
**                        of the coefficients in two adjacent iterations.
**
**            _qrcatnm    2x1 character vector of names of outcome categories
**                          - or -
**                        default, scalar 0 in which case names CAT1, CAT2 are
**                        used.
**
**            _qrev       global scalar, default 0.
**
**                        If 0, the parameters are in the order of comparing the
**                        first category to the second category.
**
**                        If 1, the above order is reversed. That is, the second
**                        category compares to the first category.
**
**            _qrfit      global scalar, default 0.
**
**                        if 1, print detailed goodness of fit measures,
**                        including table of observed and predicted outcomes.
**
**                        if 0, only print chi-square, -2*log-likehood and
**                        percent correctly predicted.
**
**            _qriter     global scalar, default 0.
**
**                        if 0, do not print information on iterations.
**
**                        if 1, send detailed information on iterations to
**                        the screen but not to the output device.
**
**                        if 2, send detailed information on iterations to
**                        the output device.
**
**            _qrpred     global scalar, default 0.
**
**                        if 0, predicted values will not be written to disk.
**
**                        if not 0, predicted probabilities for each outcome
**                        category are written to file ^_qrpnam with
**                        5 variables.  The first variable is the observed
**                        dependent variable.  The second variable, PBT_CDF,
**                        is the probability of the outcome occuring.  If X is
**                        the matrix of independent variables and B is the
**                        vector of parameters, this is defined as the c.d.f.
**                        of the normal distribution for X*B.  The third
**                        variable, PBT_PDF, is the p.d.f. of the normal
**                        distribution for X*B.  The fourth variable, PBT_XT,
**                        is the value of X*B.  The fifth variable, PBT_HAZ,
**                        is the ratio of PBT_PDF and PBT_CDF.  This is useful
**                        for two stage estimators of models with
**                        censored dependent variables.
**
**            _qrpredn    string name of dataset for predicted values.  The
**                        default name is "_qrpred".
**
**            _qrstat     global scalar, default 0.
**
**                        if 0, do not print descriptive statistics.
**
**                        if 1, print descriptive statistics.
**
**            _qrnewt     global scalar controlling method of optimization.
**
**                        if 1 (default) use Newton-Raphson
**                        if 0, use method of scoring
**
**           _qrmiter     maximum number of iterations, default = 1000.
**
**
**
**  Output:   vnames    a (K+2)x1 character vector containing the names of
**                      the variables in the model.  The order is:
**                      depvar|"CONSTANT"|indvars.
**
**            b       an NPARM=*(K+1) vector of parameter estimates in
**                    the order:  intercept|var1|var2|...varK.
**
**                    If errors are encountered a message will be sent to the
**                    error log.  Also, b will contain a scalar error code.  This
**                    code appears as missing unless it is translated with
**                    the command scalerr(b).  The codes are defined as:
**
**                          1  data file not found
**                          2  found undefined variables
**                          30  system singular
**                          31  too few nonmissing observations.
**                          71  number of categories of dependent variable not
**                              equal to 2
**                          72  one of the outcome categories has no cases
**                          73  an independent variable has no variation
**                          74  can't open file for predicted values
**                          75  out of disk space
**                          77  all cases were deleted
**
**            vc      NPARMxNPARM variance covariance matrix for the parameters
**                    in b.
**
**            ndtran  2x1 vector of observations.  Element 1 contains
**                    number of cases read from dataset; element 2
**                    contains number of cases left after deletion of
**                    missing cases controlled by __miss, it is the
**                    number of cases used in the analysis.
**
**            pct     the percent of cases in each of the outcome categories.
**                    Arranged in order lowest to highest.
**
**            meanx   the means based on nused cases of the independent
**                    variables in the order in indvars.
**
**            sdx     the standard deviations based on nused cases of the
**                    independent variables in the order in indvars.
**
**            fit     4x1 vector of goodness of fit measures.  Element 1 is
**                    the likelihood ratio chi-square assessing the overall
**                    fit of the model; element 2 is -2 times the log
**                    likelihood function evaluated at the estimated values;
**                    element 3 is -2 times the log likelihood function
**                    evaluated with the slopes fixed to zero; element 4 is
**                    the percentage of correct predictions from the model.
**
**            df      the degrees of freedom associated with lrx2.
**
**            tol     the tolerance reached.  If convergence was obtained,
**                    tol must be less than __tol.
**
**  Remarks:    See the manual for details on the model.
**
**  Library:    quantal
**
**  See Also:   logit, ordered, dtran
*/

#include gauss.ext
#include quantal.ext

proc (10) = probit(dataset,depvar,indvars);

    local fin,errmsg,nivar,sttime,obspred,nr, ii, dta,ylbl,y,x,minx,maxx,
        xx,xy,pct1,pct0,nused,meanx, ttl,i1,i2, pct,f,sdx,mask,omat,fmt,vv,
        bstrt,vc,dl,bchng,iter,tm,tzu, fpred,w, q,ppm,v,r,z,vcml,bnew,sdml,
        tml,pml,pbtmil,pbtxb,pbtcdf, pbtpdf, nout,tol,ll,prednm,llrest,
        lrx2,df,lrxp,mcfr2, madr2,nd,count, cur2,success,err,llfull,
        oldtrap,xlbl,nm,vtype,varindx,readdisk, ndtran,ynum,wx,delta,iswt,
        windx,wlbl,wt,wnused,cv, lastobs,counter,start,range,ky;
    delta = 6.045055872363794;

    dataset = "" $+ dataset;
    open fin = ^dataset;
    if fin == -1;
        errmsg = "ERROR:  Can't open file " $+ dataset $+ ".";
        goto errout(error(1));
    endif;

    sttime = hsec;
    obspred = zeros(4,1);           /* table of obs and predicted outcomes  */

    { nm, varindx } = indices(dataset, depvar|indvars);
    if scalerr(nm);
        goto errout(nm);
    endif;
    ylbl = nm[1];
    xlbl = "CONSTANT"|trimr(nm,1,0);
    nivar = rows(nm);
    if __weight $== 0;
        wlbl = "";
        windx = 0;
        iswt = 0;
    else;           /* INCLUDES __WEIGHT IF DEFINED */
        { wlbl, windx } = indices(dataset, __weight);
        if scalerr(wlbl);
            goto errout(wlbl);
        endif;
        iswt = 1;
        nm = nm|wlbl;
        varindx = varindx|windx;
    endif;

    { nm,vtype } = nametype(nm,__vtype);
    vtype = vtype[1];
    minx = (1E+50)*ones(nivar,1);
    maxx = -minx;
    _qrycat = miss(0,0);
    clear readdisk,nused;

    { nr,start,counter,lastobs } = _rngchk(dataset,__range);
    range = lastobs - start + 1;
    ndtran = zeros(2,1);
    ndtran[1] = range;
    nd = 0;

    call seekr(fin,start);
    count = counter;
    do while count < lastobs;
        readdisk = readdisk+1;
        dta = readr(fin,nr);
        count = count + rows(dta);