maxdens.src

没有说明

/*
** maxdens.src    MAXDensity  -  Kernel Density Estimation
**
**
** (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.
**
**   Adapted from a procedure written by Gary King, Department of Government,
**   Harvard University
**
**   Reference:  B.W. Silverman. 1986. _Density Estimation for Statistics
**      and Data Analysis_.  London: Chapman and Hall.
**
**   An Application: Gary King. "Constituency Service and Incumbency
**      Advantage,"  _British Journal of Political Science_, 21, 1
**      (January, 1991): 119--128.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**> MAXdensity
**
**  Purpose:   To compute kernel density estimate and plot.
**
**
**  Format:    ( px,py,sm } = MAXdensity(dataset,pars);
**
**  Input:     dataset     string, name of GAUSS dataset
**                         containing data.
**
**             pars        Kx1 vector, selected columns for
**                         estimation and display.
**
**
**  Output:    px          _max_NumPointsx1 vector, abscissae.
**             py          _max_NumPointsx1 vector, ordinates.
**             sm          Kx1, or Nxk, or Nx1 smoothing coefficients.
**
**  Remarks:
**
**       kernel density plots of the selected parameters are
**       generated.
**
**  Globals:
**
**     _max_Kernel     Kx1 character vector, type of kernel:
**
**                    NORMAL - normal kernel
**                      EPAN - Epanechnikov kernel
**                     BIWGT - biweight kernel
**                    TRIANG - triangular kernel
**                   RECTANG - rectangular kernel
**                   TNORMAL - truncated normal kernel
**
**            If _max_Kernel is scalar, the kernel is the same
**            for all parameters.  Default = { NORMAL };
**
**  _max_NumPoints    scalar, number of points to be computed for plots
**
**  _max_EndPoints    Kx2 matrix, lower (in first column) and upper
**                    (in second column) endpoints of density.  Default is
**                    minimum and maximum, respectively, of the parameter
**                    values.  If 1x2 matrix, endpoints will be the same
**                    for all parameters.
**
**  _max_Smoothing    Kx1 vector or Nx1 vector or NxK matrix, smoothing
**                    coefficients for each plot.  If scalar, smoothing
**                    coefficient will be the same for each plot.  If zero,
**                    smoothing coefficient will be computed by CMLdensity.
**                    If matrix, smoothing coefficient will be different for
**                    each observation.
**                    Default = 0;
**
**  _max_Truncate     Kx2 matrix, lower (in first column) and upper (in
**                    second column) truncation limits for truncated normal
**                    kernel. If 1x2 matrix, truncations limits will be the
**                    same for all plots.  Default is minimum and maximum,
**                    respectively.
*/

#include maxlik.ext

proc 3 = MAXDensity(dataset,pars);
    local strt,endd,h,kern,pts,px,py,smth,std,Tleft,Tright,i0;
    local fhandle,k1,vnames,vindx,t1,t2,sq5,kernm,kern0,y0,
         nobs,ivar,i,px0,py0,smth0,z,oaw,vv,str,scalarsm,obsvecsm,
         varsvecsm,matrixsm;

    kernm = { NORMAL, EPAN, BIWGT, TRIANG, RECTANG, TNORMAL };
    kern0 = { 1, 2, 3, 4, 5, 1 };
    sq5 = sqrt(5);

    if type(pars) == 13;
        pars = stof(pars);
    endif;

    if type(dataset) == 13 and dataset $/= "";
        fhandle = -1;
        open fhandle = ^dataset;
        if fhandle == -1;
            if not trapchk(4);
                errorlog dataset $+ " could not be opened";
            endif;
            retp(error(0),error(0));
        endif;

        call seekr(fhandle,1);
        if pars $== 0;
            vindx = 0;
            vnames = getname(dataset);
        else;
            { vnames,vindx } = indices(dataset,pars);
        endif;

        dataset = {};
        k1 = getnr(6,colsf(fhandle));
        do until eof(fhandle);
            y0 = readr(fhandle,k1);
            dataset = dataset | y0[.,vindx];
        endo;
        clear y0;
        if fhandle > 0;
           fhandle = close(fhandle);
        endif;
    else;
        if not (pars $== 0);
            dataset = dataset[.,pars];
            vnames = "PAR_"$+pars;
        else;
            vnames = "PAR_"$+seqa(1,1,cols(dataset));
        endif;
    endif;
    nobs = rows(dataset);

    pts = int(_max_NumPoints[1,1]);
    if pts <= 2;
        if trapchk(4);
            errorlog "ERROR: _max_NumPoints must be greater than 2 - "\
                       "Try _max_NumPoints = 100";
            retp(error(0),error(0));
        endif;
    endif;

    px0 = zeros(pts,cols(dataset));
    py0 = zeros(pts,cols(dataset));
    ScalarSm = rows(_max_Smoothing) == 1 and cols(_max_Smoothing) == 1;
    ObsVecSm = rows(_max_Smoothing) == rows(dataset) and
                cols(_max_Smoothing) == 1;
    VarsVecSm = rows(_max_Smoothing) == cols(dataset) and
                cols(_max_Smoothing) == 1;
    MatrixSm = rows(_max_Smoothing) == rows(dataset) and
                      cols(_max_Smoothing) == cols(dataset);


    if ScalarSm or VarsVecSm;
        smth0 = zeros(cols(dataset),1);
    elseif ObsVecSm or MatrixSm;
        smth0 = zeros(rows(dataset),cols(_max_Smoothing));
    else;
        errorlog "MAXDENSITY: _max_Smoothing not conformable to data";
        if not trapchk(4);
           end;
        else;
           retp(error(0),error(0),error(0));
        endif;
    endif;

#IFUNIX
    if __output;
        if sysstate(26,0) == 2;
            oaw = WinGetActive;
            vv = { 100,100,640,480,40,80,1,6,15,0,0,2,2 };
        else;
            if not trapchk(4);
                errorlog "MAXDensity:  graph not produced - not in windows"\
                    " environment";
            endif;
        endif;
    endif;
#ENDIF


    ivar = 1;
    do until ivar > cols(dataset);
        if rows(_max_EndPoints) == 1;
            i0 = 1;
        else;
            i0 = ivar;
        endif;
        strt = _max_EndPoints[i0,1];
        endd = _max_EndPoints[i0,2];

        if strt > endd;
            if trapchk(4);
                 if __output;
                     errorlog "ERROR:  endpoints improperly defined";
                     retp(error(0),error(0));
                 endif;
            endif;
        elseif strt==endd;
            strt = minc(dataset[.,ivar]);
            endd = maxc(dataset[.,ivar]);
        endif;

        if ScalarSm or ObsVecSm;
            smth = _max_Smoothing;
        elseif VarsVecSm;
            smth = _max_Smoothing[ivar];
        elseif MatrixSm;
            smth = _max_Smoothing[.,ivar];
        endif;
        if (smth < 0) and (smth /= -1);
            if trapchk(4);
                errorlog "ERROR: _max_Smoothing must be -1 or > than 0";
                retp(error(0),error(0));
            endif;
        endif;
        py = sortc(dataset[.,ivar],1);
        if smth == 0;
            std = minc((py[int(3*rows(py) / 4)] -
                            py[int(rows(py)/4)]) / 1.34 | stdc(py));
            smth = 0.9 * std * (rows(py)^(-0.2));
        endif;
        if rows(_max_Kernel) == 1;
            kern = _ml_check(_max_Kernel,1,kernm,kern0,1);
        else;
            kern = _ml_check(_max_Kernel[ivar],1,kernm,kern0,1);
        endif;

        Tleft = 0;
        Tright = 0;

        if kern == 2;
            if rows(_max_Truncate) == 1;
                i0 = 1;
            else;
                i0 = ivar;
            endif;
            if _max_Truncate[i0,1] == _max_Truncate[i0,2];
                Tleft = strt;
                Tright = endd;
            else;
                Tleft = _max_Truncate[i0,1];
                Tright = _max_Truncate[i0,2];
            endif;
        endif;

        t1 = (endd - strt) / pts;
        px = seqa(strt+0.4*t1,t1,pts);
        py = px;

        i = 1;
        do while i <= pts;
            if MatrixSm;
                z = (py[i,1] - dataset[.,ivar])./smth[i];
            else;
                z = (py[i,1] - dataset[.,ivar])./smth;
            endif;
            if (kern == 1);         /* NORMAL kernel */
                if MatrixSm;
                    py[i,1] = sumc((1 / sqrt(2*pi)) *
                             exp(-0.5*(z.*z))./smth[i])/nobs;
                else;
                    py[i,1] = sumc((1 / sqrt(2*pi)) *
                             exp(-0.5*(z.*z))./smth)/nobs;
                endif;
            elseif (kern == 2);     /* EPANECHNIKOV kernel */
                t1 = (abs(z) .< sq5);
                t2 = code(t1,sq5|1);
                if MatrixSm;
                    py[i,1] = sumc( t1.*(.75*(1-.2*(z.*z))./t2)./smth[i])/nobs;
                else;
                    py[i,1] = sumc(t1.*(.75*(1-.2*(z.*z))./t2)./smth)/nobs;
                endif;
            elseif (kern == 3);     /* BIWEIGHT kernel */
                if MatrixSm;
                    py[i,1] = sumc(.9375*(abs(z).<1).*
                                    ((1-(z.*z))^2)./smth[i])/nobs;
                else;
                    py[i,1] = sumc(.9375*(abs(z).<1).*
                                    ((1-(z.*z))^2)./smth)/nobs;
                endif;
            elseif (kern == 4);     /* TRIANGULAR kernel */
                if MatrixSm;
                    py[i,1] = sumc((abs(z).<1).*(1 - abs(z))./smth[i])/nobs;
                else;
                    py[i,1] = sumc((abs(z).<1).*(1 - abs(z))./smth)/nobs;
                endif;
            elseif (kern == 5);     /* RECTANGULAR kernel */
                if MatrixSm;
                    py[i,1] = sumc(0.5*(abs(z).<1)./smth[i])/nobs;
                else;
                    py[i,1] = sumc(0.5*(abs(z).<1)./smth)/nobs;
                endif;
            elseif (kern == 6);     /* TRUNCATED NORMAL kernel */
                if py[i,1] > tleft and z < tright;
                    if MatrixSm;
                        py[i,1] = sumc(pdfn(z)./(cdfn((tright -
                        dataset[.,ivar])/smth[i])-cdfn((tleft -
                        dataset[.,ivar])/smth[i]))./smth[i])/nobs;
                    else;
                        py[i,1] = sumc(pdfn(z)./(cdfn((tright -
                        dataset[.,ivar])/smth)-cdfn((tleft -
                        dataset[.,ivar])/smth))./smth)/nobs;
                    endif;
                else;
                    py[i,1] = 0;
                endif;
            else;
                errorlog "_max_Kernel incorrectly specified";
                if not trapchk(4);
                   end;
                else;
                    retp(error(0),error(0));
                endif;
            endif;

            i = i + 1;
        endo;

        px0[.,ivar] = px;
        py0[.,ivar] = py;
        if ScalarSm or VarsVecSm;
            smth0[ivar] = smth;
        elseif ObsVecSm or MatrixSm;
            smth0[.,ivar] = smth;
        endif;
        str = ""$+vnames[ivar];


#ifunix
        if __output and sysstate(26,0) == 2;
           call WinSetActive(WinOpenPQG(vv,str,"Density"));
            xlabel(str);
            xy(px,py);
            vv[1:2] = vv[1:2] + 20;
            call WinSetActive(oaw);
        endif;
#else
        if __output;
            xlabel(str);
            xy(px,py);
        endif;
#endif
        ivar = ivar + 1;
    endo;

    retp(px0,py0,smth0);
endp;