psnreg.src

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

    endif;

    if iswt;
        cv = seqa(1,1,rows(minx)-1);
    else;
        cv = seqa(1,1,rows(minx));
    endif;
    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;
    olsst = solpd(xy,xx);           /* OLS approximation */
    trap oldtrap,1;
    if scalerr(olsst);
        errmsg = "ERROR:  X'X matrix is singular.";
        goto errout(olsst);
    endif;
    if readdisk;    /* only read if data won't fit in memory */
        call seekr(fin,1);
    endif;

    /* ML Iterations for intercept only */
    bstrt = meanlny;
    if _qriter == 0 and __output;
        print;
        print "Restricted Model Iterations: ";;
    endif;
    bchng = 1000;
    iter = 1;
    clear f;
    ky = upper(chrs(key)) $== "C";
    do until abs(bchng) < __tol or iter > _qrmiter or ky;
        tm = 0;
        tzu = 0;
        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;
                y = trunc(dta[.,1]);
                xi = ones(rows(dta),1);
                if iswt;
                    wt = dta[.,cols(dta)];
                else;
                    wt = ones(rows(dta),1);
                endif;
                wx = xi.*wt;
                gosub comput1;
            endo;
            clear dta;
        else;
            xi = ones(rows(y),1);
            gosub comput1;
        endif;
        goto updat1;

    comput1:

        nparm = 1;
        w = xi*bstrt;
        lam = exp(w);
        tzu = tzu+ (sumc((y-lam).*(wt.*xi)));
        tm = tm+(sumc(wt.*xi.*xi.*lam));
        return;     /* end of compute subroutine */

    updat1:

/*        vcml = -invpd(tm);*/
        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;
        bnew = bstrt + bchng;

/*        bnew = bstrt - vcml*tzu;*/
/*        bchng = bnew-bstrt;       */

        if __output;
            if _qriter == 0;
                print " " ftos(iter,"*.*lf",1,0);;
            elseif _qriter == 1;
                output off;
                call psnreg_iter(iter,bstrt,bnew,"CONSTANT");
                output on;
            elseif _qriter == 2;
                call psnreg_iter(iter,bstrt,bnew,"CONSTANT");
            endif;
        endif;
        bstrt = bnew;
        iter = iter + 1;
    endo;           /* end of ML iterations loop */
    bonly = bstrt;

    /* ML Iterations with indep variables */
    bstrt = olsst;
    if _qriter == 0 and __output;
        print "\n      Full Model Iterations: ";;
    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;
                y = trunc(dta[.,1]);
                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:

        nparm = rows(bstrt);
        w = x*bstrt;
        lam = exp(w);
        tzu = tzu+ (sumc((y-lam).*wx));
        tm = tm+((lam.*x)'wx);
        return;     /* end of compute subroutine */

    update:

        vcml = -invpd(tm);
        bnew = bstrt - vcml*tzu;
        bchng = bnew-bstrt;
        if __output;
            if _qriter == 0;
                print " " ftos(iter,"*.*lf",1,0);;
            elseif _qriter == 1;
                output off;
                call psnreg_iter(iter,bstrt,bnew,"CONSTANT");
                output on;
            elseif _qriter == 2;
                call psnreg_iter(iter,bstrt,bnew,"CONSTANT");
            endif;
        endif;
        bstrt = bnew;
        iter = iter + 1;
        ky = upper(chrs(key)) $== "C";
    endo;           /* end of ML iterations loop */
    vcml = -vcml;
    sdml = sqrt(diag(vcml));
    tml = bstrt./sdml;
    pml = 2*cdfnc(abs(tml));
    tol = maxc(bchng);      /* tolerance at convergence */
    print;
    print;
    if __output;
        print ("ESTIMATES FROM POISSON REGRESSION ON:  " $+ ylbl);
        print;
        if iter > _qrmiter;
            print "   WARNING:  Maximum iterations - " ftos(iter-1,"*.*lf",
                1,0) " - exceeded";
        elseif ky;
            print "   WARNING:  Convergence forced at " ftos(iter-1,"*.*lf"
                ,1,0);;
            print " Newton-Raphson iterations.";
        else;
            print "   Convergence after " ftos(iter-1,"*.*lf",1,0);;
            print " Newton-Raphson iterations.";
        endif;
        print "   Tolerance of " ftos(tol,"*.*lf",8,7) " achieved.";
        print;
        print "               Poisson       Std.";
        print "  Variable     Estimate     Error      t-value     p>|t|";
        print "  -----------------------------------------------------";
        mask = 0~1~1~1~1;
        let fmt[5,3] = "*.*lf " 10 8 "*.*lf " 12 5 "*.*lf " 10 4 "*.*lf "
            10 2 "*.*lf " 10 3;
        omat = xlbl~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 */

    /* log-likelihood for all parameters */
    llfull = 0;
    /* log-likelihood with only intercepts */
    llrest = 0;

    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;
            y = trunc(dta[.,1]);
            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 getll;
        endo;
    else;
        gosub getll;
    endif;

    /* 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 = abs(1-(llfull/llrest));

    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",3,4));
        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 ("  Madalla's pseudo R-square:              " $+ ftos(mcfr2,""\
            "*.*lf",10,4));
    endif;

    goto done;

getll:

    w = x*bstrt;
    lam = exp(w);
    llfull = llfull+sumc(wt.*(-lam + y.*w - ln(y!)));
    w = ones(rows(x),1)*bonly;
    lam = exp(w);
    llrest = llrest+sumc(wt.*(-lam + y.*w - ln(y!)));

    return;

done:

    if fin > 0;
        fin = close(fin);
    endif;
    ndpclex;
    lrx2 = lrx2|-2*llfull|-2*llrest;
    retp(ylbl|xlbl,bstrt,vcml,ndtran,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,99999);

endp;       /* psnreg */

proc(0) = psnreg_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;       /* psnreg_iter */