rndp.src

来自「没有说明」· SRC 代码 · 共 78 行

/*
** rndp.src - random poisson numbers
**
**
** (C) Copyright 1988-1995 by 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.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**> rndp
**
**  Format      y = rndp(r,c,l);
**
**  Input       r   scalar, number of rows.
**              c   scalar, number of columns.
**
**              l   expected value(s) of Poisson variates.
**                  if 1x1 scalar, used for all variates.
**                  if rx1 vector, row of l used for each row of result.
**                  if cx1 vector, col of l used for each col of result.
**                  if rxc matrix, element of l used for each element of result.
**       
**   Output     y   RxC matrix of random poisson numbers with parameters l.
*/


proc rndp(r,c,l);
    local i,j,z,l0;
    z = zeros(r,c);
    i = 1;
    do until i > r;
        j = 1;
        do until j > c;
            if cols(l) == c and rows(l) == r;
                l0 = l[i,j];
            elseif rows(l) == r;
                l0 = l[i,.];
            elseif cols(l) == c;
                l0 = l[.,j];
            else;
                l0 = l[1,1];
            endif;
            z[i,j] = _rndp(l0);
            j = j + 1;
       endo;
       i = i + 1;
    endo;
    retp(z);
endp;

proc _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;