probit.src

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

        endif;
        ky = upper(chrs(key)) $== "C";
    endo;           /* end of ML iterations loop */

    if _qrev;
        bstrt = -bstrt;
    endif;
    vcml = solpd(eye(rows(tm)),tm);
    sdml = sqrt(diag(vcml));
    tml = bstrt./sdml;
    pml = 2*cdfnc(abs(tml));
    tol = maxc(bchng);      /* tolerance at convergence */
    if _qriter == 0;
        print;
    endif;
    if __output;
        print;
        print ("ESTIMATES FROM PROBIT ANALYSIS OF VARIABLE:  " $+ ylbl);
        print;
        if iter > _qrmiter;
            print "   WARNING:  Maximum iterations - " ftos(iter-1,"*.*lf",
                1,0) " - exceeded.";
        elseif ky;
            print "   WARNING:  Convergence forced after " ftos(iter-1,"*.*"\
                "lf",1,0);;
            if _qrnewt;
                print " Newton-Raphson iterations.";
            else;
                print " Method of Scoring iterations.";
            endif;
        else;
            print "   Convergence after " ftos(iter-1,"*.*lf",1,0);;
            if _qrnewt;
                print " Newton-Raphson iterations.";
            else;
                print " Method of Scoring iterations.";
            endif;
        endif;
        print ("   Tolerance of " $+ ftos(tol,"*.*lf ",8,7) $+ "achieved af"\
            "ter" $+ ftos((hsec-sttime)/6000,"*.*lf ",5,2) $+ "minutes.");
        print;
        print "                         Probit       Std.";
        print "  Variable  Comparison  Estimate     Error      t-value     "\
            "p>|t|";
        print "  ----------------------------------------------------------"\
            "-----";
        mask = 0~0~1~1~1~1;
        let fmt[6,3] = "*.*lf " 10 8 "*.*s" 8 8 "*.*lf " 14 5 "*.*lf " 10
            4 "*.*lf " 10 2 "*.*lf " 10 3;
        if vtype;
            if _qrev;
                omat = xlbl~reshape(ftocv(_qrycat[2],1,0) $+ "/" $+
                    ftocv(_qrycat[1],1,0),nivar,1);
            else;
                omat = xlbl~reshape(ftocv(_qrycat[1],1,0) $+ "/" $+
                    ftocv(_qrycat[2],1,0),nivar,1);
            endif;
        else;
            if _qrev;
                omat = xlbl~reshape(_qrtmp[2] $+ "/" $+ _qrtmp[1],nivar,1);
            else;
                omat = xlbl~reshape(_qrtmp[1] $+ "/" $+ _qrtmp[2],nivar,1);
            endif;
        endif;
        omat = omat~bstrt~sdml~tml~pml;
        call printfm(omat,mask,fmt);
    elseif iter > _qrmiter;
        errorlog "WARNING:  Maximum iterations exceeded";
    elseif ky;
        errorlog "WARNING:  convergence forced";
    endif;

    /* Compute likelihood ratio statistics */

    ll = 0;
    if _qrpred;

    /* Open file for predicted probabilities */

        prednm = ylbl|"PBTCDF"|"PBTPDF"|"PBTXB"|"PBTHAZ";
        create fpred = ^_qrpredn with ^prednm,rows(prednm),8;
        if fpred == -1;
            errmsg = "ERROR:  Can't open file for predicted values: " $+
                _qrpredn $+ ".";
            goto errout(error(74));
        endif;
    endif;

    if readdisk;
        call seekr(fin,start);
        count = counter;
        do while count < lastobs;
            dta = readr(fin,nr);
            count = count + rows(dta);
            if count > lastobs;
                dta = trimr(dta,0,count-lastobs);
            endif;
            if __miss == 1;
                dta = packr(dta[.,varindx]);
            else;
                dta = dta[.,varindx];
            endif;
            if scalmiss(dta);
                continue;
            endif;
            if vtype;
                y = trunc(dta[.,1]) .== _qrycat[2];
            else;
                y = dta[.,1] .$== _qrycat[2];
            endif;
            x = ones(rows(dta),1)~trimr(dta',1,0)';
            if iswt;
                wt = dta[.,cols(dta)];
                x = trimr(x',0,1)';
            else;
                wt = ones(rows(dta),1);
            endif;
            clear dta;
            gosub pred;
            if _qrpred;     /* write observed outcome and probabilities  */
                if writer(fpred,y~pbtcdf~pbtpdf~pbtxb~pbtmil) /= rows(y);
                    errmsg = "ERROR:  Disk full.  Predicted file not comple"\
                        "ted.";
                    goto errout(error(75));
                endif;
            endif;
        endo;
    else;
        gosub pred;
        if _qrpred;         /* write observed outcome and probabilities  */
            if writer(fpred,y~pbtcdf~pbtpdf~pbtxb~pbtmil) /= rows(y);
                errmsg = "ERROR:  disk full, predicted file not completed.";
                goto errout(error(75));
            endif;
        endif;
    endif;

    /* log-likelihood for all parameters */
    llfull = ll;
    /* log-likelihood with only intercepts */
    llrest = sumc((wnused*pct).*ln(pct));
    /* likelihood ratio test statistics - Maddalla, p.40 */
    lrx2 = 2*(llfull - llrest);
    df = (nivar-1);
    lrxp = cdfchic(lrx2,df);
    /* McFadden's pseudo-R2 - Madalla, p.40 */
    mcfr2 = 1-(llfull/llrest);
    ii = 2/wnused;
    /* Madalla's R2 - Madalla, p. 39 */
    madr2 = 1-exp(-lrx2/wnused);
    /* Cragg & Uhler's pseudo-R2 - Madalla, p.40 */
    ii = 2/wnused;
    cur2 = (exp(ii*llfull)-exp(ii*llrest))/(1 - exp(ii*llrest));
    obspred = rotater(reshape(obspred,2,2),1)./wnused;
    /* percent correct predictions - Maddall, p.70 */
    success = 100*(obspred[1,1]+obspred[2,2]);
    obspred = wnused*(obspred~sumc(obspred'));

    if __output;
        print;
        print "MEASURES OF FIT:";
        print;
        print ("  Likelihood Ratio Chi-square:            " $+ ftos(lrx2,"*"\
            ".*lf ",10,4));
        print ("       with " $+ ftos(df,"*.*lf ",1,0) $+ " d.f., prob = "
            $+ ftos(lrxp,"*.*lf ",4,3));
        print ("  -2 Log Likelihood for full model:       " $+
            ftos(-2*llfull,"*.*lf ",10,4));
        print ("  -2 Log Likelihood for restricted model: " $+
            ftos(-2*llrest,"*.*lf ",10,4));
        print ("  Percent Correctly Predicted:            " $+
            ftos(success,"*.*lf ",10,4));
        if _qrfit;
            print ("  Madalla's pseudo R-square:              " $+
                ftos(madr2,"*.*lf ",10,4));
            print ("  McFadden's pseudo R-square:             " $+
                ftos(mcfr2,"*.*lf ",10,4));
            print ("  Cragg and Uhler's pseudo R-square:      " $+
                ftos(cur2,"*.*lf ",10,4));
            print;
            print "TABLE OF OBSERVED AND PREDICTED OUTCOMES:";
            print;
            print "           |      Predicted";
            call formatcv("*.*lf "~8~8);
            call formatnv("*.*lf "~8~0);
            print "  Observed |";;
            call printfmt(_qrtmp'~"Total",0);
            print;
            print "  ------------------------------------";
            call formatcv("*.*lf "~10~8);
            call printfmt(_qrtmp[1],0);
            print "|";;
            call printfmt(obspred[1,.],1);
            print;
            call printfmt(_qrtmp[2],0);
            print "|";;
            call printfmt(obspred[2,.],1);
            print;
            print "  ------------------------------------";
            print "     Total |";;
            call printfmt(sumc(obspred)',1);
            print;
        endif;
        if _qrpred;
            nout = rowsf(fpred);
            fpred = close(fpred);
            if not fpred and __output;
                print;
                print "PREDICTED VALUES SUCCESSFULLY WRITTEN TO DISK:";
                print;
                print ("  The file " $+ _qrpredn $+ " was written with " $+
                    ftos(nout,"*.*lf ",1,0) $+ "cases.");
                print;
                print "  The following variables are in the file:";
                print;
                print ("    Dependent variable (Y):       " $+ ylbl);
                print "    P(Y=1) = c.d.f normal of XB:  PBTCDF";
                print "    p.d.f. normal of XB:          PBTPDF";
                print "    X*B:                          PBTXB";
                print "    pdfn(xb)./cdfn(xb):           PBTHAZ";
            endif;
        endif;
    endif;

    goto done;

pred:

    /* Compute predicted probabilities */
    pbtxb = x*bstrt;
    pbtcdf = cdfn(pbtxb);           /* prob y=1 */
    pbtpdf = pdfn(pbtxb);
    pbtmil = (pbtxb >= -delta).*(pbtpdf./pbtcdf) - (pbtxb < -delta)*pbtxb;
    ll = ll+sumc(wt.*((1-y).*ln(1-pbtcdf)+y.*ln(pbtcdf)));
    obspred = obspred+ countwts((10*(y+1))+(maxindc((pbtcdf~(1-pbtcdf))')),
        11|12|21|22,wt);
    return;

done:

    if fin > 0;
        fin = close(fin);
    endif;
    ndpclex;
    lrx2 = lrx2|-2*llfull|-2*llrest|success;
    print;
    retp(ylbl|xlbl,bstrt,vcml,ndtran,pct,meanx,sdx,lrx2,df,tol);

errout:

    pop err;
    cls;
    if fin > 0;
        fin = close(fin);
    endif;
    ndpclex;
    if not trapchk(1);
        errorlog errmsg;
        print;
        end;
    endif;
    retp(0,err,0,0,0,0,0,0,0,99999);

endp;       /* probit */

proc(0) = probit_iter(iter,st,en,lbl);
    local omat,fmt,mask,st0,i;
    print;
    print "ITERATION " ftos(iter,"*.*lf",1,0) ":";
    print;
    print "  Parameter    Start Value        End Value    % Change";
    print "  -----------------------------------------------------";
    mask = 0~1~1~1;
    let fmt[4,3] = "*.*lf " 10 8 "*.*lf " 16 8 "*.*lf " 16 8 "*.*lf " 11 5;
    i = 0;
    st0 = st;
    do until i == rows(st);
        i = i+1;
        if st[i]/=0;
            st0[i] = 100*(en[i]-st[i])/st[i];
        else;
            st0[i] = -999;
        endif;
    endo;
    omat = lbl~st~en~st0;
    call printfm(omat,mask,fmt);
    if rows(omat) == 1;
        print;
    endif;
endp;       /* probit_iter */