statistics: Fixed bug #7768: for cdfgam, the Scale parameter was, in fact, the Rate
[scilab.git] / scilab / modules / statistics / tests / unit_tests / cdfgam.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2010 - DIGITEO - Michael Baudin
4 //
5 //  This file is distributed under the same license as the Scilab package.
6 // =============================================================================
7
8
9 // <-- JVM NOT MANDATORY -->
10 // <-- ENGLISH IMPOSED -->
11
12 //
13 // assert_close --
14 //   Returns 1 if the two real matrices computed and expected are close,
15 //   i.e. if the relative distance between computed and expected is lesser than epsilon.
16 // Arguments
17 //   computed, expected : the two matrices to compare
18 //   epsilon : a small number
19 //
20 function flag = assert_close ( computed, expected, epsilon )
21   if expected==0.0 then
22     shift = norm(computed-expected);
23   else
24     shift = norm(computed-expected)/norm(expected);
25   end
26   if shift < epsilon then
27     flag = 1;
28   else
29     flag = 0;
30   end
31   if flag <> 1 then pause,end
32 endfunction
33 //
34 // assert_equal --
35 //   Returns 1 if the two real matrices computed and expected are equal.
36 // Arguments
37 //   computed, expected : the two matrices to compare
38 //   epsilon : a small number
39 //
40 function flag = assert_equal ( computed , expected )
41   if computed==expected then
42     flag = 1;
43   else
44     flag = 0;
45   end
46   if flag <> 1 then pause,end
47 endfunction
48
49 function d = assert_computedigits ( computed , expected )
50   nre = size(expected,"r")
51   nce = size(expected,"c")
52   // Update shape
53   expected = expected (:)
54   computed = computed (:)
55   //
56   dmin = 0
57   dmax = -log10(2^(-53))
58   //
59   d = zeros(expected)
60   //
61   n = size(expected,"*")
62   for i = 1 : n
63     if ( isnan(expected(i)) & isnan(computed(i)) ) then
64       d(i) = dmax
65     elseif ( isnan(expected(i)) & ~isnan(computed(i)) ) then
66       d(i) = dmin
67     elseif ( ~isnan(expected(i)) & isnan(computed(i)) ) then
68       d(i) = dmin
69       // From now, both expected and computed are non-nan
70     elseif ( expected(i) == 0 & computed(i) == 0 ) then
71       d(i) = dmax
72     elseif ( expected(i) == 0 & computed(i) <> 0 ) then
73       d(i) = dmin
74       // From now, expected(i) is non-zero
75     elseif ( expected(i) == computed(i) ) then
76       d(i) = dmax
77       // From now, expected and computed are different
78     elseif ( expected(i) == %inf & computed(i) <> %inf ) then
79       d(i) = dmin
80     elseif ( expected(i) == -%inf & computed(i) <> -%inf ) then
81       d(i) = dmin
82       // From now, neither of computed, nor expected is infinity
83     else
84       d(i) = max ( -log10 ( abs(computed(i)-expected(i)) / abs(expected(i)) ) , dmin )
85     end
86   end
87   //
88   // Reshape
89   d = matrix(d,nre,nce)
90 endfunction
91
92 //
93 // Assessing the quality of the Normal distribution function
94 // References
95 //   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 
96 //   Computation of Statistical Distributions (ELV), Leo Kn├╝sel 
97 // Table 5
98 // Check Gamma distribution with parameters (x, alpha, beta = 1, Sigma = 1)
99 //
100
101
102 // Table of inputs from Yalta, 2008
103 // [x shape scale P ]
104 table = [
105  0.1 , 0.1 , 1 , 0.827552 
106  0.2 , 0.1 , 1 , 0.879420 
107  0.2 , 0.2 , 1 , 0.764435 
108  0.3 , 0.2 , 1 , 0.816527 
109  0.3 , 0.3 , 1 , 0.726957 
110  0.4 , 0.3 , 1 , 0.776381 
111  0.4 , 0.4 , 1 , 0.701441 
112  0.5 , 0.4 , 1 , 0.748019 
113  0.5 , 0.5 , 1 , 0.682689 
114  0.6 , 0.5 , 1 , 0.726678 
115 ];
116
117 precision = 1.e-5;
118 ntests = size(table,"r");
119 for i = 1 : ntests
120   x = table(i,1);
121   shape = table(i,2);
122   scale = table(i,3);
123   expected = table(i,4);
124   // Caution: this is the rate !
125   rate = 1/scale;
126   [computed,Q]=cdfgam("PQ",x,shape,rate);
127   assert_close ( computed , expected , precision );
128   assert_close ( Q , 1 - expected , precision );
129 end
130
131 // Table of inputs computed from R-2.8.1
132 // [x shape scale PDF-P CDF-P CDF-Q]
133 table = [
134 1.000000000000000056D-01 1.000000000000000056D-01 1.000000000000000000D+00 7.554920138253073958D-01 8.275517595858505882D-01 1.724482404141494951D-01
135 2.000000000000000111D-01 1.000000000000000056D-01 1.000000000000000000D+00 3.663307993056703071D-01 8.794196267900568076D-01 1.205803732099432063D-01
136 2.000000000000000111D-01 2.000000000000000111D-01 1.000000000000000000D+00 6.462857778271943188D-01 7.644345975029189777D-01 2.355654024970809945D-01
137 2.999999999999999889D-01 2.000000000000000111D-01 1.000000000000000000D+00 4.227875047602157044D-01 8.165267943336527168D-01 1.834732056663473110D-01
138 2.999999999999999889D-01 2.999999999999999889D-01 1.000000000000000000D+00 5.752117576599179438D-01 7.269573437103662439D-01 2.730426562896338116D-01
139 4.000000000000000222D-01 2.999999999999999889D-01 1.000000000000000000D+00 4.255407854753925911D-01 7.763805810166358734D-01 2.236194189833642099D-01
140 4.000000000000000222D-01 4.000000000000000222D-01 1.000000000000000000D+00 5.236648604477927016D-01 7.014412706419403953D-01 2.985587293580597157D-01
141 5.000000000000000000D-01 4.000000000000000222D-01 1.000000000000000000D+00 4.144555659263016167D-01 7.480185547260104206D-01 2.519814452739895239D-01
142 5.000000000000000000D-01 5.000000000000000000D-01 1.000000000000000000D+00 4.839414490382866751D-01 6.826894921370858516D-01 3.173105078629140929D-01
143 5.999999999999999778D-01 5.000000000000000000D-01 1.000000000000000000D+00 3.997355278034666060D-01 7.266783217077018575D-01 2.733216782922981980D-01
144 5.000000000000000000D-01 5.000000000000000000D-01 2.000000000000000000D+00 4.393912894677223790D-01 5.204998778130465187D-01 4.795001221869534258D-01
145 5.000000000000000000D-01 5.000000000000000000D-01 3.000000000000000000D+00 3.899393114454822729D-01 4.362971383492270094D-01 5.637028616507729906D-01
146 5.000000000000000000D-01 5.000000000000000000D-01 4.000000000000000000D+00 3.520653267642995243D-01 3.829249225480261809D-01 6.170750774519737636D-01
147 1.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 2.075537487102973866D-01 8.427007929497148941D-01 1.572992070502851891D-01
148 2.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 5.399096651318804896D-02 9.544997361036415828D-01 4.550026389635838248D-02
149 4.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 5.166746338523012412D-03 9.953222650189527121D-01 4.677734981047261889D-03
150 1.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 8.099910956089122777D-06 9.999922557835689840D-01 7.744216431044085842D-06
151 2.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 2.600281868827196957D-10 9.999999997460371493D-01 2.539628589470869077D-10
152 4.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 3.789795640412981196D-19 1.000000000000000000D+00 3.744097384202895045D-19
153 //1.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 2.098828115677222045D-45 1.000000000000000000D+00 2.088487583762558879D-45
154 3.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 1.67694904029982009D-132 1.000000000000000000D+00 1.67416798469182012D-132
155 //5.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 1.79762504374667411D-219 1.000000000000000000D+00 1.79583278480075297D-219
156 //1.000000000000000000D+03 5.000000000000000000D-01 1.000000000000000000D+00 0.000000000000000000D+00 1.000000000000000000D+00 0.000000000000000000D+00
157 1.000000000000000021D-02 5.000000000000000000D-01 1.000000000000000000D+00 5.585758033944684620D+00 1.124629160182848975D-01 8.875370839817151580D-01
158 1.000000000000000048D-04 5.000000000000000000D-01 1.000000000000000000D+00 5.641331674102550409D+01 1.128341555584961957D-02 9.887165844441503371D-01
159 1.000000000000000021D-08 5.000000000000000000D-01 1.000000000000000000D+00 5.641895779058606422D+03 1.128379163334249004D-04 9.998871620836665697D-01
160 9.999999999999999452D-21 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477570534D+09 1.128379167095512970D-10 9.999999998871620388D-01
161 9.999999999999999293D-41 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477568717D+19 1.128379167095512972D-20 1.000000000000000000D+00
162 1.00000000000000002D-100 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477541988D+49 1.128379167095513082D-50 1.000000000000000000D+00
163 9.99999999999999982D-201 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477511468D+99 1.12837916709551300D-100 1.000000000000000000D+00
164 //1.00000000000000002D-300 5.000000000000000000D-01 1.000000000000000000D+00 5.64189583547731891D+149 1.12837916709551298D-150 1.000000000000000000D+00
165 //0.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00                     %inf 0.000000000000000000D+00 1.000000000000000000D+00
166 ];
167
168 // For the inversion of Shape, require only 8 digits, as 
169 // a consequence of bug #7569: http://bugzilla.scilab.org/show_bug.cgi?id=7569
170 //
171 // Some tests do not pass:
172 // http://bugzilla.scilab.org/show_bug.cgi?id=8031
173 // http://bugzilla.scilab.org/show_bug.cgi?id=8030
174 //
175 // Prints the number of accurate digits.
176
177 precision = 1.e-12;
178 precinverse = 1.e-8;
179
180 ntests = size(table,"r");
181 for i = 1 : ntests
182   x = table(i,1);
183   shape = table(i,2);
184   scale = table(i,3);
185   p = table(i,5);
186   q = table(i,6);
187   // Caution: this is the rate !
188   rate = 1/scale;
189   [p1,q1] = cdfgam("PQ",x,shape,rate);
190   x1      = cdfgam("X",shape,rate,p,q);
191   shape1  = cdfgam("Shape",rate,p,q,x);
192   rate1   = cdfgam("Rate",p,q,x,shape);
193   if ( %t ) then
194     assert_close ( p1 , p , precision );
195     assert_close ( q1 , q , precision );
196     assert_close ( x1 , x , precision );
197     assert_close ( shape1 , shape , precinverse );
198     assert_close ( rate1 , rate , precinverse );
199   end
200   if ( %f ) then
201     dp = assert_computedigits ( p1 , p );
202     dq = assert_computedigits ( q1 , q );
203     dx = assert_computedigits ( x1 , x );
204     ds = assert_computedigits ( shape1 , shape );
205     dr = assert_computedigits ( rate1 , rate );
206     mprintf("Test #%3d/%3d: Digits p1= %.1f, q1=%.1f, X= %.1f, S= %.1f, R= %.1f\n",i,ntests,dp,dq,dx,ds,dr);
207   end
208 end
209