statistics: Bug #7569: increased accuracy of inversion of several CDF 21/1921/4
Michael Baudin [Tue, 14 Sep 2010 16:06:59 +0000 (18:06 +0200)]
Change-Id: I1c64bd7ca81e72f98c4cf83372c7e862b8abfe69

18 files changed:
scilab/CHANGES_5.3.X
scilab/RELEASE_NOTES_5.3.X
scilab/modules/statistics/src/dcdflib/cdfbet.f
scilab/modules/statistics/src/dcdflib/cdfbin.f
scilab/modules/statistics/src/dcdflib/cdfchi.f
scilab/modules/statistics/src/dcdflib/cdfchn.f
scilab/modules/statistics/src/dcdflib/cdff.f
scilab/modules/statistics/src/dcdflib/cdffnc.f
scilab/modules/statistics/src/dcdflib/cdfgam.f
scilab/modules/statistics/src/dcdflib/cdfnbn.f
scilab/modules/statistics/src/dcdflib/cdfpoi.f
scilab/modules/statistics/src/dcdflib/cdft.f
scilab/modules/statistics/tests/unit_tests/cdfgam.dia.ref
scilab/modules/statistics/tests/unit_tests/cdfgam.tst
scilab/modules/statistics/tests/unit_tests/cdfnor.dia.ref
scilab/modules/statistics/tests/unit_tests/cdfnor.tst
scilab/modules/statistics/tests/unit_tests/cdfpoi.dia.ref
scilab/modules/statistics/tests/unit_tests/cdfpoi.tst

index 417725e..b2fb06e 100644 (file)
@@ -104,6 +104,10 @@ Bug Fixes:
 
 * bug 7231 fixed - mtlb_num2str did not manage second input argument.
 
+* bug 7569 fixed - the number of accurate digits during inversion of cdfbet, 
+                   cdfgam, cdfbin, cdfchi, cdfchin, cdff, cdffnc, cdfnbn, cdfpoi 
+                   was only 8.
+
 * bug 7640 fixed - xs2pdf, xs2eps, xs2emf crashed if filename prefix had less
                    of three characters.
 
index a2b3199..f41927d 100644 (file)
@@ -17,5 +17,8 @@
 - Under Linux, the binary requires at least the GLIBC version 2.7. This should
   be fixed in the beta 3.
 
+- The accuracy of inversion of most distribution functions is increased: see
+  bug #7569 for details.
+
 Besides these points, do not hesitate to report bugs on:
 http://bugzilla.scilab.org/
index a50da86..a886de6 100644 (file)
@@ -101,7 +101,7 @@ C
 C**********************************************************************
 C     .. Parameters ..
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION zero,inf
index e728e69..1c8c17a 100644 (file)
@@ -95,7 +95,7 @@ C     .. Parameters ..
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION zero,inf
       PARAMETER (zero=1.0D-300,inf=1.0D300)
       DOUBLE PRECISION one
index 68f8818..b0325e5 100644 (file)
@@ -81,7 +81,7 @@ C
 C**********************************************************************
 C     .. Parameters ..
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION zero,inf
index ee4c71c..9ed7225 100644 (file)
@@ -95,7 +95,7 @@ C     .. Parameters ..
       DOUBLE PRECISION tent4
       PARAMETER (tent4=1.0D4)
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION zero,one,inf
index 17fca81..90435cd 100644 (file)
@@ -90,7 +90,7 @@ C
 C**********************************************************************
 C     .. Parameters ..
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION zero,inf
index 57f62cf..e4d7cee 100644 (file)
@@ -106,7 +106,7 @@ C     .. Parameters ..
       DOUBLE PRECISION tent4
       PARAMETER (tent4=1.0D4)
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION zero,one,inf
index aab74e9..d690bfa 100644 (file)
@@ -103,7 +103,7 @@ C     variable. See below comments "CSS".
 C**********************************************************************
 C     .. Parameters ..
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION zero,inf
index bc8644f..2dd86ac 100644 (file)
@@ -102,7 +102,7 @@ C
 C**********************************************************************
 C     .. Parameters ..
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION inf
index cab3bd3..cacbf5c 100644 (file)
@@ -79,7 +79,7 @@ C
 C**********************************************************************
 C     .. Parameters ..
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION inf
index ddeb9d6..26c6475 100644 (file)
@@ -77,7 +77,7 @@ C
 C**********************************************************************
 C     .. Parameters ..
       DOUBLE PRECISION tol
-      PARAMETER (tol=1.0D-8)
+      PARAMETER (tol=1.0D-13)
       DOUBLE PRECISION atol
       PARAMETER (atol=1.0D-50)
       DOUBLE PRECISION zero,inf
index 826d481..69eecc3 100644 (file)
 // =============================================================================
-// Tests for cdfgam() function
-//
-// Scilab Team
-// Copyright INRIA
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Michael Baudin
 //
+//  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-prec  = 1.e-5;
-shape = 2;
-scale = 3;
-bn    = 1;
-deff('[y]=Gamma(x)','y=bn*(x^(shape-1) * exp(-scale*x))');
-bn     = intg(0,10^3,Gamma);
-bn     = 1/bn;
-x      = 1:10;
-[P,Q]  = cdfgam("PQ",x,shape*ones(x),scale*ones(x));
-P1     = [];
-for xx=x
-       P1=[P1,intg(0,xx,Gamma)];
+// <-- JVM NOT MANDATORY -->
+// <-- ENGLISH IMPOSED -->
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+function d = assert_computedigits ( computed , expected )
+  nre = size(expected,"r")
+  nce = size(expected,"c")
+  // Update shape
+  expected = expected (:)
+  computed = computed (:)
+  //
+  dmin = 0
+  dmax = -log10(2^(-53))
+  //
+  d = zeros(expected)
+  //
+  n = size(expected,"*")
+  for i = 1 : n
+    if ( isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmax
+    elseif ( isnan(expected(i)) & ~isnan(computed(i)) ) then
+      d(i) = dmin
+    elseif ( ~isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmin
+      // From now, both expected and computed are non-nan
+    elseif ( expected(i) == 0 & computed(i) == 0 ) then
+      d(i) = dmax
+    elseif ( expected(i) == 0 & computed(i) <> 0 ) then
+      d(i) = dmin
+      // From now, expected(i) is non-zero
+    elseif ( expected(i) == computed(i) ) then
+      d(i) = dmax
+      // From now, expected and computed are different
+    elseif ( expected(i) == %inf & computed(i) <> %inf ) then
+      d(i) = dmin
+    elseif ( expected(i) == -%inf & computed(i) <> -%inf ) then
+      d(i) = dmin
+      // From now, neither of computed, nor expected is infinity
+    else
+      d(i) = max ( -log10 ( abs(computed(i)-expected(i)) / abs(expected(i)) ) , dmin )
+    end
+  end
+  //
+  // Reshape
+  d = matrix(d,nre,nce)
+endfunction
+//
+// Assessing the quality of the Normal distribution function
+// References
+//   Yalta, A. T. 2008. The accuracy of statistical distributions in Microsoft®Excel 2007. Comput. Stat. Data Anal. 52, 10 (Jun. 2008), 4579-4586. DOI= http://dx.doi.org/10.1016/j.csda.2008.03.005 
+//   Computation of Statistical Distributions (ELV), Leo Knüsel 
+// Table 5
+// Check Gamma distribution with parameters (x, alpha, beta = 1, Sigma = 1)
+//
+// Table of inputs from Yalta, 2008
+// [x shape scale P ]
+table = [
+ 0.1 , 0.1 , 1 , 0.827552 
+ 0.2 , 0.1 , 1 , 0.879420 
+ 0.2 , 0.2 , 1 , 0.764435 
+ 0.3 , 0.2 , 1 , 0.816527 
+ 0.3 , 0.3 , 1 , 0.726957 
+ 0.4 , 0.3 , 1 , 0.776381 
+ 0.4 , 0.4 , 1 , 0.701441 
+ 0.5 , 0.4 , 1 , 0.748019 
+ 0.5 , 0.5 , 1 , 0.682689 
+ 0.6 , 0.5 , 1 , 0.726678 
+];
+precision = 1.e-5;
+ntests = size(table,"r");
+for i = 1 : ntests
+  x = table(i,1);
+  shape = table(i,2);
+  scale = table(i,3);
+  expected = table(i,4);
+  // Caution: this is the rate !
+  rate = 1/scale;
+  [computed,Q]=cdfgam("PQ",x,shape,rate);
+  assert_close ( computed , expected , precision );
+  assert_close ( Q , 1 - expected , precision );
+end
+// Table of inputs computed from R-2.8.1
+// [x shape scale PDF-P CDF-P CDF-Q]
+table = [
+1.000000000000000056D-01 1.000000000000000056D-01 1.000000000000000000D+00 7.554920138253073958D-01 8.275517595858505882D-01 1.724482404141494951D-01
+2.000000000000000111D-01 1.000000000000000056D-01 1.000000000000000000D+00 3.663307993056703071D-01 8.794196267900568076D-01 1.205803732099432063D-01
+2.000000000000000111D-01 2.000000000000000111D-01 1.000000000000000000D+00 6.462857778271943188D-01 7.644345975029189777D-01 2.355654024970809945D-01
+2.999999999999999889D-01 2.000000000000000111D-01 1.000000000000000000D+00 4.227875047602157044D-01 8.165267943336527168D-01 1.834732056663473110D-01
+2.999999999999999889D-01 2.999999999999999889D-01 1.000000000000000000D+00 5.752117576599179438D-01 7.269573437103662439D-01 2.730426562896338116D-01
+4.000000000000000222D-01 2.999999999999999889D-01 1.000000000000000000D+00 4.255407854753925911D-01 7.763805810166358734D-01 2.236194189833642099D-01
+4.000000000000000222D-01 4.000000000000000222D-01 1.000000000000000000D+00 5.236648604477927016D-01 7.014412706419403953D-01 2.985587293580597157D-01
+5.000000000000000000D-01 4.000000000000000222D-01 1.000000000000000000D+00 4.144555659263016167D-01 7.480185547260104206D-01 2.519814452739895239D-01
+5.000000000000000000D-01 5.000000000000000000D-01 1.000000000000000000D+00 4.839414490382866751D-01 6.826894921370858516D-01 3.173105078629140929D-01
+5.999999999999999778D-01 5.000000000000000000D-01 1.000000000000000000D+00 3.997355278034666060D-01 7.266783217077018575D-01 2.733216782922981980D-01
+5.000000000000000000D-01 5.000000000000000000D-01 2.000000000000000000D+00 4.393912894677223790D-01 5.204998778130465187D-01 4.795001221869534258D-01
+5.000000000000000000D-01 5.000000000000000000D-01 3.000000000000000000D+00 3.899393114454822729D-01 4.362971383492270094D-01 5.637028616507729906D-01
+5.000000000000000000D-01 5.000000000000000000D-01 4.000000000000000000D+00 3.520653267642995243D-01 3.829249225480261809D-01 6.170750774519737636D-01
+1.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 2.075537487102973866D-01 8.427007929497148941D-01 1.572992070502851891D-01
+2.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 5.399096651318804896D-02 9.544997361036415828D-01 4.550026389635838248D-02
+4.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 5.166746338523012412D-03 9.953222650189527121D-01 4.677734981047261889D-03
+1.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 8.099910956089122777D-06 9.999922557835689840D-01 7.744216431044085842D-06
+2.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 2.600281868827196957D-10 9.999999997460371493D-01 2.539628589470869077D-10
+4.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 3.789795640412981196D-19 1.000000000000000000D+00 3.744097384202895045D-19
+//1.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 2.098828115677222045D-45 1.000000000000000000D+00 2.088487583762558879D-45
+3.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 1.67694904029982009D-132 1.000000000000000000D+00 1.67416798469182012D-132
+//5.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 1.79762504374667411D-219 1.000000000000000000D+00 1.79583278480075297D-219
+//1.000000000000000000D+03 5.000000000000000000D-01 1.000000000000000000D+00 0.000000000000000000D+00 1.000000000000000000D+00 0.000000000000000000D+00
+1.000000000000000021D-02 5.000000000000000000D-01 1.000000000000000000D+00 5.585758033944684620D+00 1.124629160182848975D-01 8.875370839817151580D-01
+1.000000000000000048D-04 5.000000000000000000D-01 1.000000000000000000D+00 5.641331674102550409D+01 1.128341555584961957D-02 9.887165844441503371D-01
+1.000000000000000021D-08 5.000000000000000000D-01 1.000000000000000000D+00 5.641895779058606422D+03 1.128379163334249004D-04 9.998871620836665697D-01
+9.999999999999999452D-21 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477570534D+09 1.128379167095512970D-10 9.999999998871620388D-01
+9.999999999999999293D-41 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477568717D+19 1.128379167095512972D-20 1.000000000000000000D+00
+1.00000000000000002D-100 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477541988D+49 1.128379167095513082D-50 1.000000000000000000D+00
+9.99999999999999982D-201 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477511468D+99 1.12837916709551300D-100 1.000000000000000000D+00
+//1.00000000000000002D-300 5.000000000000000000D-01 1.000000000000000000D+00 5.64189583547731891D+149 1.12837916709551298D-150 1.000000000000000000D+00
+//0.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00                     %inf 0.000000000000000000D+00 1.000000000000000000D+00
+];
+// For the inversion of Shape, require only 8 digits, as 
+// a consequence of bug #7569: http://bugzilla.scilab.org/show_bug.cgi?id=7569
+//
+// Some tests do not pass:
+// http://bugzilla.scilab.org/show_bug.cgi?id=8031
+// http://bugzilla.scilab.org/show_bug.cgi?id=8030
+//
+// Prints the number of accurate digits.
+precision = 1.e-12;
+precinverse = 1.e-8;
+ntests = size(table,"r");
+for i = 1 : ntests
+  x = table(i,1);
+  shape = table(i,2);
+  scale = table(i,3);
+  p = table(i,5);
+  q = table(i,6);
+  // Caution: this is the rate !
+  rate = 1/scale;
+  [p1,q1]=cdfgam("PQ",x,shape,rate);
+  x1     = cdfgam("X",shape,rate,p,q);
+  shape1 = cdfgam("Shape",rate,p,q,x);
+  rate1 = cdfgam("Scale",p,q,x,shape);
+  if ( %t ) then
+    assert_close ( p1 , p , precision );
+    assert_close ( q1 , q , precision );
+    assert_close ( x1 , x , precision );
+    assert_close ( shape1 , shape , precinverse );
+    assert_close ( rate1 , rate , precinverse );
+  end
+  if ( %f ) then
+    dp = assert_computedigits ( p1 , p );
+    dq = assert_computedigits ( q1 , q );
+    dx = assert_computedigits ( x1 , x );
+    ds = assert_computedigits ( shape1 , shape );
+    dr = assert_computedigits ( rate1 , rate );
+    mprintf("Test #%3d/%3d: Digits p1= %.1f, q1=%.1f, X= %.1f, S= %.1f, R= %.1f\n",i,ntests,dp,dq,dx,ds,dr);
+  end
 end
-if norm(P1-P) > prec then bugmes();quit;end
-shape    = 2*x;
-scale    = 3*x;
-[P,Q]    = cdfgam("PQ",x,shape,scale);
-[x1]     = cdfgam("X",shape,scale,P,Q);
-[shape1] = cdfgam("Shape",scale,P,Q,x);
-[scale1] = cdfgam("Scale",P,Q,x,shape);
-if norm(shape1-shape) > 10*prec then bugmes();quit;end
-if norm(scale1-scale) > 10*prec then bugmes();quit;end
-if norm(x1-x)         > prec    then bugmes();quit;end
index 3058061..9c4d1f7 100644 (file)
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - INRIA - Michael Baudin
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
 
-// =============================================================================
-// Tests for cdfgam() function
-// =============================================================================
 
-prec  = 1.e-5;
+// <-- JVM NOT MANDATORY -->
+// <-- ENGLISH IMPOSED -->
 
-shape = 2;
-scale = 3;
-bn    = 1;
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
 
-deff('[y]=Gamma(x)','y=bn*(x^(shape-1) * exp(-scale*x))');
+function d = assert_computedigits ( computed , expected )
+  nre = size(expected,"r")
+  nce = size(expected,"c")
+  // Update shape
+  expected = expected (:)
+  computed = computed (:)
+  //
+  dmin = 0
+  dmax = -log10(2^(-53))
+  //
+  d = zeros(expected)
+  //
+  n = size(expected,"*")
+  for i = 1 : n
+    if ( isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmax
+    elseif ( isnan(expected(i)) & ~isnan(computed(i)) ) then
+      d(i) = dmin
+    elseif ( ~isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmin
+      // From now, both expected and computed are non-nan
+    elseif ( expected(i) == 0 & computed(i) == 0 ) then
+      d(i) = dmax
+    elseif ( expected(i) == 0 & computed(i) <> 0 ) then
+      d(i) = dmin
+      // From now, expected(i) is non-zero
+    elseif ( expected(i) == computed(i) ) then
+      d(i) = dmax
+      // From now, expected and computed are different
+    elseif ( expected(i) == %inf & computed(i) <> %inf ) then
+      d(i) = dmin
+    elseif ( expected(i) == -%inf & computed(i) <> -%inf ) then
+      d(i) = dmin
+      // From now, neither of computed, nor expected is infinity
+    else
+      d(i) = max ( -log10 ( abs(computed(i)-expected(i)) / abs(expected(i)) ) , dmin )
+    end
+  end
+  //
+  // Reshape
+  d = matrix(d,nre,nce)
+endfunction
 
-bn     = intg(0,10^3,Gamma);
-bn     = 1/bn;
-x      = 1:10;
-[P,Q]  = cdfgam("PQ",x,shape*ones(x),scale*ones(x));
-P1     = [];
+//
+// Assessing the quality of the Normal distribution function
+// References
+//   Yalta, A. T. 2008. The accuracy of statistical distributions in Microsoft®Excel 2007. Comput. Stat. Data Anal. 52, 10 (Jun. 2008), 4579-4586. DOI= http://dx.doi.org/10.1016/j.csda.2008.03.005 
+//   Computation of Statistical Distributions (ELV), Leo Knüsel 
+// Table 5
+// Check Gamma distribution with parameters (x, alpha, beta = 1, Sigma = 1)
+//
 
-for xx=x
-       P1=[P1,intg(0,xx,Gamma)];
+
+// Table of inputs from Yalta, 2008
+// [x shape scale P ]
+table = [
+ 0.1 , 0.1 , 1 , 0.827552 
+ 0.2 , 0.1 , 1 , 0.879420 
+ 0.2 , 0.2 , 1 , 0.764435 
+ 0.3 , 0.2 , 1 , 0.816527 
+ 0.3 , 0.3 , 1 , 0.726957 
+ 0.4 , 0.3 , 1 , 0.776381 
+ 0.4 , 0.4 , 1 , 0.701441 
+ 0.5 , 0.4 , 1 , 0.748019 
+ 0.5 , 0.5 , 1 , 0.682689 
+ 0.6 , 0.5 , 1 , 0.726678 
+];
+
+precision = 1.e-5;
+ntests = size(table,"r");
+for i = 1 : ntests
+  x = table(i,1);
+  shape = table(i,2);
+  scale = table(i,3);
+  expected = table(i,4);
+  // Caution: this is the rate !
+  rate = 1/scale;
+  [computed,Q]=cdfgam("PQ",x,shape,rate);
+  assert_close ( computed , expected , precision );
+  assert_close ( Q , 1 - expected , precision );
 end
 
-if norm(P1-P) > prec then pause,end
+// Table of inputs computed from R-2.8.1
+// [x shape scale PDF-P CDF-P CDF-Q]
+table = [
+1.000000000000000056D-01 1.000000000000000056D-01 1.000000000000000000D+00 7.554920138253073958D-01 8.275517595858505882D-01 1.724482404141494951D-01
+2.000000000000000111D-01 1.000000000000000056D-01 1.000000000000000000D+00 3.663307993056703071D-01 8.794196267900568076D-01 1.205803732099432063D-01
+2.000000000000000111D-01 2.000000000000000111D-01 1.000000000000000000D+00 6.462857778271943188D-01 7.644345975029189777D-01 2.355654024970809945D-01
+2.999999999999999889D-01 2.000000000000000111D-01 1.000000000000000000D+00 4.227875047602157044D-01 8.165267943336527168D-01 1.834732056663473110D-01
+2.999999999999999889D-01 2.999999999999999889D-01 1.000000000000000000D+00 5.752117576599179438D-01 7.269573437103662439D-01 2.730426562896338116D-01
+4.000000000000000222D-01 2.999999999999999889D-01 1.000000000000000000D+00 4.255407854753925911D-01 7.763805810166358734D-01 2.236194189833642099D-01
+4.000000000000000222D-01 4.000000000000000222D-01 1.000000000000000000D+00 5.236648604477927016D-01 7.014412706419403953D-01 2.985587293580597157D-01
+5.000000000000000000D-01 4.000000000000000222D-01 1.000000000000000000D+00 4.144555659263016167D-01 7.480185547260104206D-01 2.519814452739895239D-01
+5.000000000000000000D-01 5.000000000000000000D-01 1.000000000000000000D+00 4.839414490382866751D-01 6.826894921370858516D-01 3.173105078629140929D-01
+5.999999999999999778D-01 5.000000000000000000D-01 1.000000000000000000D+00 3.997355278034666060D-01 7.266783217077018575D-01 2.733216782922981980D-01
+5.000000000000000000D-01 5.000000000000000000D-01 2.000000000000000000D+00 4.393912894677223790D-01 5.204998778130465187D-01 4.795001221869534258D-01
+5.000000000000000000D-01 5.000000000000000000D-01 3.000000000000000000D+00 3.899393114454822729D-01 4.362971383492270094D-01 5.637028616507729906D-01
+5.000000000000000000D-01 5.000000000000000000D-01 4.000000000000000000D+00 3.520653267642995243D-01 3.829249225480261809D-01 6.170750774519737636D-01
+1.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 2.075537487102973866D-01 8.427007929497148941D-01 1.572992070502851891D-01
+2.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 5.399096651318804896D-02 9.544997361036415828D-01 4.550026389635838248D-02
+4.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 5.166746338523012412D-03 9.953222650189527121D-01 4.677734981047261889D-03
+1.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 8.099910956089122777D-06 9.999922557835689840D-01 7.744216431044085842D-06
+2.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 2.600281868827196957D-10 9.999999997460371493D-01 2.539628589470869077D-10
+4.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 3.789795640412981196D-19 1.000000000000000000D+00 3.744097384202895045D-19
+//1.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 2.098828115677222045D-45 1.000000000000000000D+00 2.088487583762558879D-45
+3.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 1.67694904029982009D-132 1.000000000000000000D+00 1.67416798469182012D-132
+//5.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 1.79762504374667411D-219 1.000000000000000000D+00 1.79583278480075297D-219
+//1.000000000000000000D+03 5.000000000000000000D-01 1.000000000000000000D+00 0.000000000000000000D+00 1.000000000000000000D+00 0.000000000000000000D+00
+1.000000000000000021D-02 5.000000000000000000D-01 1.000000000000000000D+00 5.585758033944684620D+00 1.124629160182848975D-01 8.875370839817151580D-01
+1.000000000000000048D-04 5.000000000000000000D-01 1.000000000000000000D+00 5.641331674102550409D+01 1.128341555584961957D-02 9.887165844441503371D-01
+1.000000000000000021D-08 5.000000000000000000D-01 1.000000000000000000D+00 5.641895779058606422D+03 1.128379163334249004D-04 9.998871620836665697D-01
+9.999999999999999452D-21 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477570534D+09 1.128379167095512970D-10 9.999999998871620388D-01
+9.999999999999999293D-41 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477568717D+19 1.128379167095512972D-20 1.000000000000000000D+00
+1.00000000000000002D-100 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477541988D+49 1.128379167095513082D-50 1.000000000000000000D+00
+9.99999999999999982D-201 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477511468D+99 1.12837916709551300D-100 1.000000000000000000D+00
+//1.00000000000000002D-300 5.000000000000000000D-01 1.000000000000000000D+00 5.64189583547731891D+149 1.12837916709551298D-150 1.000000000000000000D+00
+//0.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00                     %inf 0.000000000000000000D+00 1.000000000000000000D+00
+];
 
-shape    = 2*x;
-scale    = 3*x;
-[P,Q]    = cdfgam("PQ",x,shape,scale);
-[x1]     = cdfgam("X",shape,scale,P,Q);
-[shape1] = cdfgam("Shape",scale,P,Q,x);
-[scale1] = cdfgam("Scale",P,Q,x,shape);
+// For the inversion of Shape, require only 8 digits, as 
+// a consequence of bug #7569: http://bugzilla.scilab.org/show_bug.cgi?id=7569
+//
+// Some tests do not pass:
+// http://bugzilla.scilab.org/show_bug.cgi?id=8031
+// http://bugzilla.scilab.org/show_bug.cgi?id=8030
+//
+// Prints the number of accurate digits.
+
+precision = 1.e-12;
+precinverse = 1.e-8;
+
+ntests = size(table,"r");
+for i = 1 : ntests
+  x = table(i,1);
+  shape = table(i,2);
+  scale = table(i,3);
+  p = table(i,5);
+  q = table(i,6);
+  // Caution: this is the rate !
+  rate = 1/scale;
+  [p1,q1]=cdfgam("PQ",x,shape,rate);
+  x1     = cdfgam("X",shape,rate,p,q);
+  shape1 = cdfgam("Shape",rate,p,q,x);
+  rate1 = cdfgam("Scale",p,q,x,shape);
+  if ( %t ) then
+    assert_close ( p1 , p , precision );
+    assert_close ( q1 , q , precision );
+    assert_close ( x1 , x , precision );
+    assert_close ( shape1 , shape , precinverse );
+    assert_close ( rate1 , rate , precinverse );
+  end
+  if ( %f ) then
+    dp = assert_computedigits ( p1 , p );
+    dq = assert_computedigits ( q1 , q );
+    dx = assert_computedigits ( x1 , x );
+    ds = assert_computedigits ( shape1 , shape );
+    dr = assert_computedigits ( rate1 , rate );
+    mprintf("Test #%3d/%3d: Digits p1= %.1f, q1=%.1f, X= %.1f, S= %.1f, R= %.1f\n",i,ntests,dp,dq,dx,ds,dr);
+  end
+end
 
-if norm(shape1-shape) > 10*prec then pause,end
-if norm(scale1-scale) > 10*prec then pause,end
-if norm(x1-x)         > prec    then pause,end
index 0f3cbd5..ab3334e 100644 (file)
 // =============================================================================
-// Tests for cdfnor() function
-//
-// Scilab Team
-// Copyright INRIA
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Michael Baudin
 //
+//  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-prec  = 1.e-5;
-x     = 0:0.1:3;
-[p,q] = cdfnor("PQ",x,0.0*ones(x),1.0*ones(x));
-if norm( p - 1/2*(1+erf(x/sqrt(2))))> 100*%eps then bugmes();quit;end
-yy    = cdfnor("PQ",0.5,0.0,1.0);
-deff('y=f(t)','y=exp(-t^2/2)');
-if norm(yy-(1/2+ 1/sqrt(2*%pi)*intg(0,0.5,f)))> %eps then bugmes();quit;end
-[x1] = cdfnor("X",0.0*ones(x),1.0*ones(x),p,q);
-if norm(x-x1) > prec then bugmes();quit;end
-M     = 2*x;
-[p,q] = cdfnor("PQ",x,M,1.0*ones(x));
-M1    = cdfnor("Mean",1.0*ones(x),p,q,x);
-if norm(M-M1) > prec then bugmes();quit;end
-// avoid M=x
-M      = 2*x+1;
-Std    = x+1;
-[p,q]  = cdfnor("PQ",x,M,Std);
-Std1   = cdfnor("Std",p,q,x,M);
-if norm(Std-Std1) > prec then bugmes();quit;end
+// <-- JVM NOT MANDATORY -->
+// <-- ENGLISH IMPOSED -->
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+function d = assert_computedigits ( computed , expected )
+  nre = size(expected,"r")
+  nce = size(expected,"c")
+  // Update shape
+  expected = expected (:)
+  computed = computed (:)
+  //
+  dmin = 0
+  dmax = -log10(2^(-53))
+  //
+  d = zeros(expected)
+  //
+  n = size(expected,"*")
+  for i = 1 : n
+    if ( isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmax
+    elseif ( isnan(expected(i)) & ~isnan(computed(i)) ) then
+      d(i) = dmin
+    elseif ( ~isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmin
+      // From now, both expected and computed are non-nan
+    elseif ( expected(i) == 0 & computed(i) == 0 ) then
+      d(i) = dmax
+    elseif ( expected(i) == 0 & computed(i) <> 0 ) then
+      d(i) = dmin
+      // From now, expected(i) is non-zero
+    elseif ( expected(i) == computed(i) ) then
+      d(i) = dmax
+      // From now, expected and computed are different
+    elseif ( expected(i) == %inf & computed(i) <> %inf ) then
+      d(i) = dmin
+    elseif ( expected(i) == -%inf & computed(i) <> -%inf ) then
+      d(i) = dmin
+      // From now, neither of computed, nor expected is infinity
+    else
+      d(i) = max ( -log10 ( abs(computed(i)-expected(i)) / abs(expected(i)) ) , dmin )
+    end
+  end
+  //
+  // Reshape
+  d = matrix(d,nre,nce)
+endfunction
+// 
+// References
+// Assessing the quality of the Normal distribution function
+// References
+//   Yalta, A. T. 2008. The accuracy of statistical distributions in Microsoft®Excel 2007. Comput. Stat. Data Anal. 52, 10 (Jun. 2008), 4579-4586. DOI= http://dx.doi.org/10.1016/j.csda.2008.03.005 
+//   Computation of Statistical Distributions (ELV), Leo Knüsel 
+// Table 6
+// Check inverse of normal law
+//
+// Table from Yalta, 2008.
+table = [
+  0.5 , 0.0 
+  1.e-1 , -1.28155 
+  1.e-2 , -2.32635 
+  1.e-3 , -3.09023 
+  1.e-4 , -3.71902 
+  1.e-5 , -4.26489 
+  1.e-6 , -4.75342 
+  1.e-7 , -5.19934 
+  1.e-15 , -7.94135 
+  1.e-16 , -8.22208 
+  1.e-100 , -21.2735 
+  1.e-197 , -29.9763 
+  1.e-198 , -30.0529 
+  1.e-300 , -37.0471 
+];
+precision = 1.e-5;
+nt = size(table,"r");
+for k = 1 : nt
+  p = table(k,1);
+  expected = table(k,2);
+  q = 1 - p;
+  Mean = 0;
+  Std = 1;
+  computed = cdfnor ( "X" , Mean , Std , p , q );
+  assert_close ( computed , expected , precision );
+end
+//
+// Values from R-2.8.1
+// table = [x mu sigma PDF-P CDF-P CDF-Q]
+// Some tests do not pass with Scilab.
+//
+// See : http://bugzilla.scilab.org/show_bug.cgi?id=8032
+//
+// Prints the number of accurate digits.
+table = [
+//                    %inf 1.000000000000000000D+00 2.000000000000000000D+00 0.000000000000000000D+00 1.000000000000000000D+00 0.000000000000000000D+00
+//1.630146181031128094D+01 1.000000000000000000D+00 2.000000000000000000D+00 3.885547484725481156D-14 9.999999999999900080D-01 9.992007221626602924D-15
+1.372268177939483991D+01 1.000000000000000000D+00 2.000000000000000000D+00 3.255794261707647101D-10 9.999999998999999917D-01 1.000000082740384053D-10
+1.139867516458132002D+01 1.000000000000000000D+00 2.000000000000000000D+00 2.689766239123555948D-07 9.999999000000000526D-01 9.999999994736501079D-08
+9.529781587847679702D+00 1.000000000000000000D+00 2.000000000000000000D+00 2.239366490571142939D-05 9.999900000000000455D-01 9.999999999954530395D-06
+8.438032970911416797D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.979239833799469940D-04 9.999000000000000110D-01 9.999999999998900014D-05
+7.180464612335624608D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.683545038532002075D-03 9.989999999999999991D-01 1.000000000000003924D-03
+5.652695748081682403D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.332607110172901940D-02 9.899999999999999911D-01 1.000000000000000021D-02
+3.563103131089199849D+00 1.000000000000000000D+00 2.000000000000000000D+00 8.774916596624346421D-02 8.999999999999999112D-01 1.000000000000001027D-01
+2.683242467145829036D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.399809602039041034D-01 8.000000000000001554D-01 1.999999999999999001D-01
+1.506694206271600001D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.931712667484301871D-01 6.000000000000000888D-01 3.999999999999999112D-01
+//1.000000000000000000D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.994711402007164069D-01 5.000000000000000000D-01 5.000000000000000000D-01
+-6.832424671458280363D-01 1.000000000000000000D+00 2.000000000000000000D+00 1.399809602039041867D-01 2.000000000000000111D-01 8.000000000000000444D-01
+-1.563103131089200959D+00 1.000000000000000000D+00 2.000000000000000000D+00 8.774916596624339482D-02 1.000000000000000056D-01 9.000000000000000222D-01
+-3.652695748081681959D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.332607110172904022D-02 1.000000000000001062D-02 9.899999999999999911D-01
+-7.529781587845650215D+00 1.000000000000000000D+00 2.000000000000000000D+00 2.239366490580828930D-05 9.999999999999960160D-06 9.999900000000000455D-01
+-1.172268180480810962D+01 1.000000000000000000D+00 2.000000000000000000D+00 3.255793998537787058D-10 1.000000000000008954D-10 9.999999998999999917D-01
+-1.752468017959679969D+01 1.000000000000000000D+00 2.000000000000000000D+00 4.683961267403037244D-20 1.000000000000072017D-20 1.000000000000000000D+00
+-2.886667506957698137D+01 1.000000000000000000D+00 2.000000000000000000D+00 7.499857137115675432D-50 9.999999999999744883D-51 1.000000000000000000D+00
+//                   -%inf 1.000000000000000000D+00 2.000000000000000000D+00 0.000000000000000000D+00 0.000000000000000000D+00 1.000000000000000000D+00
+];
+precision = 1.e-12;
+precinv = 1.e-8;
+nt = size(table,"r");
+for k = 1 : nt
+  x = table(k,1);
+  mu = table(k,2);
+  std = table(k,3);
+  p = table(k,5);
+  q = table(k,6);
+  [ p1 , q1 ] = cdfnor("PQ",x,mu,std);
+  x1 = cdfnor("X",mu,std,p,q);
+  mu1 = cdfnor("Mean",std,p,q,x);
+  std1 = cdfnor("Std",p,q,x,mu);
+  if ( %t ) then
+    assert_close ( p1 , p , precision );
+    assert_close ( q1 , q , precision );
+    assert_close ( x1 , x , precinv );
+    assert_close ( mu1 , mu , precinv );
+    assert_close ( std1 , std , precinv );
+  end
+  if ( %f ) then
+    dP = assert_computedigits ( p1 , p );
+    dQ = assert_computedigits ( q1 , q );
+    dx = assert_computedigits ( x1 , x );
+    dmu = assert_computedigits ( mu1 , mu );
+    dstd = assert_computedigits ( std1 , std );
+    mprintf("Test #%3d/%3d: Digits p1= %.1f, q1=%.1f, X=%.1f, M=%.1f, S=%.1f\n",k,nt,dP,dQ,dx,dmu,dstd);
+  end
+end
index 2b5f8e2..5b9a689 100644 (file)
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - INRIA - Michael Baudin
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
 
-// =============================================================================
-// Tests for cdfnor() function
-// =============================================================================
 
-prec  = 1.e-5;
+// <-- JVM NOT MANDATORY -->
+// <-- ENGLISH IMPOSED -->
+
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+function d = assert_computedigits ( computed , expected )
+  nre = size(expected,"r")
+  nce = size(expected,"c")
+  // Update shape
+  expected = expected (:)
+  computed = computed (:)
+  //
+  dmin = 0
+  dmax = -log10(2^(-53))
+  //
+  d = zeros(expected)
+  //
+  n = size(expected,"*")
+  for i = 1 : n
+    if ( isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmax
+    elseif ( isnan(expected(i)) & ~isnan(computed(i)) ) then
+      d(i) = dmin
+    elseif ( ~isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmin
+      // From now, both expected and computed are non-nan
+    elseif ( expected(i) == 0 & computed(i) == 0 ) then
+      d(i) = dmax
+    elseif ( expected(i) == 0 & computed(i) <> 0 ) then
+      d(i) = dmin
+      // From now, expected(i) is non-zero
+    elseif ( expected(i) == computed(i) ) then
+      d(i) = dmax
+      // From now, expected and computed are different
+    elseif ( expected(i) == %inf & computed(i) <> %inf ) then
+      d(i) = dmin
+    elseif ( expected(i) == -%inf & computed(i) <> -%inf ) then
+      d(i) = dmin
+      // From now, neither of computed, nor expected is infinity
+    else
+      d(i) = max ( -log10 ( abs(computed(i)-expected(i)) / abs(expected(i)) ) , dmin )
+    end
+  end
+  //
+  // Reshape
+  d = matrix(d,nre,nce)
+endfunction
 
-x     = 0:0.1:3;
-[p,q] = cdfnor("PQ",x,0.0*ones(x),1.0*ones(x));
-if norm( p - 1/2*(1+erf(x/sqrt(2))))> 100*%eps then pause,end
-yy    = cdfnor("PQ",0.5,0.0,1.0);
+// 
+// References
+// Assessing the quality of the Normal distribution function
+// References
+//   Yalta, A. T. 2008. The accuracy of statistical distributions in Microsoft®Excel 2007. Comput. Stat. Data Anal. 52, 10 (Jun. 2008), 4579-4586. DOI= http://dx.doi.org/10.1016/j.csda.2008.03.005 
+//   Computation of Statistical Distributions (ELV), Leo Knüsel 
+// Table 6
+// Check inverse of normal law
+//
+  
+// Table from Yalta, 2008.
+table = [
+  0.5 , 0.0 
+  1.e-1 , -1.28155 
+  1.e-2 , -2.32635 
+  1.e-3 , -3.09023 
+  1.e-4 , -3.71902 
+  1.e-5 , -4.26489 
+  1.e-6 , -4.75342 
+  1.e-7 , -5.19934 
+  1.e-15 , -7.94135 
+  1.e-16 , -8.22208 
+  1.e-100 , -21.2735 
+  1.e-197 , -29.9763 
+  1.e-198 , -30.0529 
+  1.e-300 , -37.0471 
+];
+
+precision = 1.e-5;
+nt = size(table,"r");
+for k = 1 : nt
+  p = table(k,1);
+  expected = table(k,2);
+  q = 1 - p;
+  Mean = 0;
+  Std = 1;
+  computed = cdfnor ( "X" , Mean , Std , p , q );
+  assert_close ( computed , expected , precision );
+end
+
+//
+// Values from R-2.8.1
+// table = [x mu sigma PDF-P CDF-P CDF-Q]
+// Some tests do not pass with Scilab.
+//
+// See : http://bugzilla.scilab.org/show_bug.cgi?id=8032
+//
+// Prints the number of accurate digits.
+table = [
+//                    %inf 1.000000000000000000D+00 2.000000000000000000D+00 0.000000000000000000D+00 1.000000000000000000D+00 0.000000000000000000D+00
+//1.630146181031128094D+01 1.000000000000000000D+00 2.000000000000000000D+00 3.885547484725481156D-14 9.999999999999900080D-01 9.992007221626602924D-15
+1.372268177939483991D+01 1.000000000000000000D+00 2.000000000000000000D+00 3.255794261707647101D-10 9.999999998999999917D-01 1.000000082740384053D-10
+1.139867516458132002D+01 1.000000000000000000D+00 2.000000000000000000D+00 2.689766239123555948D-07 9.999999000000000526D-01 9.999999994736501079D-08
+9.529781587847679702D+00 1.000000000000000000D+00 2.000000000000000000D+00 2.239366490571142939D-05 9.999900000000000455D-01 9.999999999954530395D-06
+8.438032970911416797D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.979239833799469940D-04 9.999000000000000110D-01 9.999999999998900014D-05
+7.180464612335624608D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.683545038532002075D-03 9.989999999999999991D-01 1.000000000000003924D-03
+5.652695748081682403D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.332607110172901940D-02 9.899999999999999911D-01 1.000000000000000021D-02
+3.563103131089199849D+00 1.000000000000000000D+00 2.000000000000000000D+00 8.774916596624346421D-02 8.999999999999999112D-01 1.000000000000001027D-01
+2.683242467145829036D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.399809602039041034D-01 8.000000000000001554D-01 1.999999999999999001D-01
+1.506694206271600001D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.931712667484301871D-01 6.000000000000000888D-01 3.999999999999999112D-01
+//1.000000000000000000D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.994711402007164069D-01 5.000000000000000000D-01 5.000000000000000000D-01
+-6.832424671458280363D-01 1.000000000000000000D+00 2.000000000000000000D+00 1.399809602039041867D-01 2.000000000000000111D-01 8.000000000000000444D-01
+-1.563103131089200959D+00 1.000000000000000000D+00 2.000000000000000000D+00 8.774916596624339482D-02 1.000000000000000056D-01 9.000000000000000222D-01
+-3.652695748081681959D+00 1.000000000000000000D+00 2.000000000000000000D+00 1.332607110172904022D-02 1.000000000000001062D-02 9.899999999999999911D-01
+-7.529781587845650215D+00 1.000000000000000000D+00 2.000000000000000000D+00 2.239366490580828930D-05 9.999999999999960160D-06 9.999900000000000455D-01
+-1.172268180480810962D+01 1.000000000000000000D+00 2.000000000000000000D+00 3.255793998537787058D-10 1.000000000000008954D-10 9.999999998999999917D-01
+-1.752468017959679969D+01 1.000000000000000000D+00 2.000000000000000000D+00 4.683961267403037244D-20 1.000000000000072017D-20 1.000000000000000000D+00
+-2.886667506957698137D+01 1.000000000000000000D+00 2.000000000000000000D+00 7.499857137115675432D-50 9.999999999999744883D-51 1.000000000000000000D+00
+//                   -%inf 1.000000000000000000D+00 2.000000000000000000D+00 0.000000000000000000D+00 0.000000000000000000D+00 1.000000000000000000D+00
+];
 
-deff('y=f(t)','y=exp(-t^2/2)');
-if norm(yy-(1/2+ 1/sqrt(2*%pi)*intg(0,0.5,f)))> %eps then pause,end
+precision = 1.e-12;
+precinv = 1.e-8;
 
-[x1] = cdfnor("X",0.0*ones(x),1.0*ones(x),p,q);
-if norm(x-x1) > prec then pause,end 
+nt = size(table,"r");
+for k = 1 : nt
+  x = table(k,1);
+  mu = table(k,2);
+  std = table(k,3);
+  p = table(k,5);
+  q = table(k,6);
+  [ p1 , q1 ] = cdfnor("PQ",x,mu,std);
+  x1 = cdfnor("X",mu,std,p,q);
+  mu1 = cdfnor("Mean",std,p,q,x);
+  std1 = cdfnor("Std",p,q,x,mu);
+  if ( %t ) then
+    assert_close ( p1 , p , precision );
+    assert_close ( q1 , q , precision );
+    assert_close ( x1 , x , precinv );
+    assert_close ( mu1 , mu , precinv );
+    assert_close ( std1 , std , precinv );
+  end
+  if ( %f ) then
+    dP = assert_computedigits ( p1 , p );
+    dQ = assert_computedigits ( q1 , q );
+    dx = assert_computedigits ( x1 , x );
+    dmu = assert_computedigits ( mu1 , mu );
+    dstd = assert_computedigits ( std1 , std );
+    mprintf("Test #%3d/%3d: Digits p1= %.1f, q1=%.1f, X=%.1f, M=%.1f, S=%.1f\n",k,nt,dP,dQ,dx,dmu,dstd);
+  end
+end
 
-M     = 2*x;
-[p,q] = cdfnor("PQ",x,M,1.0*ones(x));
-M1    = cdfnor("Mean",1.0*ones(x),p,q,x);
-if norm(M-M1) > prec then pause,end
 
-// avoid M=x
-M      = 2*x+1;
-Std    = x+1;
-[p,q]  = cdfnor("PQ",x,M,Std);
-Std1   = cdfnor("Std",p,q,x,M);
-if norm(Std-Std1) > prec then pause,end
index bc1e366..7e6588b 100644 (file)
 // =============================================================================
-// Tests for cdfpoi() function
-//
-// Scilab Team
-// Copyright INRIA
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Michael Baudin
 //
+//  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-prec  = 1.e-5;
-xlam  = 0.34;
-deff('k=fact(n)','if n<=1;k=1;else k=n*fact(n-1);end');
-deff('s=S(k)','s=0;for j=0:k,s=s+ exp(-xlam)*(xlam)^j/fact(j);end');
-X     = 0:10;
-C     = feval(X,S);
-[P,Q] = cdfpoi("PQ",X,xlam*ones(X));
-if norm(P-C) > 100*%eps then bugmes();quit;end
-xlam   = 0.1*(1:11);
-X      = 0:10;
-[P,Q]  = cdfpoi("PQ",X,xlam);
-X1     = cdfpoi("S",xlam,P,Q);
-xlam1  = cdfpoi("Xlam",P,Q,X);
-if norm(X1-X) > prec then bugmes();quit;end
-if norm(xlam1-xlam) > prec then bugmes();quit;end
+// <-- JVM NOT MANDATORY -->
+// <-- ENGLISH IMPOSED -->
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+function d = assert_computedigits ( computed , expected )
+  nre = size(expected,"r")
+  nce = size(expected,"c")
+  // Update shape
+  expected = expected (:)
+  computed = computed (:)
+  //
+  dmin = 0
+  dmax = -log10(2^(-53))
+  //
+  d = zeros(expected)
+  //
+  n = size(expected,"*")
+  for i = 1 : n
+    if ( isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmax
+    elseif ( isnan(expected(i)) & ~isnan(computed(i)) ) then
+      d(i) = dmin
+    elseif ( ~isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmin
+      // From now, both expected and computed are non-nan
+    elseif ( expected(i) == 0 & computed(i) == 0 ) then
+      d(i) = dmax
+    elseif ( expected(i) == 0 & computed(i) <> 0 ) then
+      d(i) = dmin
+      // From now, expected(i) is non-zero
+    elseif ( expected(i) == computed(i) ) then
+      d(i) = dmax
+      // From now, expected and computed are different
+    elseif ( expected(i) == %inf & computed(i) <> %inf ) then
+      d(i) = dmin
+    elseif ( expected(i) == -%inf & computed(i) <> -%inf ) then
+      d(i) = dmin
+      // From now, neither of computed, nor expected is infinity
+    else
+      d(i) = max ( -log10 ( abs(computed(i)-expected(i)) / abs(expected(i)) ) , dmin )
+    end
+  end
+  //
+  // Reshape
+  d = matrix(d,nre,nce)
+endfunction
+//
+// Assessing the quality of the Normal distribution function
+// References
+//   Yalta, A. T. 2008. The accuracy of statistical distributions in Microsoft®Excel 2007. Comput. Stat. Data Anal. 52, 10 (Jun. 2008), 4579-4586. DOI= http://dx.doi.org/10.1016/j.csda.2008.03.005 
+//   Computation of Statistical Distributions (ELV), Leo Knüsel 
+// Table 4
+// Check Poisson distribution with parameters (lambda, k, Sigma)
+// If Sigma = 1, the cumulated distribution function is to be computed.
+//
+// table = [x lambda p precision]
+//
+table = [
+  1e+03 , 1e+03 , 0.508409 , 1.e-5
+  1e+05 , 1e+05 , 0.500841 , 1.e-5
+  1e+07 , 1e+07 , 0.500084 , 1.e-5
+  1e+09 , 1e+09 , 0.500008 , 1.e-5
+];
+nt = size(table,"r");
+for k = 1 : nt
+  Xk = table(k,1);
+  lambda = table(k,2);
+  expected = table(k,3);
+  precision = table(k,4);
+  [computed,Q]=cdfpoi("PQ",Xk,lambda);
+  assert_close ( computed , expected , precision );
+end
+//
+// Values from R-2.8.1
+// table = [x lambda PDF-P CDF-P CDF-Q]
+// Some tests do not pass with Scilab.
+//
+//
+// Prints the number of accurate digits.
+table = [
+1.000000000000000000D+03 1.000000000000000000D+03 1.261461134872150086D-02 5.084093671685060434D-01 4.915906328314940121D-01
+1.000000000000000000D+05 1.000000000000000000D+05 1.261565209705300949D-03 5.008410430993400775D-01 4.991589569006599225D-01
+1.000000000000000000D+07 1.000000000000000000D+07 1.261566250497027949D-04 5.000841044163260030D-01 4.999158955836739970D-01
+1.000000000000000000D+09 1.000000000000000000D+09 1.261566260904949930D-05 5.000084104417390485D-01 4.999915895582610070D-01
+//0.000000000000000000D+00 2.000000000000000000D+02 1.383896526736738008D-87 1.383896526736738008D-87 1.000000000000000000D+00
+5.000000000000000000D+01 2.000000000000000000D+02 5.123049239702292812D-37 6.815847235588002180D-37 1.000000000000000000D+00
+8.000000000000000000D+01 2.000000000000000000D+02 2.337628804093846817D-22 3.875088395176848769D-22 1.000000000000000000D+00
+1.030000000000000000D+02 2.000000000000000000D+02 1.417198708132908064D-14 2.891647405032629846D-14 9.999999999999711342D-01
+1.040000000000000000D+02 2.000000000000000000D+02 2.725382131024812095D-14 5.617029536057430897D-14 9.999999999999438227D-01
+1.330000000000000000D+02 2.000000000000000000D+02 1.013218285038530591D-07 2.943900186784097476D-07 9.999997056099813042D-01
+1.340000000000000000D+02 2.000000000000000000D+02 1.512266097072429104D-07 4.456166283856515992D-07 9.999995543833716249D-01
+2.000000000000000000D+02 2.000000000000000000D+02 2.819772768592081896D-02 5.187943096786845620D-01 4.812056903213154380D-01
+2.500000000000000000D+02 2.000000000000000000D+02 7.744905800132999957D-05 9.997153785997120456D-01 2.846214002883620114D-04
+2.800000000000000000D+02 2.000000000000000000D+02 1.602914502548084948D-08 9.999999615206405235D-01 3.847935945037056327D-08
+3.140000000000000000D+02 2.000000000000000000D+02 2.235682172206060878D-14 9.999999999999616973D-01 3.832269096110103145D-14
+3.150000000000000000D+02 2.000000000000000000D+02 1.419480744257838997D-14 9.999999999999759082D-01 2.412788351852324890D-14
+4.000000000000000000D+02 2.000000000000000000D+02 5.580687539454764804D-36 1.000000000000000000D+00 5.525962083726706326D-36
+6.000000000000000000D+02 2.000000000000000000D+02 4.53747297550031877D-115 1.000000000000000000D+00 2.26028138789988093D-115
+//9.000000000000000000D+02 2.000000000000000000D+02 1.73230170612912589D-286 1.000000000000000000D+00 4.94036667462315110D-287
+//1.000000000000000000D+03 2.000000000000000000D+02 0.000000000000000000D+00 1.000000000000000000D+00 0.000000000000000000D+00
+];
+precision = 1.e-12;
+precinv = 1.e-8;
+nt = size(table,"r");
+for k = 1 : nt
+  x = table(k,1);
+  lambda = table(k,2);
+  p = table(k,4);
+  q = table(k,5);
+  [ p1 , q1 ] = cdfpoi("PQ",x,lambda);
+  x1 = cdfpoi("S",lambda,p,q);
+  lambda1 = cdfpoi("Xlam",p,q,x);
+  if ( %t ) then
+    assert_close ( p1 , p , precision );
+    assert_close ( q1 , q , precision );
+    assert_close ( x1 , x , precinv );
+    assert_close ( lambda1 , lambda , precinv );
+  end
+  if ( %f ) then
+    dP = assert_computedigits ( p1 , p );
+    dQ = assert_computedigits ( q1 , q );
+    dx = assert_computedigits ( x1 , x );
+    dl = assert_computedigits ( lambda , lambda1 );
+    mprintf("Test #%3d/%3d: Digits p1= %.1f, q1=%.1f, X=%.1f, Lambda=%.1f\n",k,nt,dP,dQ,dx,dl);
+  end
+end
index 9445504..d0cd817 100644 (file)
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - INRIA - Michael Baudin
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
 
-// =============================================================================
-// Tests for cdfpoi() function
-// =============================================================================
 
-prec  = 1.e-5;
+// <-- JVM NOT MANDATORY -->
+// <-- ENGLISH IMPOSED -->
 
-xlam  = 0.34;
 
-deff('k=fact(n)','if n<=1;k=1;else k=n*fact(n-1);end');
-deff('s=S(k)','s=0;for j=0:k,s=s+ exp(-xlam)*(xlam)^j/fact(j);end');
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+function d = assert_computedigits ( computed , expected )
+  nre = size(expected,"r")
+  nce = size(expected,"c")
+  // Update shape
+  expected = expected (:)
+  computed = computed (:)
+  //
+  dmin = 0
+  dmax = -log10(2^(-53))
+  //
+  d = zeros(expected)
+  //
+  n = size(expected,"*")
+  for i = 1 : n
+    if ( isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmax
+    elseif ( isnan(expected(i)) & ~isnan(computed(i)) ) then
+      d(i) = dmin
+    elseif ( ~isnan(expected(i)) & isnan(computed(i)) ) then
+      d(i) = dmin
+      // From now, both expected and computed are non-nan
+    elseif ( expected(i) == 0 & computed(i) == 0 ) then
+      d(i) = dmax
+    elseif ( expected(i) == 0 & computed(i) <> 0 ) then
+      d(i) = dmin
+      // From now, expected(i) is non-zero
+    elseif ( expected(i) == computed(i) ) then
+      d(i) = dmax
+      // From now, expected and computed are different
+    elseif ( expected(i) == %inf & computed(i) <> %inf ) then
+      d(i) = dmin
+    elseif ( expected(i) == -%inf & computed(i) <> -%inf ) then
+      d(i) = dmin
+      // From now, neither of computed, nor expected is infinity
+    else
+      d(i) = max ( -log10 ( abs(computed(i)-expected(i)) / abs(expected(i)) ) , dmin )
+    end
+  end
+  //
+  // Reshape
+  d = matrix(d,nre,nce)
+endfunction
+
+
+//
+// Assessing the quality of the Normal distribution function
+// References
+//   Yalta, A. T. 2008. The accuracy of statistical distributions in Microsoft®Excel 2007. Comput. Stat. Data Anal. 52, 10 (Jun. 2008), 4579-4586. DOI= http://dx.doi.org/10.1016/j.csda.2008.03.005 
+//   Computation of Statistical Distributions (ELV), Leo Knüsel 
+// Table 4
+// Check Poisson distribution with parameters (lambda, k, Sigma)
+// If Sigma = 1, the cumulated distribution function is to be computed.
+//
+// table = [x lambda p precision]
+//
+
+table = [
+  1e+03 , 1e+03 , 0.508409 , 1.e-5
+  1e+05 , 1e+05 , 0.500841 , 1.e-5
+  1e+07 , 1e+07 , 0.500084 , 1.e-5
+  1e+09 , 1e+09 , 0.500008 , 1.e-5
+];
+
+nt = size(table,"r");
+for k = 1 : nt
+  Xk = table(k,1);
+  lambda = table(k,2);
+  expected = table(k,3);
+  precision = table(k,4);
+  [computed,Q]=cdfpoi("PQ",Xk,lambda);
+  assert_close ( computed , expected , precision );
+end
+
+//
+// Values from R-2.8.1
+// table = [x lambda PDF-P CDF-P CDF-Q]
+// Some tests do not pass with Scilab.
+//
+//
+// Prints the number of accurate digits.
 
-X     = 0:10;
-C     = feval(X,S);
-[P,Q] = cdfpoi("PQ",X,xlam*ones(X));
+table = [
+1.000000000000000000D+03 1.000000000000000000D+03 1.261461134872150086D-02 5.084093671685060434D-01 4.915906328314940121D-01
+1.000000000000000000D+05 1.000000000000000000D+05 1.261565209705300949D-03 5.008410430993400775D-01 4.991589569006599225D-01
+1.000000000000000000D+07 1.000000000000000000D+07 1.261566250497027949D-04 5.000841044163260030D-01 4.999158955836739970D-01
+1.000000000000000000D+09 1.000000000000000000D+09 1.261566260904949930D-05 5.000084104417390485D-01 4.999915895582610070D-01
+//0.000000000000000000D+00 2.000000000000000000D+02 1.383896526736738008D-87 1.383896526736738008D-87 1.000000000000000000D+00
+5.000000000000000000D+01 2.000000000000000000D+02 5.123049239702292812D-37 6.815847235588002180D-37 1.000000000000000000D+00
+8.000000000000000000D+01 2.000000000000000000D+02 2.337628804093846817D-22 3.875088395176848769D-22 1.000000000000000000D+00
+1.030000000000000000D+02 2.000000000000000000D+02 1.417198708132908064D-14 2.891647405032629846D-14 9.999999999999711342D-01
+1.040000000000000000D+02 2.000000000000000000D+02 2.725382131024812095D-14 5.617029536057430897D-14 9.999999999999438227D-01
+1.330000000000000000D+02 2.000000000000000000D+02 1.013218285038530591D-07 2.943900186784097476D-07 9.999997056099813042D-01
+1.340000000000000000D+02 2.000000000000000000D+02 1.512266097072429104D-07 4.456166283856515992D-07 9.999995543833716249D-01
+2.000000000000000000D+02 2.000000000000000000D+02 2.819772768592081896D-02 5.187943096786845620D-01 4.812056903213154380D-01
+2.500000000000000000D+02 2.000000000000000000D+02 7.744905800132999957D-05 9.997153785997120456D-01 2.846214002883620114D-04
+2.800000000000000000D+02 2.000000000000000000D+02 1.602914502548084948D-08 9.999999615206405235D-01 3.847935945037056327D-08
+3.140000000000000000D+02 2.000000000000000000D+02 2.235682172206060878D-14 9.999999999999616973D-01 3.832269096110103145D-14
+3.150000000000000000D+02 2.000000000000000000D+02 1.419480744257838997D-14 9.999999999999759082D-01 2.412788351852324890D-14
+4.000000000000000000D+02 2.000000000000000000D+02 5.580687539454764804D-36 1.000000000000000000D+00 5.525962083726706326D-36
+6.000000000000000000D+02 2.000000000000000000D+02 4.53747297550031877D-115 1.000000000000000000D+00 2.26028138789988093D-115
+//9.000000000000000000D+02 2.000000000000000000D+02 1.73230170612912589D-286 1.000000000000000000D+00 4.94036667462315110D-287
+//1.000000000000000000D+03 2.000000000000000000D+02 0.000000000000000000D+00 1.000000000000000000D+00 0.000000000000000000D+00
+];
 
-if norm(P-C) > 100*%eps then pause,end
+precision = 1.e-12;
+precinv = 1.e-8;
+nt = size(table,"r");
+for k = 1 : nt
+  x = table(k,1);
+  lambda = table(k,2);
+  p = table(k,4);
+  q = table(k,5);
+  [ p1 , q1 ] = cdfpoi("PQ",x,lambda);
+  x1 = cdfpoi("S",lambda,p,q);
+  lambda1 = cdfpoi("Xlam",p,q,x);
+  if ( %t ) then
+    assert_close ( p1 , p , precision );
+    assert_close ( q1 , q , precision );
+    assert_close ( x1 , x , precinv );
+    assert_close ( lambda1 , lambda , precinv );
+  end
+  if ( %f ) then
+    dP = assert_computedigits ( p1 , p );
+    dQ = assert_computedigits ( q1 , q );
+    dx = assert_computedigits ( x1 , x );
+    dl = assert_computedigits ( lambda , lambda1 );
+    mprintf("Test #%3d/%3d: Digits p1= %.1f, q1=%.1f, X=%.1f, Lambda=%.1f\n",k,nt,dP,dQ,dx,dl);
+  end
+end
 
-xlam   = 0.1*(1:11);
-X      = 0:10;
-[P,Q]  = cdfpoi("PQ",X,xlam);
-X1     = cdfpoi("S",xlam,P,Q);
-xlam1  = cdfpoi("Xlam",P,Q,X);
 
-if norm(X1-X) > prec then pause,end
-if norm(xlam1-xlam) > prec then pause,end