cmlclim.src

GAUSS软件的CML模块

/*
** cmlclim.src     CMLClimits - confidence limits by inversion of
**                              Wald statistic
**
** (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.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**> CMLClimits
**
**  Purpose:    Computes confidence limits by inversion of the Wald statistic.
**
**  Format:     cl = CMLClimits(b,vc)
**
**  Input:      b     Kx1 vector, parameter estimates.
**
**              vc    KxK matrix, covariance matrix of estimates.
**
**  Output:     cl    Kx2 matrix, lower (first column) and upper (second
**                    column) limits of the selected parameters.
**
**
**  Globals:   _cml_Alpha    (1-_cml_Alpha)% two-tailed limits are computed.
**                           Default = .95
**
**             _cml_NumObs   scalar, number of observations.  Must be set.
**
**             _cml_Select   Lx1 vector, selection of columns of dataset.
**                           For example, _cml_Select = { 1, 3, 4 }.
*/

#include cml.ext

proc (1) = cmlclimits(coefs,h);
    local oldt,cl,np,tv,lv,lv1,i,k,phi0,phi1,phi2,f0,f1,f2,iter,in,
        im,ff,ff0,select,outp,LLoutput,ltst,hi,sde;

    outp = 0;
    LLoutput = 0;
#ifUNIX
    if __output;
        output off;
        print;
        print "Calculating confidence limits.....";
        print;
        print "To toggle output to screen, press O";
        output on;
        outp = 1;
        LLoutput = 1;
    endif;
#ELSE
    if __output;
        cls;
        locate 2,1;
        printdos "Calculating confidence limits.....";
        locate 3,1;
        printdos "To toggle output to screen, press O";
        if __output > 2;
            outp = 1;
            LLoutput = 1;
        else;
            outp = __output;
            LLoutput = __output;
        endif;
    endif;
#ENDIF

    if _cml_NumObs == 0;
        errorlog "ERROR: _cml_NumObs not set";
        if not trapchk(4);
            end;
        else;
            retp(error(0));
        endif;
    endif;
    coefs = vec(real(coefs));
    h = real(h);
    if rows(coefs) /= rows(h);
        errorlog "parameter vector not conformable to covariance matrix";
        if not trapchk(4);
            end;
        else;
            retp(error(0));
        endif;
    endif;
    if cols(h) /= rows(h);
        errorlog "covariance matrix not square";
        if not trapchk(4);
            end;
        else;
            retp(error(0));
        endif;
    endif;
    if scalmiss(_cml_Select) or _cml_Select == 0;
        select = seqa(1,1,rows(coefs));
    else;
        if maxc(_cml_Select) > rows(coefs) or minc(_cml_Select) < 1;
            errorlog "error in selection vector";
            if not trapchk(4);
                end;
            else;
                retp(error(0));
            endif;
        else;
            select = _cml_Select;
        endif;
    endif;

    im = abs(diag(h));
    im = im .< 1e-8 .or im .== error(0);
    sde = sqrt(diag(h).*(1-im));
    if sumc(im);
        np = rows(coefs);
        in = packr(miss((1-im).*seqa(1,1,rows(im)),0));     /* good rows in h */
        im = packr(miss(im.*seqa(1,1,rows(im)),0));         /* rows =
                                                            :: = zeros in h
                                                            */
        if scalmiss(in);
            errorlog "ERROR:  no parameters left after deleting zero ro"\
                    "ws and columns";
            if not trapchk(4);
                end;
            else;
                retp(error(0));
            endif;
        endif;
        h = h[in,in];
    else;
        in = 0;     /* all rows good */
        np = 0;
    endif;

    oldt = trapchk(1);
    trap 1,1;
    hi = solpd(eye(rows(h)),h);
    trap oldt,1;
    if scalmiss(hi);
        local r0,e0,a0,a1,a2,k1;
        { r0,e0 } = qre(h);
        k1 = sumc(abs(diag(r0)) .> 1e-8);
        if k1 == rows(h);
            hi = cholsol(eye(k1),r0);
        else;
            a0 = utrisol(r0[1:k1,k1+1:rows(r0)],r0[1:k1,1:k1]);
            if not scalmiss(a0);
                a1 = cholsol(eye(k1),r0[1:k1,1:k1]);
                if not scalmiss(a1);
                    a2 = a1*a0;
                    h = (a1~a2)|(a2'~a0'*a2);
                    hi = h[e0[e0],e0[e0]];
                else;
                    hi = error(0);
                endif;
            else;
                hi = error(0);
            endif;
        endif;
    endif;

    if scalmiss(hi);
        errorlog "ERROR: covariance matrix of parameters failed to invert";
        if not trapchk(4);
            end;
        else;
            retp(error(0));
        endif;
    endif;
    if not(in == 0);
        h = hi;
        hi = zeros(np,np);
        hi[in,in] = h;
    endif;

    cl = miss(zeros(rows(coefs),2),0);

    tv = cdftci(0.5*_cml_Alpha,_cml_NumObs-rows(coefs));
    lv = tv * tv;

    lv1 = cdftci(_cml_Alpha,_cml_NumObs-rows(coefs))^2;

    i = 0;
    do until i == rows(select);
        i = i + 1;

        if not in $== 0;
            if scalmiss(indnv(i,in));
                continue;
            endif;
        endif;

/* lower bound */

        { ff, phi0 } = _cl_set(select[i],coefs,tv*sde[i],-1);
        if scalmiss(phi0);
                errorlog "ERROR: starting value calculation failed";
            if not trapchk(4);
                end;
            else;
                goto A1;
            endif;
        endif;

        f0 = _cl_quad(hi,coefs,phi0,select[i],1,outp);
        if scalmiss(f0);
            errorlog "lower confidence limit for parameter "$+ftos(i,"%"\
                    "*.*lf",1,0)$+" failed";
            if not trapchk(4);
                end;
            else;
                goto A1;
            endif;
        endif;
        ltst = _cl_test(select[i],coefs,tv*sde[i],phi0,-1,lv,lv1);
        if abs(f0-ltst) < __tol;
            cl[select[i],1] = phi0;
            goto A1;
        endif;
        if coefs[select[i]] /= 0;
            k = .01 * 10^trunc(log(abs(coefs[select[i]])));
        else;
            k = .01;
        endif;
        phi1 = phi0 + maxc(k|k*abs(phi0));

        f1 = _cl_quad(hi,coefs,phi1,select[i],1,outp);
        if scalmiss(f1);
            errorlog "lower confidence limit for parameter "$+ftos(i,"%"\
                    "*.*lf",1,0)$+" failed";
            if not trapchk(4);
                end;
            else;
                goto A1;
            endif;
        endif;
        ltst = _cl_test(select[i],coefs,tv*sde[i],phi0,-1,lv,lv1);
        if abs(f0-ltst) < __tol;
            cl[select[i],1] = phi1;
            goto A1;
        endif;

        ff0 = 0;
        iter = 0;
        do until iter > _cml_MaxIters;
            iter = iter + 1;

            if f1 == f0;
                cl[select[i],1] = phi0;
                break;
            endif;

            phi2 = (phi1 * (lv - f0) - phi0 * (lv - f1)) / (f1 - f0);
            do while phi2 > coefs[select[i]];
                phi2 = 0.5 * (minc(phi1|phi0) + phi2);
            endo;

            ltst = _cl_test(select[i],coefs,tv*sde[i],phi2,-1,lv,lv1);
            if ltst /= lv;
                phi2 = (phi1 * (ltst - f0) - phi0 * (ltst - f1)) / (f1 - f0);
            endif;

            { ff, phi2 } = _cl_set(select[i],coefs,coefs[select[i]]-phi2,-1);
            if ff0 and ff;
                cl[select[i],1] = phi2;
                break;
            else;
                ff0 = 1;
                ff = 0;
            endif;

            f2 = _cl_quad(hi,coefs,phi2,select[i],1,outp);
            if scalmiss(f2);
                errorlog "lower confidence limit for parameter " $+
                        ftos(i,"%*.*lf",1,0)$+" failed";
                if not trapchk(4);
                    end;
                else;
                    break;
                endif;
            endif;
            ltst = _cl_test(select[i],coefs,tv*sde[i],phi2,-1,lv,lv1);
            if abs(f2-ltst) < __tol;
                cl[select[i],1] = phi2;
                break;
            endif;
            f0 = f1;
            phi0 = phi1;
            f1 = f2;
            phi1 = phi2;
        endo;

    A1:

/* upper bound */

        { ff, phi0 } = _cl_set(select[i],coefs,tv*sde[i],1);
        if scalmiss(phi0);
            errorlog "ERROR: starting value calculation failed";
            if not trapchk(4);
                end;
            else;
                continue;
            endif;
        endif;

        f0 = _cl_quad(hi,coefs,phi0,select[i],-1,outp);
        if scalmiss(f0);
            errorlog "upper confidence limit for parameter "$+ftos(i,"%"\
                    "*.*lf",1,0)$+" failed";
            if not trapchk(4);
                end;
            else;
                continue;
            endif;
        endif;
        ltst = _cl_test(select[i],coefs,tv*sde[i],phi0,1,lv,lv1);
        if abs(f0-ltst) < __tol;
            cl[select[i],2] = phi0;
            continue;
        endif;

        if coefs[select[i]] /= 0;
            k = .01 * 10^trunc(log(abs(coefs[select[i]])));
        else;
            k = .01;
        endif;
        phi1 = phi0 - maxc(k|k*abs(phi0));
        f1 = _cl_quad(hi,coefs,phi1,select[i],-1,outp);
        if scalmiss(f1);
            errorlog "upper confidence limit for parameter "$+ftos(i,"%"\
                    "*.*lf",1,0)$+" failed";
            if not trapchk(4);
                end;
            else;
                continue;
            endif;
        endif;
        ltst = _cl_test(select[i],coefs,tv*sde[i],phi1,1,lv,lv1);
        if abs(f1-ltst) < __tol;
            cl[select[i],2] = phi1;
            continue;
        endif;

        ff0 = 0;
        iter = 1;
        do until iter > _cml_MaxIters;
            iter = iter + 1;

            if f0 == f1;
                cl[select[i],2] = phi0;
                break;
            endif;
            phi2 = (phi1 * (lv - f0) - phi0 * (lv - f1)) / (f1 - f0);
            do while phi2 < coefs[select[i]];
                phi2 = 0.5 * (minc(phi1|phi0) + phi2);
            endo;

            ltst = _cl_test(select[i],coefs,tv*sde[i],phi2,1,lv,lv1);
            if ltst /= lv;
                phi2 = (phi1 * (ltst - f0) - phi0 * (ltst - f1)) / (f1 - f0);
            endif;

            { ff, phi2 } = _cl_set(select[i],coefs,phi2-coefs[select[i]],1);
            if ff0 and ff;
                cl[select[i],2] = phi2;
                break;
            else;
                ff0 = 1;
                ff = 0;
            endif;
            f2 = _cl_quad(hi,coefs,phi2,select[i],-1,outp);
            if scalmiss(f2);
                errorlog "upper confidence limit for parameter " $+
                        ftos(i,"%*.*lf",1,0)$+" failed";
                if not trapchk(4);
                    end;
                else;
                    break;
                endif;
            endif;
            ltst = _cl_test(select[i],coefs,tv*sde[i],phi2,1,lv,lv1);
            if abs(f2-ltst) < __tol;
                cl[select[i],2] = phi2;
                break;
            endif;
            f0 = f1;
            phi0 = phi1;
            f1 = f2;
            phi1 = phi2;
        endo;
    endo;

    retp(cl);
endp;


proc _cl_test(n,x,d,phi,a,lv,lv1);
    local ineqproc;

if _cml_InfCorr;

    x[n] = phi + a * d;
    if not scalmiss(_cml_C);
        if not ((_cml_C * x - _cml_D) >= 0);
            retp(lv1);
        endif;
    endif;

    if not scalmiss(_cml_Bounds);
        if rows(_cml_Bounds) == 1;
            if not(x[n] >= _cml_Bounds[1,1]);
               retp(lv1);
            endif;
            if not (x[n] <= _cml_Bounds[1,2]);
                retp(lv1);
            endif;
        else;
            if not(x[n] >= _cml_Bounds[n,1]);
                retp(lv1);
            endif;
            if not(x[n] <= _cml_Bounds[n,2]);
                retp(lv1);
            endif;
        endif;
    endif;

    if not scalmiss(_cml_IneqProc);
        IneqProc = _cml_IneqProc;
        local ineqproc:proc;
        if not(ineqproc(x) >= 0);
            retp(lv1);
        endif;
    endif;

endif;

    retp(lv);
endp;


proc(2) = _cl_set(n,x,d,a);
    local t,ineqproc,ineqjacob,g,ff,i,j,iter,df;
    clear ff;
    x[n] = x[n] + a * d;
    if not scalmiss(_cml_C);
        t = _cml_C * x - _cml_D;
        if not (t >= 0);
            i = minindc(t);
            j = packr(miss(seqa(1,1,rows(x)),n));
            x[n] = (_cml_D[i] - _cml_C[i,j]*x[j]) / _cml_C[i,n];
            t = _cml_C * x - _cml_D;
            ff = 1;
        endif;
    endif;
    if not scalmiss(_cml_Bounds);
        if rows(_cml_Bounds) == 1;
            if x[n] < _cml_Bounds[1,1];
                x[n] = _cml_Bounds[1,1];
                ff = 1;
            endif;
            if x[n] > _cml_Bounds[1,2];
                x[n] = _cml_Bounds[1,2];
                ff = 1;
            endif;
        else;
            if x[n] < _cml_Bounds[n,1];
                x[n] = _cml_Bounds[n,1];
                ff = 1;
            endif;
            if x[n] > _cml_Bounds[n,2];
                x[n] = _cml_Bounds[n,2];
                ff = 1;
            endif;
        endif;
    endif;
    if not scalmiss(_cml_IneqProc);
        IneqProc = _cml_IneqProc;