maxbayes.src

没有说明

/*
**   maxbayes.src     MAXBAYES - Bayesian Inference from Weighted
**                    Likelihood Bootstrap Likelihood Bootstrap
**
**
** (C) Copyright 1995-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.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**> MaxBayes
**
**  Purpose  generates simulated posterior distribution of parameters
**
**  Format: { x,f,g,cov,retcode } = MAXBayes(dataset,vars,&fct,start)
**
**  Input:    dataset      string containing name of GAUSS data set, or
**                         name of data matrix stored in memory
**
**            vars         Kx1 vector or scalar zero.  If Kx1, character
**                         vector of labels selected for analysis,
**                         or numeric vector of column numbers in data
**                         set of variables selected for analysis.
**                         If scalar zero, all columns are selected.
**
**            fct          the name of a procedure that returns either
**                         the log-likelihood for one observation or a
**                         vector of log-likelihoods for a matrix of
**                         observations
**
**            start        Kx1 vector, start values
**
**
**  Output:   b            Kx1 vector of mean of posterior distribution
**                         of parameters
**
**            f            scalar, mean weighted bootstrap log-likelihood
**
**            g            Kx1 vector, mean gradient of weighted bootstrap
**
**            cov          KxK covariance matrix of simulated posterior
**
**            retcode      return code:
**
**                            0   normal convergence
**                            1   forced exit
**                            2   maximum number of iterations exceeded
**                            3   function calculation failed
**                            4   gradient calculation failed
**                            5   Hessian calculation failed
**                            6   step length calculation failed
**                            7   function cannot be evaluated at initial
**                                parameter values
**                            8   number of elements in the gradient vector
**                                inconsistent with number of starting values
**                            9   gradient function returned a column vector
**                                rather than the required row vector
**                           10   secant update failed
**                           11   maximum time exceeded
**                           12   weights could not be found
**                           20   Hessian failed to invert
**                           34   data set could not be opened
**                           99   termination condition unknown
**
**
**  Globals:
**
**   _max_BayesAlpha   exponent of the Dirichlet random variates used in
**                     the weights for the weighted bootstrap. See
**                     Newton and Raftery, "Approximate Bayesian Inference
**                     with the Weighted Likelihood Bootstrap",  J.R.Statist.
**                     Soc. B (1994), 56:3-48.  Default = 1.4.
**
**    _max_PriorProc   pointer to proc for computing prior.  This proc
**                     takes the parameter vector as its only argument,
**                     are returns a scalar probability.  If a proc is not
**                     provided, a uniform prior is assumed.
**
**    _max_NumSample   scalar, number of re-samples in the weighted likelihood
**                     bootstrap.
**
**    _max_BootFname   string, file name of GAUSS dataset (do not include
**                     the .DAT extension) containing simulated posterior
**                     of the parameters.  If not specified, MAXBayes
**                     will select the file name, BOOTxxxx where xxxx is
**                     0000 incremented by 1 until a name is found that
**                     doesn't exist on the current directory.
**
**   MAXLIK globals are relevant.  See MAXLIK.SRC for their description.
*/

#include maxlik.ext


proc (5) = maxbayes(dataset,var,lfct,start);
    local x,f,g,h,retcode,prior,isPrior,alpha,oldt,lg,ii,jj,tme,itdta,it0,
          nobs,tt0,fhandle,fout,mm,mn,gg,iter,title0,i,j,f0,ncase,ttime,
          np,fhat,terrell,num,ch,datx,datf,ofname,L1,wgt;

    alpha = _max_BayesAlpha;

    if not scalmiss(_max_PriorProc);
        prior = _max_PriorProc;
        local prior:proc;
        isPrior = 1;
    else;
        isPrior = 0;
    endif;


    if type(dataset) == 13;
        if dataset $/= "";
            open fhandle = ^dataset;
            if fhandle == -1;
                errorlog dataset $+ " could not be opened";
                retp(start,error(0),error(0),error(0),error(34));
            endif;
            nobs = rowsf(fhandle);
            fhandle = close(fhandle);
        else;
            errorlog "data set not specified";
            retp(start,error(0),error(0),error(0),error(34));
        endif;
    else;
        nobs = rows(dataset);
    endif;

    tme = 0;
    ttime = date;
    iter = 1;
    datx = zeros(_max_NumSample,rows(start));
    datf = zeros(_max_NumSample,1);
    tt0 = ftos(_max_NumSample,"%0*.*lf",1,0);
    title0 = __title;
    clear mn,gg,itdta,it0;
    do until iter > _max_NumSample or tme > _max_MaxTime;

        if __weight == 0;
            wgt = (-ln(rndu(nobs,1)))^alpha;
        else;
           wgt = __weight .* (-ln(rndu(nobs,1)))^alpha;
        endif;

        title0 = __title $+ " - " $+ ftos(iter,"%0*.*lf",1,0) $+
                        " of " $+ tt0 $+ " -";

    { x,f,g,L1,retcode,L1,it0,L1,L1,_max_dat,L1,_max_row,
         _max_dsn,_max_Diagnostic } = _Max(dataset,var,lfct,start,
    _max_Algorithm,
    _max_Diagnostic,
    _max_GradCheckTol,
    _max_LineSearch,
                  0,
    _max_GradMethod,
    _max_GradStep,
    _max_Delta,
    _max_Extrap,
    _max_GradProc,
    _max_GradTol,
    _max_HessProc,
    _max_Interp,
    _max_Key,
    _max_Lag,
    _max_MaxIters,
    _max_MaxTime,
    _max_MaxTry,
    _max_NumObs,
    _max_ParNames,
    _max_RandRadius,
    _max_Options,
    _max_UserSearch,
    _max_UserNumGrad,
    _max_UserNumHess,
    _max_Active,
    _max_dat,
    _max_dsn,
    _max_row,
    __altnam,
    __output,
    __row,
    title0,
     wgt
    );

        if retcode == 0;
            datx[iter,.] = x';
            datf[iter] = f;
            mn = mn + x;
            gg = gg + g;
            start = mn / iter;
            itdta = itdta + it0;
        else;
            datx[iter,.] = miss(zeros(1,cols(datx)),0);
            datf[iter] = error(0);
        endif;

        tme = ethsec(ttime,date)/6000;
        iter = iter + 1;
    endo;
    datx = packr(datx);
    datf = packr(datf);
    mm = vcx(datx);
    f0 = meanc(datf);
    gg = g / rows(datf);

    itdta[1:2] = itdta[1:2] / iter;
    itdta[3] = "BAYES";
    _max_IterData = itdta;
    if tme >= _max_MaxTime;
        retcode = 11;
    else;
        retcode = 0;
    endif;


    np = rows(start);
    terrell = ( ((np+8)^((np+6)/2)) /
                (16*rows(datx)*(np+2)*gamma((np+8)/2)) )^(2/(np+4));

    oldt = trapchk(1);
    trap 1,1;
    ch = inv(chol(mm*terrell));
    trap oldt,1;
    if scalmiss(ch);
         errorlog "estimated covariance matrix not positive definite";
         retp(start,error(0),error(0),error(0),error(20));
    endif;

#IFUNIX
#ELSE
    if __output;
        cls;
        locate 1,1;
        ii = 19;
        jj = 1;
        printdos "Computing weights ";
    endif;
#ENDIF
    fhat = zeros(rows(datx),1);
    i = 1;
    do until i > rows(datx);
#IFUNIX
#ELSE
        if __output;
            if ii > 78;
                ii = 1;
                jj = jj + 1;
            else;
                ii = ii + 1;
            endif;
            locate jj,ii;
            printdos ".";
        endif;
#ENDIF
        if isPrior;
            f0 = rows(datx)*datf[i]*prior(datx[i,.]');
        else;
            f0 = rows(datx)*datf[i];
        endif;
        if abs(f0) > 1e-16;
            f0 = sumc(exp(sumc((-.5 * ((datx[i,.] - datx)*ch)^2)'))/f0);
            if abs(f0) > 1e-16;
                fhat[i] = 1 / f0;
            else;
                fhat[i] = 0;
            endif;
        else;
            fhat[i] = 0;
        endif;

        i = i + 1;
    endo;
    fhat = fhat / sumc(fhat);

    if _max_BootFname $== "";
        ofname = tempname("","boot",".dat");
        if ofname $== "";
            errorlog "ERROR:  file name selection for bootstrap dataset"\
                               " failed";
            retp(start,error(0),error(0),error(0),error(34));
        endif;
        _max_BootFname = ofname;
    else;
        ofname = _max_BootFname;
    endif;

    create fout = ^ofname with PAR_,rows(start),8;

#IFUNIX
#ELSE
    if __output;
        cls;
        locate 1,1;
        ii = 18;
        jj = 1;
        printdos "Writing data set ";
    endif;
#ENDIF


    mm = 0;
    mn = 0;
    ncase = 0;
    i = 1;
    do until i > rows(fhat);
#IFUNIX
#ELSE
        if __output;
            if ii > 78;
                ii = 1;
                jj = jj + 1;
            else;
                ii = ii + 1;
            endif;
            locate jj,ii;
            printdos ".";
        endif;
#ENDIF
        num = _max_rndp(fhat[i]*rows(datx));
        j = 1;
        do until j > num;
            call writer(fout,datx[i,.]);
            ncase = ncase + 1;
            mn = mn + datx[i,.];
            mm = mm + moment(datx[i,.],0);
            j = j + 1;
        endo;
        i = i + 1;
    endo;
    mn = mn' / ncase;
    mm = mm / ncase - mn * mn';
    fout = close(fout);

#IFUNIX
#ELSE
    if __output;
        cls;
    endif;
#ENDIF
    retp(mn,f0,gg,mm,0);
endp;



proc _max_rndp(l);
    local x,r,r0;
    x = 0;
    r0 = rndu(1,1);
    r = exp(-l);
    do while r < R0;
       r0 = r0 - r;
       x = x + 1;
       r = r*l/x;
    endo;
    retp(x);
endp;