LazStats: Better numerical stability of PoissonPDF.

git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7721 8e941d3f-bd1b-0410-a28a-d453659cc2b4

applications/lazstats/source/units/functionslib.pas

@@ -52,13 +52,14 @@ function KolmogorovProb(z: double): double;
 function KolmogorovTest(na: integer; const a: DblDyneVec; nb: integer;
   const b: DblDyneVec; option: String; AReport: TStrings): double;
 
-procedure poisson_cdf ( x : integer; a : double; VAR cdf : double );
+//procedure poisson_cdf ( x : integer; a : double; VAR cdf : double );
 procedure poisson_cdf_values (VAR n : integer; VAR a : double; VAR x : integer;
                               VAR fx : double );
 procedure poisson_cdf_inv (VAR cdf : double; VAR a : double; VAR x : integer );
 procedure poisson_check ( a : double );
 function factorial(x : integer) : integer;
-procedure poisson_pdf ( x : integer; VAR a : double; VAR pdf : double );
+
+//function PoissonPDF(x: integer; a: double): Double;
 
 
 implementation
@@ -1932,6 +1933,7 @@ begin
      result := prob;
 end;
 
+(*     wp: moved to MathUnit for easier testing
 
 procedure poisson_cdf ( x : integer; a : double; VAR cdf : double );
 VAR
@@ -1982,7 +1984,7 @@ begin
     cdf := sum2;
   end;
 end;
-
+   *)
 procedure poisson_cdf_values (VAR n : integer; VAR a : double; VAR x : integer;
                               VAR fx : double );
 VAR
@@ -2202,24 +2204,22 @@ begin
      ShowMessage('POISSON_CHECK - Fatal error.  A <= 0.');
 end;
 
-function factorial(x : integer) : longint; //integer;
-VAR
-   decx : longint; // integer;
-   product : longint; //integer;
+function Factorial(x: integer): longint; //integer;
+var
+  decx: longint; // integer;
+  product: longint; //integer;
 begin
-   decx := x;
-   product := 1;
-   while (decx > 0) do
-   begin
-        product := decx * product;
-        decx := decx - 1;
-   end;
-   result := product;
+  decx := x;
+  product := 1;
+  while (decx > 0) do
+  begin
+    product := decx * product;
+    decx := decx - 1;
+  end;
+  result := product;
 end;
 
-
-procedure poisson_pdf ( x : integer; VAR a : double; VAR pdf : double );
-begin
+(*    wp: moved to MathUnit for easier testing
 //
 //*******************************************************************************
 //
@@ -2261,11 +2261,14 @@ begin
 //
 //    Output, real PDF, the value of the PDF.
 //
-  if ( x < 0 ) then pdf := 0.0E+00
+function PoissonPDF(x: integer; a: double): Double;
+begin
+  if (x < 0) then
+    Result := 0.0
   else
-    pdf := exp ( - a ) * power(a,x) / factorial ( x );
+    Result := exp(-a) * power(a, x) / factorial(x);
 //      pdf := exp ( - a ) * power(a,x) / exp(logfactorial( x ));
 end;
-
+          *)
 end.
 

applications/lazstats/source/units/mathunit.pas

@@ -38,6 +38,10 @@ function CalcC4(n: Integer): Double;
 
 procedure MantisseAndExponent(x: Double; out Mantisse: Double; out Exponent: Integer);
 
+function FactorialLn(n: Integer): Double;
+
+function PoissonPDF(n: integer; a: double): Double;
+function PoissonCDF(n: Integer; a: double): Double;
 
 implementation
 
@@ -365,5 +369,113 @@ begin
 end;
 
 
+var
+  FactLnArray: array[1..100] of Double;
+
+procedure InitFactLn;
+var
+  i: Integer;
+begin
+  for i := Low(FactLnArray) to High(FactLnArray) do FactLnArray[i] := -1.0;
+end;
+
+{ Returns ln(n!) }
+function FactorialLn(n: Integer): Double;
+begin
+  if n < 0 then
+    raise Exception.Create('Negative factorial.');
+
+  if n <= 99 then
+  begin
+    if FactLnArray[n+1] < 0.0 then
+      FactLnArray[n+1] := GammaLn(n + 1.0);
+    Result := FactLnArray[n+1];
+  end else
+    Result := GammaLn(n + 1.0);
+end;
+
+
+{ POISSON_PDF evaluates the Poisson probability distribution function (PDF).
+
+  Formula:
+     PDF(n, A) = EXP (- A) * A*^n / n!
+
+  Discussion:
+    PDF(n, A) is the probability that the number of events observed
+    in a unit time period will be n, given the expected number
+    of events in a unit time.
+
+    The parameter A is the expected number of events per unit time.
+
+    The Poisson PDF is a discrete version of the Exponential PDF.
+
+    The time interval between two Poisson events is a random
+    variable with the Exponential PDF.
+
+  Modified:
+    01 February 1999
+
+  Author:
+    John Burkardt
+
+  Parameters:
+    Input, integer n, the argument of the PDF:  0 <= n
+    Input, real A, the parameter of the PDF.:   0.0E+00 < A.
+
+  Output, real PDF, the value of the PDF.
+}
+function PoissonPDF(n: integer; a: double): Double;
+// wp: modified for better numerical stability by calculating with the logs
+var
+  arg: Double;
+begin
+  if n < 0 then
+    raise exception.Create('Negative argument in PoissonCDF');
+  arg := -a + n * ln(a) - FactorialLn(n);
+  Result := exp(arg);
+end;
+
+
+{  POISSON_CDF evaluates the Poisson cumulative distribution function (CDF)
+
+   Definition:
+     CDF(X,A) is the probability that the number of events observed
+     in a unit time period will be no greater than X, given that the
+     expected number of events in a unit time period is A.
+
+  Modified:
+    28 January 1999
+
+  Author:
+    John Burkardt
+
+  Parameters:
+    Input, integer N, the argument of the CDF.  N >= 0.
+    Input, real A, the parameter of the PDF. 0.0 < A.
+    Output, real CDF, the value of the CDF.
+}
+function PoissonCDF(n: integer; a: double): Double;
+var
+  i: integer;
+  last, new1, sum2: double;
+begin
+  if (n < 0) then
+    raise Exception.Create('Negative argument in PoissonCDF');
+
+  new1 := exp(-a);
+  sum2 := new1;
+  for i := 1 to n do
+  begin
+    last := new1;
+    new1 := last * a /  i ;
+    sum2 := sum2 + new1;
+  end;
+  Result := sum2;
+end;
+
+
+initialization
+  InitFactLn();
+
 end.