probit.src

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

        if count > lastobs;
            dta = trimr(dta,0,count-lastobs);
        endif;
        if __miss == 1;
            dta = packr(dta[.,varindx]);
        else;
            dta = dta[.,varindx];
        endif;
        nd = rows(dta);
        if nd == 0;
            continue;
        endif;
        ndtran[2] = ndtran[2] + nd;
        ynum = dta[.,1];
        if vtype;
            ynum = trunc(ynum);
        endif;
        if scalmiss(_qrycat);
            _qrycat = unique(ynum,vtype);
        else;
            _qrycat = union(ynum,_qrycat,vtype);
        endif;
    endo;
    readdisk = readdisk>1;
    if vtype == 0;
        _qrtmp = _qrycat;
    elseif rows(_qrcatnm) .ne 2;
        _qrtmp = 0 $+ "Cat_"$+ftocv(_qrycat,1,0);
    else;
        _qrtmp = _qrcatnm;
    endif;
    nused = ndtran[2];
    if nused == 0;
        errmsg = "ERROR:  No observations left after deleting missing values.";
        goto errout(error(77));
    endif;
    if rows(_qrycat) /= 2;
        errmsg = "ERROR:  Range of dependent categories is not two.";
        goto errout(error(71));
    endif;

    /* OLS Approximation and Descriptive Statistics */
    clear xx,xy,pct1,meanx,wnused,nd;
    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;
        nd = rows(dta);
        if nd == 0;
            continue;
        endif;
        x = ones(nd,1)~trimr(dta',1,0)';
        if iswt;
            wt = dta[.,cols(dta)];
            wnused = wnused+sumc(wt);
            x = trimr(x',0,1)';
        else;
            wt = ones(nd,1);
        endif;
        wx = x.*wt;
        if vtype;
            y = trunc(dta[.,1]) .== _qrycat[2];
        else;
            y = dta[.,1] .$== _qrycat[2];
        endif;
        xx = xx + wx'*x;
        xy = xy + wx'*y;
        pct1 = pct1+sumc(wt.*y);
        meanx = meanx+sumc(wx);
        minx = minc((minx~minc(x))');
        maxx = maxc((maxx~maxc(x))');
    endo;
    clear dta;
    if not iswt;
        wnused = nused;
    endif;
    meanx = meanx/wnused;
    sdx = sqrt((diag(xx)./(wnused))-(meanx.*meanx));
    pct1 = pct1/wnused;
    pct0 = 1-pct1;
    pct = pct0|pct1;

    if __output;
        call header("PROBIT",dataset,_qr_ver);
        print;
        print "CASES PROCESSED BY PROBIT:";
        print;
        print "  " ftos(nused,"%-*.*lf",1,0) " cases were kept out of "
            ftos(range,"%-*.*lf",1,0) " in file.";
        if iswt;
            print;
            print "WEIGHTING IS IN EFFECT:  ";
            print;
            print ("  Computations have used the weight variable:  " $+ wlbl);
        endif;
        call formatcv("*.*lf "~10~8);
        call formatnv("*.*lf "~10~4);
        if __output;
            print;
            print "DEPENDENT CATEORIES ARE DESIGNATED AS:";
            print;
            i1 = 1;
            i2 = 0;
            do until i1 > 2;
                if vtype ==1;
                    print ("  " $+ ftocv(_qrycat[i1],1,0) $+ " - " $+
                        _qrtmp[i1,1] $+ "    ");
                else;
                    print ("  " $+ ftocv(i1,1,0) $+ " - " $+ _qrtmp[i1,1]
                        $+ "    ");
                endif;
                if i2 == 3;
                    i2 = 0;
                    print;
                else;
                    i2 = i2 + 1;
                endif;
                i1 = i1 + 1;
            endo;
            print;
        endif;

        if __output and _qrstat;
            print ("\nDISTRIBUTION AMONG OUTCOME CATEGORIES FOR " $+ ylbl);
            print;
            call printfmt(" "~_qrtmp',zeros(1,rows(_qrtmp)+1));
            print;
            print "PROPORTION ";;
            call printfmt(pct',1);
            print;
            if minc(pct) < .00000001;
                errmsg = "ERROR:  An outcome category has no cases.";
                goto errout(error(72));
            endif;
            print;
            print "DESCRIPTIVE STATISTICS (N="$+ftos(nused,"*.*lf",1,0)$+"):";
/*
            if iswt;
                print ("                       (N="$+ftos(wnused,"*.*lf ",1,0)
                   $+ " weighted):");
            endif;
*/
            print;
            let ttl[1,5] = "  " "Mean" "Std Dev" "Minimum" "Maximum";
            call printfmt(ttl,0);
            print;
            mask = 0~1~1~1~1;
            omat = trimr(xlbl~meanx~sdx~minx~maxx,1,0);
            let fmt[5,3] = "*.*lf " 10 8 "*.*lf " 10 4 "*.*lf " 10 4 "*.*lf"\
                " " 10 4 "*.*lf " 10 4;
            call printfm(omat,mask,fmt);
        endif;
    endif;

    cv = seqa(2,1,rows(minx)-1);
    if minc(minx[cv]-maxx[cv]) == 0;
        errmsg = "ERROR:  An independent variable has no variation.";
        goto errout(error(73));
    endif;

    if wnused < nivar+1;
        errmsg = "ERROR:  Too few nonmissing observations.";
        goto errout(error(31));
    endif;

    /* Initialize variables for iterations */

    oldtrap = trapchk(1);
    trap 1,1;
    bstrt = solpd(xy,xx);           /* OLS approximation */
    trap oldtrap,1;
    if scalerr(bstrt);
        errmsg = "ERROR:  X'X matrix is singular.";
        goto errout(bstrt);
    endif;

    /* Adjust the OLS start values to approximate probit coefficients  */
    vc = (pct0).*pct1';
    dl = pct1-vc*(ln(pct1/pct0));
    vv = zeros(nivar,1);
    vv[1,.] = dl';
/*    bstrt = ((bstrt-vv)*invpd(vc));*/
    bstrt = solpd((bstrt-vv)',vc)';
    clear vv,dl;

    /* ML Iterations */

    if _qriter == 0 and __output;
        print;
        print "Iteration: ";;
    endif;
    bchng = 1000;
    iter = 1;
    clear f;
    ky = upper(chrs(key)) $== "C";
    do until abs(bchng) < __tol or iter > _qrmiter or ky;
        tm = zeros(nivar,nivar);
        tzu = zeros(nivar,1);
        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;
                wx = x.*wt;
                gosub compute;
            endo;
            clear dta;
        else;
            gosub compute;
        endif;
        goto update;

    compute:

        w = x*bstrt;
        q = pdfn(w);
        f = cdfn(w);
        ppm = f.*(1-f);
        tzu = -(wx)'((q./(ppm)).*(y-f))+tzu;

        if _qrnewt;
            v = q+w.*f;
            r = q.*(y.*(v./(f.*f))+(1-y).*((v-w)./((1-f).*(1-f))));
            tm = tm + wx'*(x.*r);
        else;
            z = q./sqrt(ppm);
            tm = tm + (wx.*z)'*(x.*z);
        endif;
        return;     /* end of compute subroutine */

    update:

        oldtrap = trapchk(1);
        trap 1,1;
        bchng = solpd(tzu,tm);
        trap oldtrap,1;
        if scalerr(tm);
            errmsg = "\nERROR:  Iterations failed.";
            goto errout(error(35));
        endif;

/*        vcml = invpd(tm);*/
/*        bnew = bstrt - vcml*tzu;*/
/*        bchng = bnew-bstrt;*/
        bnew = bstrt - bchng;
        if __output;
            if _qriter == 0;
                print " " ftos(iter,"*.*lf",1,0);;
            elseif _qriter == 1;
                output off;
                call probit_iter(iter,bstrt,bnew,xlbl);
                output on;
            elseif _qriter == 2;
                call probit_iter(iter,bstrt,bnew,xlbl);
            endif;
        endif;
        bstrt = bnew;
        iter = iter + 1;
        if readdisk;
            call seekr(fin,start);