cmlexpon.src

GAUSS软件的CML模块

/*
** cmlexpon.src - Constrained Exponential Regression Model
**                (For the analysis of duration data)
**
** (C) Copyright 1994-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.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**  FORMAT:    { b,vc,llik } = CMLExpon(dataset,dep,ind);
**
**  INPUT:
**      dataset = name of Gauss dataset on disk or matrix in memory
**      dep     = dependent variable name, or column number of input matrix;
**                elements should be durations in time so that 0<y<infinity.
**      ind     = vector of independent variable names, or column numbers
**                of input matrix
**  OUTPUT:
**      b    = vector of effect parameters that maximize the likelihood
**      vc   = variance-covariance matrix of b
**      llik = value of the log-likelihood at the maximum
**
**  GLOBALS:
**
**  _cmlc_Inference    = CML for constrained maximum likelihood estimates
**                     = BOOT for bootstrapped estimates
**                     = BAYES for Bayesian inference
**
**  _cmlc_Censor scalar,
**
**                        0 = exponential gamma model, no censoring (default),
**                   symbol = use this variable from dataset to censor.
**                            Code each row as 0 if censored or 1 if not.
**                  integer = use this column of input matrix to censor.
**                           Code each row as 0 if censored or 1 if not.
**
**   _cmlc_Start   choose method of calculating starting values.
**                     0 = LS (default),
**                     1 = set to vector stored in _cnsval,
**                     2 = rndu-0.5,
**                     3 = zeros, or set to vector
**
**      __output    1  =  print output to screen (default),
**                  0  =  do not print to screen
**
**  OTHER GLOBALS:
**      see CML.
**
**  EXAMPLE:
**      dataset = "coal";
**      dep = { months };
**      ind = { nparties, unemp, polariz };
**      _cmlc_Censor = { ciep };
**      call CMLExpon(dataset,dep,ind);
**
**  REFERENCE:  Gary King, Jim Alt, Nancy Burns, and Michael Laver. "A
**              Unified Model of Coalition Duration in Parlimentary
**              Democracies," photocopy, Department of Government, Harvard U.
*/

#include gauss.ext
#include cmlcount.ext
#include cml.ext

@ starting values @

proc _cmlc_svexp(dataset,dep,ind);
    local pars,b;
    if ind == 0;
        pars = 1;
    else;
        pars = rows(ind)+1;
    endif;
    if _cmlc_Start == 0;
        if pars == 1;
            b = 0;
        else;
            b = -clols(dataset,dep,ind);
        endif;
    elseif _cmlc_Start == 1;
        b = _cmlc_StartValues;
        if rows(b)/=pars;
            errorlog "ERROR:  b is the wrong size for _cmlc_Start.";
            end;
        endif;
    elseif _cmlc_Start == 2;
        b = rndu(pars,1)-0.5;
    elseif _cmlc_Start == 3;
        b = zeros(pars,1);
    else;
        b = _cmlc_Start;
        if rows(b)/=pars;
            errorlog "ERROR:  ROWS(_cmlc_Start) is wrong.";
            end;
        endif;
    endif;
    retp(b);
endp;


@ likelihood function @

proc _cmlc_liexp(b,dta);
    local y,x,xb,n,res,c,dd;
    y = dta[.,1];
    n = rows(dta);
    c = cols(dta)-(_cmlc_Censor/=0);
    x = ones(n,1);
    if c>1;
        x = x~dta[.,2:c];
    endif;

    if _cmlc_Censor/=0;
        dd = dta[.,cols(dta)];
    else;
        dd = 1;
    endif;

    xb = x*b;
    res = -(dd.*xb)-exp(-xb).*y;
    retp(res);
endp;

@ the main proc @
proc 3 = CMLExpon(dataset,dep,ind);
    local vars,se,b,logl,g,vc,st,ret;
    _cml_CovPar = 3;
    _cmlc_fn = dataset;
    if dep $== 0;
        errorlog "ERROR:  DEP must be a variable name or number.";
        end;
    endif;
    if type(dataset) /= 13;
        if _cmlc_Censor > cols(dataset);
            errorlog "ERROR:  If dataset is a matrix in memory, "\
             "_cmlc_Censor must be the number of a column in the "\
             "matrix.";
            end;
        endif;
    endif;
    if type(dataset) /= 13;
        if (maxc(ind) > cols(dataset)) or (dep > cols(dataset));
            if not trapchk(4);
                errorlog "\nERROR:  If dataset is a matrix, the DEP "\
                    "and IND must be column numbers.\n";
                end;
            endif;
            retp(error(0),error(0),error(0));
        endif;
    endif;
    if (type(dataset) == 13) and (type(_cmlc_Censor) == 13);
        _cmlc_Censor = indcv(_cmlc_Censor,getname(dataset));
    endif;
    st = _cmlc_svexp(dataset,dep,ind);
    if __title $==  "";
      if _cmlc_Censor/=0;
          __title = "Censored ";
      endif;
          __title = __title $+ "Exponential Regression Model";
    endif;
    if ind == 0;
        vars = dep;
    else;
        vars = dep|ind;
    endif;
    if _cmlc_Censor/=0;
        vars = vars|_cmlc_Censor;
    endif;
    _cml_GradProc = &_cmlc_der1exp;
    _cml_HessProc = &_cmlc_der2exp;
    local infm,inf0,lcInf;
    infm = { CML, BOOT, BAYES };
    inf0 = { 1, 2, 3 };
    LcInf = _cml_check(_cmlc_Inference,1,infm,inf0,1);
    if LcInf == 1;
        { b,logl,g,vc,ret } = cml(dataset,vars,&_cmlc_liexp,st);
    elseif LcInf == 2;
        { b,logl,g,vc,ret } = cmlboot(dataset,vars,&_cmlc_liexp,st);
    elseif LcInf == 3;
        { b,logl,g,vc,ret } = cmlbayes(dataset,vars,&_cmlc_liexp,st);
    endif;
    if ret /= 0;
        errorlog "ERROR: Model estimation failed.";
        end;
    endif;
    if type(dataset) == 13;
        vars = "beta0";
        if ind/=0;
            vars = vars|ind;
        endif;
    else;
        vars = "beta0";
        if ind/=0;
            vars = vars|
            ((0 $+ "Col." $+ zeros(rows(ind),1))$+_cmlc_ftosm(ind,2));
        endif;
    endif;
    _cmlc_vr = vars;
    _cmlc_dp = dep;
    ndpclex;
    retp(b,vc,logl*_cml_NumObs);
endp;



proc _cmlc_der1exp(b,dta);

@ First derivatives @

    local y,x,xb,n,res,c,dd;
    y = dta[.,1];
    n = rows(dta);
    c = cols(dta)-(_cmlc_Censor/=0);
    x = ones(n,1);
    if c>1;
        x = x~dta[.,2:c];
    endif;

    if _cmlc_Censor/=0;
        dd = dta[.,cols(dta)];
    else;
        dd = 1;
    endif;

    xb = x*b;
    res = (-dd+exp(-xb).*y).*x;
    retp(res);
endp;

proc _cmlc_der2exp(b,dta);

@ Second derivatives @

    local y,x,xb,n,res,c,dd;
    y = dta[.,1];
    n = rows(dta);
    c = cols(dta)-(_cmlc_Censor/=0);
    x = ones(n,1);
    if c>1;
        x = x~dta[.,2:c];
    endif;

    if _cmlc_Censor/=0;
        dd = dta[.,cols(dta)];
    else;
        dd = 1;
    endif;

    xb = x*b;
    res = sqrt(exp(-xb).*y);
    res = -moment(res.*x,0);
    retp(res);
endp;