lnpdfn.src

没有说明

/*
**  lnpdfn.src   compute log of normal probability function
**
** (C) Copyright 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.
**
**  Format                                                        Line
** =========================================================================
**  z = lnpdfn(x);                                                 26
**  z = lnpdfmvn(x,s);                                             53
**
*/

/*
**> lnpdfn
**
**  Purpose:  Computes normal log-probabilities.
**
**  Format:   z = lnpdfn(x);
**
**  Input:    x    NxK matrix, data.
**
**  Output:   z    NxK matrix, log-probabilities.
**
**  Remarks:  This does not compute the log of the joint Normal
**            pdf.  Instead, the scalar Normal density function
**            is computed element by element.
**
**            For multivariate probabilities with covariance matrix
**            see lnpdfmvn.
**
*/

proc lnpdfn(x);
    retp( -.918938533204672741 - (x.*x) / 2 );
endp;


/*
**> lnpdfmvn
**
**  Purpose:   Computes multivariate normal log-probabilities.
**
**  Format:   z = lnpdfmvn(x,s);
**
**  Input:    x    NxK matrix, data.
**
**            s    KxK matrix, covariance matrix.
**
**  Output:   z    Nx1 vector, log-probabilities.
**
*/


proc lnpdfmvn(x,s);
      local si, oldt, lndet;

      if cols(x) /= cols(s);
          if not trapchk(1);
              errorlog "ERROR:  covariance matrix not conformable";
              end;
          else;
              retp(error(0));
          endif;
      endif;

      if cols(x) == 1;
           retp(-.918938533204672741 - 0.5 * ( ln(s) + x.*x / s ) );
      endif;

      if rows(s) /= cols(s);
          if not trapchk(1);
              errorlog "ERROR:  covariance matrix not square";
              end;
          else;
              retp(error(0));
          endif;
      endif;

      oldt = trapchk(1);
      trap 1,1;
      si = invpd(s);
      trap oldt,1;
      if scalmiss(si);
          if not trapchk(1);
             errorlog "ERROR: covariance matrix not positive definite";
             end;
          else;
              retp(error(0));
          endif;
      endif;

      lndet = ln(detl);

    retp(-0.5*(rows(s)*1.83787706640934548 + lndet + sumc(((x*si).*x)')));
endp;