Poisson distribution

 BBS: Inland Empire Archive
Date: 03-22-93 (13:16)             Number: 252
From: DIK COATES                   Refer#: NONE
  To: ALL                           Recvd: NO  
Subj: Poisson distribution           Conf: (2) Quik_Bas

Just goofin' around for the last couple of days... putting together a library
of curve fitting utilities... routine is a little kludgy for checking y>0 and
x>=0 and being an integer value... original approach was to use

if INT(a!(c%,1))+.499999 = CLNG(a!(c%,1)) THEN

as a test for x being an integer.  Is there anyone that has a better way of
determining if a value is an INTEGER (or LONG) value?  Have included the
source for the routine affected in case anyone else is interested... Have
used Sigma in comments... If BBS filters ASCII, subst value is shown...
Also, if anyone has a formula for the coef of determination, would be much
appreciated... the math is giving me real heartburn...

Entered into Public Domain by R.A. Coates, Mar 22/93


'*********************************************************** FUNCTION Poisson
' The function returns a curve that approximates the data that would
' accompany a distribution of random events.  This is represented by the
' Poisson distribution using an equation of the form:  y = (ab^x)/x!
'
' CALL: temp% = Poisson% (dt(), cd, cdc, a, b)
'
' ARG:  dt()  - floating point array of data
'       cd    - coefficient of determination
'       cdc   - corrected coef. of determination
'       a     - coef. solution to the curve approximation
'       b     - coef. solution to the curve approximation
'
' RET:  TRUE if successful, otherwise FALSE
'
' REV:  93-03-21
'
FUNCTION Poisson% (dt!(), cd!, cdc!, a!, b!)
  ns% = UBOUND(dt!, 1)

  IF ns% < 3 THEN
    EXIT FUNCTION
  END IF

  REDIM dtt%(ns%)                   'transformed x-values into integer values

  FOR c% = 1 TO ns%                                             'exit if y<=0
    IF dt!(c%, 2) <= 0 THEN
      EXIT FUNCTION
    END IF

    IF dt!(c%, 1) < 0 OR dt! > 170 THEN                 'exit if x<0 or x>170
      EXIT FUNCTION                                    'if >170 DBLE overflow
    ELSE
      dtt%(c%) = CINT(dt!(c%, 1))
    END IF
  NEXT c%

  Poisson% = TRUE

  FOR c% = 1 TO ns%
    sumx# = sumx# + dtt%(c%)                                              'Ex
    sumx2# = sumx2# + dtt%(c%) * dtt%(c%)                               'Ex^2
    dummy% = FACTORIAL%(dtt%(c%), value#)
    tempx# = LOG(value#)                                              'ln(x!)
    sumlfx# = sumlfx# + tempx#                                       'Eln(x!)
    sumxlfx# = sumxlfx# + tempx# * dtt%(c%)                      'E(x*ln(x!))
    tempy# = LOG(dt!(c%, 2))                                           'ln(y)
    sumly# = sumly# + tempy#                                          'Eln(y)
    sumxly# = sumxly# + tempy# * dtt%(c%)                         'E(x*ln(y))
  NEXT c%
                            'E -is the same as the Greek summation sign Sigma

  temp1# = ns% * sumx2# - sumx# * sumx#
  temp2# = sumlfx# + sumly#
  temp3# = sumxlfx# + sumxly#
  pwra# = (temp2# * sumx2# - temp3# * sumx#) / temp1#
  a! = EXP(pwra#)
  pwrb# = (ns% * temp3# - temp2# * sumx#) / temp1#
  b! = EXP(pwrb#)
  cd! = 0
  cdc! = 0
END FUNCTION 'Poisson%()


'********************************************************* FUNCTION FACTORIAL
' The function calculates the factorial of an integer value.  On success, it
' returns TRUE and on failure, it returns FALSE.  The integer value must be
' greater than 0 and less than 171.  The value returned is a double precision
' variable.
'
' CALL: temp% = FACTORIAL% (n, value)
'
' ARG:  n     - integer value of the factorial number (n!)
'       value - calculated value of the factorial, returned as a double
'               precision number.
'
' RET:  TRUE if successful, otherwise FALSE
'
' REV:  93-03-21
'
FUNCTION FACTORIAL% (n%, value#)
  IF n% < 0 OR n% > 170 THEN
    FACTORIAL% = FALSE
    value# = 0
    EXIT FUNCTION
  ELSE
    FACTORIAL% = TRUE
    value# = 1

    IF n% THEN
      FOR c% = 1 TO n%
        value# = value# * c%
      NEXT c%
    END IF
  END IF
END FUNCTION 'FACTORIAL%()


Regards Dik

... Sign at gynecologist's office:  I'm at your cervix!

--- Maximus 2.01wb
 * Origin: Durham Systems (ONLINE!) (1:229/110)