|
0
|
1 #include "stats.hpp"
|
|
|
2
|
|
|
3 /** The digamma function in long double precision.
|
|
|
4 * @param x the real value of the argument
|
|
|
5 * @return the value of the digamma (psi) function at that point
|
|
|
6 * @author Richard J. Mathar
|
|
|
7 * @since 2005-11-24
|
|
|
8 */
|
|
|
9 long double digammal(long double x)
|
|
|
10 {
|
|
|
11 /* force into the interval 1..3 */
|
|
|
12 if( x < 0.0L )
|
|
|
13 return digammal(1.0L-x)+M_PIl/tanl(M_PIl*(1.0L-x)) ; /* reflection formula */
|
|
|
14 else if( x < 1.0L )
|
|
|
15 return digammal(1.0L+x)-1.0L/x ;
|
|
|
16 else if ( x == 1.0L)
|
|
|
17 return -M_GAMMAl ;
|
|
|
18 else if ( x == 2.0L)
|
|
|
19 return 1.0L-M_GAMMAl ;
|
|
|
20 else if ( x == 3.0L)
|
|
|
21 return 1.5L-M_GAMMAl ;
|
|
|
22 else if ( x > 3.0L)
|
|
|
23 /* duplication formula */
|
|
|
24 return 0.5L*(digammal(x/2.0L)+digammal((x+1.0L)/2.0L))+M_LN2l ;
|
|
|
25 else
|
|
|
26 {
|
|
|
27 /* Just for your information, the following lines contain
|
|
|
28 * the Maple source code to re-generate the table that is
|
|
|
29 * eventually becoming the Kncoe[] array below
|
|
|
30 * interface(prettyprint=0) :
|
|
|
31 * Digits := 63 :
|
|
|
32 * r := 0 :
|
|
|
33 *
|
|
|
34 * for l from 1 to 60 do
|
|
|
35 * d := binomial(-1/2,l) :
|
|
|
36 * r := r+d*(-1)^l*(Zeta(2*l+1) -1) ;
|
|
|
37 * evalf(r) ;
|
|
|
38 * print(%,evalf(1+Psi(1)-r)) ;
|
|
|
39 *o d :
|
|
|
40 *
|
|
|
41 * for N from 1 to 28 do
|
|
|
42 * r := 0 :
|
|
|
43 * n := N-1 :
|
|
|
44 *
|
|
|
45 * for l from iquo(n+3,2) to 70 do
|
|
|
46 * d := 0 :
|
|
|
47 * for s from 0 to n+1 do
|
|
|
48 * d := d+(-1)^s*binomial(n+1,s)*binomial((s-1)/2,l) :
|
|
|
49 * od :
|
|
|
50 * if 2*l-n > 1 then
|
|
|
51 * r := r+d*(-1)^l*(Zeta(2*l-n) -1) :
|
|
|
52 * fi :
|
|
|
53 * od :
|
|
|
54 * print(evalf((-1)^n*2*r)) ;
|
|
|
55 *od :
|
|
|
56 *quit :
|
|
|
57 */
|
|
|
58 static long double Kncoe[] = { .30459198558715155634315638246624251L,
|
|
|
59 .72037977439182833573548891941219706L, -.12454959243861367729528855995001087L,
|
|
|
60 .27769457331927827002810119567456810e-1L, -.67762371439822456447373550186163070e-2L,
|
|
|
61 .17238755142247705209823876688592170e-2L, -.44817699064252933515310345718960928e-3L,
|
|
|
62 .11793660000155572716272710617753373e-3L, -.31253894280980134452125172274246963e-4L,
|
|
|
63 .83173997012173283398932708991137488e-5L, -.22191427643780045431149221890172210e-5L,
|
|
|
64 .59302266729329346291029599913617915e-6L, -.15863051191470655433559920279603632e-6L,
|
|
|
65 .42459203983193603241777510648681429e-7L, -.11369129616951114238848106591780146e-7L,
|
|
|
66 .304502217295931698401459168423403510e-8L, -.81568455080753152802915013641723686e-9L,
|
|
|
67 .21852324749975455125936715817306383e-9L, -.58546491441689515680751900276454407e-10L,
|
|
|
68 .15686348450871204869813586459513648e-10L, -.42029496273143231373796179302482033e-11L,
|
|
|
69 .11261435719264907097227520956710754e-11L, -.30174353636860279765375177200637590e-12L,
|
|
|
70 .80850955256389526647406571868193768e-13L, -.21663779809421233144009565199997351e-13L,
|
|
|
71 .58047634271339391495076374966835526e-14L, -.15553767189204733561108869588173845e-14L,
|
|
|
72 .41676108598040807753707828039353330e-15L, -.11167065064221317094734023242188463e-15L } ;
|
|
|
73
|
|
|
74 register long double Tn_1 = 1.0L ; /* T_{n-1}(x), started at n=1 */
|
|
|
75 register long double Tn = x-2.0L ; /* T_{n}(x) , started at n=1 */
|
|
|
76 register long double resul = Kncoe[0] + Kncoe[1]*Tn ;
|
|
|
77
|
|
|
78 x -= 2.0L ;
|
|
|
79 int n ;
|
|
|
80
|
|
|
81 for( n = 2 ; n < int( sizeof(Kncoe)/sizeof(long double) ) ; n++)
|
|
|
82 {
|
|
|
83 const long double Tn1 = 2.0L * x * Tn - Tn_1 ; /* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */
|
|
|
84 resul += Kncoe[n]*Tn1 ;
|
|
|
85 Tn_1 = Tn ;
|
|
|
86 Tn = Tn1 ;
|
|
|
87 }
|
|
|
88 return resul ;
|
|
|
89 }
|
|
|
90 }
|
|
|
91
|
|
|
92
|
|
|
93
|
|
|
94 double trigamma ( double x, int *ifault )
|
|
|
95
|
|
|
96 //****************************************************************************
|
|
|
97 // purpose:
|
|
|
98 //
|
|
|
99 // trigamma calculates trigamma(x) = d**2 log(gamma(x)) / dx**2
|
|
|
100 //
|
|
|
101 // licensing:
|
|
|
102 //
|
|
|
103 // this code is distributed under the gnu lgpl license.
|
|
|
104 //
|
|
|
105 // modified:
|
|
|
106 //
|
|
|
107 // 19 january 2008
|
|
|
108 //
|
|
|
109 // author:
|
|
|
110 //
|
|
|
111 // original fortran77 version by be schneider.
|
|
|
112 // c++ version by john burkardt.
|
|
|
113 //
|
|
|
114 // reference:
|
|
|
115 //
|
|
|
116 // be schneider,
|
|
|
117 // algorithm as 121:
|
|
|
118 // trigamma function,
|
|
|
119 // applied statistics,
|
|
|
120 // volume 27, number 1, pages 97-99, 1978.
|
|
|
121 //
|
|
|
122 // parameters:
|
|
|
123 //
|
|
|
124 // input, double x, the argument of the trigamma function.
|
|
|
125 // 0 < x.
|
|
|
126 //
|
|
|
127 // output, int *ifault, error flag.
|
|
|
128 // 0, no error.
|
|
|
129 // 1, x <= 0.
|
|
|
130 //
|
|
|
131 // output, double trigamma, the value of the trigamma function at x.
|
|
|
132 //
|
|
|
133 {
|
|
|
134 double a = 0.0001;
|
|
|
135 double b = 5.0;
|
|
|
136 double b2 = 0.1666666667;
|
|
|
137 double b4 = -0.03333333333;
|
|
|
138 double b6 = 0.02380952381;
|
|
|
139 double b8 = -0.03333333333;
|
|
|
140 double value;
|
|
|
141 double y;
|
|
|
142 double z;
|
|
|
143 //
|
|
|
144 // check the input.
|
|
|
145 //
|
|
|
146 if ( x <= 0.0 )
|
|
|
147 {
|
|
|
148 *ifault = 1;
|
|
|
149 value = 0.0;
|
|
|
150 return value;
|
|
|
151 }
|
|
|
152
|
|
|
153 *ifault = 0;
|
|
|
154 z = x;
|
|
|
155 //
|
|
|
156 // use small value approximation if x <= a.
|
|
|
157 //
|
|
|
158 if ( x <= a )
|
|
|
159 {
|
|
|
160 value = 1.0 / x / x;
|
|
|
161 return value;
|
|
|
162 }
|
|
|
163 //
|
|
|
164 // increase argument to ( x + i ) >= b.
|
|
|
165 //
|
|
|
166 value = 0.0;
|
|
|
167
|
|
|
168 while ( z < b )
|
|
|
169 {
|
|
|
170 value = value + 1.0 / z / z;
|
|
|
171 z = z + 1.0;
|
|
|
172 }
|
|
|
173 //
|
|
|
174 // apply asymptotic formula if argument is b or greater.
|
|
|
175 //
|
|
|
176 y = 1.0 / z / z;
|
|
|
177
|
|
|
178 value = value + 0.5 *
|
|
|
179 y + ( 1.0 + y * ( b2+ y * ( b4 + y * ( b6+ y * b8 )))) / z;
|
|
|
180
|
|
|
181 return value;
|
|
|
182 }
|
|
|
183
|
|
|
184
|
|
|
185 double LogGammaDensity( double x, double k, double theta )
|
|
|
186 {
|
|
|
187 return -k * log( theta ) + ( k - 1 ) * log( x ) - x / theta - lgamma( k ) ;
|
|
|
188 }
|
|
|
189
|
|
|
190 double MixtureGammaAssignment( double x, double pi, double* k, double *theta )
|
|
|
191 {
|
|
|
192 if ( pi == 1 )
|
|
|
193 return 1 ;
|
|
|
194 else if ( pi == 0 )
|
|
|
195 return 0 ;
|
|
|
196
|
|
|
197 double lf0 = LogGammaDensity( x, k[0], theta[0] ) ;
|
|
|
198 double lf1 = LogGammaDensity( x, k[1], theta[1] ) ;
|
|
|
199 return (double)1.0 / ( 1.0 + exp( lf1 + log( 1 - pi ) - lf0 - log( pi ) ) ) ;
|
|
|
200 }
|
|
|
201
|
|
|
202 //****************************************************************************80
|
|
|
203
|
|
|
204 double alnorm ( double x, bool upper )
|
|
|
205
|
|
|
206 //****************************************************************************80
|
|
|
207 //
|
|
|
208 // Purpose:
|
|
|
209 //
|
|
|
210 // ALNORM computes the cumulative density of the standard normal distribution.
|
|
|
211 //
|
|
|
212 // Licensing:
|
|
|
213 //
|
|
|
214 // This code is distributed under the GNU LGPL license.
|
|
|
215 //
|
|
|
216 // Modified:
|
|
|
217 //
|
|
|
218 // 17 January 2008
|
|
|
219 //
|
|
|
220 // Author:
|
|
|
221 //
|
|
|
222 // Original FORTRAN77 version by David Hill.
|
|
|
223 // C++ version by John Burkardt.
|
|
|
224 //
|
|
|
225 // Reference:
|
|
|
226 //
|
|
|
227 // David Hill,
|
|
|
228 // Algorithm AS 66:
|
|
|
229 // The Normal Integral,
|
|
|
230 // Applied Statistics,
|
|
|
231 // Volume 22, Number 3, 1973, pages 424-427.
|
|
|
232 //
|
|
|
233 // Parameters:
|
|
|
234 //
|
|
|
235 // Input, double X, is one endpoint of the semi-infinite interval
|
|
|
236 // over which the integration takes place.
|
|
|
237 //
|
|
|
238 // Input, bool UPPER, determines whether the upper or lower
|
|
|
239 // interval is to be integrated:
|
|
|
240 // .TRUE. => integrate from X to + Infinity;
|
|
|
241 // .FALSE. => integrate from - Infinity to X.
|
|
|
242 //
|
|
|
243 // Output, double ALNORM, the integral of the standard normal
|
|
|
244 // distribution over the desired interval.
|
|
|
245 //
|
|
|
246 {
|
|
|
247 double a1 = 5.75885480458;
|
|
|
248 double a2 = 2.62433121679;
|
|
|
249 double a3 = 5.92885724438;
|
|
|
250 double b1 = -29.8213557807;
|
|
|
251 double b2 = 48.6959930692;
|
|
|
252 double c1 = -0.000000038052;
|
|
|
253 double c2 = 0.000398064794;
|
|
|
254 double c3 = -0.151679116635;
|
|
|
255 double c4 = 4.8385912808;
|
|
|
256 double c5 = 0.742380924027;
|
|
|
257 double c6 = 3.99019417011;
|
|
|
258 double con = 1.28;
|
|
|
259 double d1 = 1.00000615302;
|
|
|
260 double d2 = 1.98615381364;
|
|
|
261 double d3 = 5.29330324926;
|
|
|
262 double d4 = -15.1508972451;
|
|
|
263 double d5 = 30.789933034;
|
|
|
264 double ltone = 7.0;
|
|
|
265 double p = 0.398942280444;
|
|
|
266 double q = 0.39990348504;
|
|
|
267 double r = 0.398942280385;
|
|
|
268 bool up;
|
|
|
269 double utzero = 18.66;
|
|
|
270 double value;
|
|
|
271 double y;
|
|
|
272 double z;
|
|
|
273
|
|
|
274 up = upper;
|
|
|
275 z = x;
|
|
|
276
|
|
|
277 if ( z < 0.0 )
|
|
|
278 {
|
|
|
279 up = !up;
|
|
|
280 z = - z;
|
|
|
281 }
|
|
|
282
|
|
|
283 if ( ltone < z && ( ( !up ) || utzero < z ) )
|
|
|
284 {
|
|
|
285 if ( up )
|
|
|
286 {
|
|
|
287 value = 0.0;
|
|
|
288 }
|
|
|
289 else
|
|
|
290 {
|
|
|
291 value = 1.0;
|
|
|
292 }
|
|
|
293 return value;
|
|
|
294 }
|
|
|
295
|
|
|
296 y = 0.5 * z * z;
|
|
|
297
|
|
|
298 if ( z <= con )
|
|
|
299 {
|
|
|
300 value = 0.5 - z * ( p - q * y
|
|
|
301 / ( y + a1 + b1
|
|
|
302 / ( y + a2 + b2
|
|
|
303 / ( y + a3 ))));
|
|
|
304 }
|
|
|
305 else
|
|
|
306 {
|
|
|
307 value = r * exp ( - y )
|
|
|
308 / ( z + c1 + d1
|
|
|
309 / ( z + c2 + d2
|
|
|
310 / ( z + c3 + d3
|
|
|
311 / ( z + c4 + d4
|
|
|
312 / ( z + c5 + d5
|
|
|
313 / ( z + c6 ))))));
|
|
|
314 }
|
|
|
315
|
|
|
316 if ( !up )
|
|
|
317 {
|
|
|
318 value = 1.0 - value;
|
|
|
319 }
|
|
|
320
|
|
|
321 return value;
|
|
|
322 }
|
|
|
323
|
|
|
324 //****************************************************************************80
|
|
|
325
|
|
|
326 double gammad ( double x, double p, int *ifault )
|
|
|
327
|
|
|
328 //****************************************************************************80
|
|
|
329 //
|
|
|
330 // Purpose:
|
|
|
331 //
|
|
|
332 // GAMMAD computes the Incomplete Gamma Integral
|
|
|
333 //
|
|
|
334 // Licensing:
|
|
|
335 //
|
|
|
336 // This code is distributed under the GNU LGPL license.
|
|
|
337 //
|
|
|
338 // Modified:
|
|
|
339 //
|
|
|
340 // 20 January 2008
|
|
|
341 //
|
|
|
342 // Author:
|
|
|
343 //
|
|
|
344 // Original FORTRAN77 version by B Shea.
|
|
|
345 // C++ version by John Burkardt.
|
|
|
346 //
|
|
|
347 // Reference:
|
|
|
348 //
|
|
|
349 // B Shea,
|
|
|
350 // Algorithm AS 239:
|
|
|
351 // Chi-squared and Incomplete Gamma Integral,
|
|
|
352 // Applied Statistics,
|
|
|
353 // Volume 37, Number 3, 1988, pages 466-473.
|
|
|
354 //
|
|
|
355 // Parameters:
|
|
|
356 //
|
|
|
357 // Input, double X, P, the parameters of the incomplete
|
|
|
358 // gamma ratio. 0 <= X, and 0 < P.
|
|
|
359 //
|
|
|
360 // Output, int IFAULT, error flag.
|
|
|
361 // 0, no error.
|
|
|
362 // 1, X < 0 or P <= 0.
|
|
|
363 //
|
|
|
364 // Output, double GAMMAD, the value of the incomplete
|
|
|
365 // Gamma integral.
|
|
|
366 //
|
|
|
367 {
|
|
|
368 double a;
|
|
|
369 double an;
|
|
|
370 double arg;
|
|
|
371 double b;
|
|
|
372 double c;
|
|
|
373 double elimit = - 88.0;
|
|
|
374 double oflo = 1.0E+37;
|
|
|
375 double plimit = 1000.0;
|
|
|
376 double pn1;
|
|
|
377 double pn2;
|
|
|
378 double pn3;
|
|
|
379 double pn4;
|
|
|
380 double pn5;
|
|
|
381 double pn6;
|
|
|
382 double rn;
|
|
|
383 double tol = 1.0E-14;
|
|
|
384 bool upper;
|
|
|
385 double value;
|
|
|
386 double xbig = 1.0E+08;
|
|
|
387
|
|
|
388 value = 0.0;
|
|
|
389 //
|
|
|
390 // Check the input.
|
|
|
391 //
|
|
|
392 if ( x < 0.0 )
|
|
|
393 {
|
|
|
394 *ifault = 1;
|
|
|
395 return value;
|
|
|
396 }
|
|
|
397
|
|
|
398 if ( p <= 0.0 )
|
|
|
399 {
|
|
|
400 *ifault = 1;
|
|
|
401 return value;
|
|
|
402 }
|
|
|
403
|
|
|
404 *ifault = 0;
|
|
|
405
|
|
|
406 if ( x == 0.0 )
|
|
|
407 {
|
|
|
408 value = 0.0;
|
|
|
409 return value;
|
|
|
410 }
|
|
|
411 //
|
|
|
412 // If P is large, use a normal approximation.
|
|
|
413 //
|
|
|
414 if ( plimit < p )
|
|
|
415 {
|
|
|
416 pn1 = 3.0 * sqrt ( p ) * ( pow ( x / p, 1.0 / 3.0 )
|
|
|
417 + 1.0 / ( 9.0 * p ) - 1.0 );
|
|
|
418
|
|
|
419 upper = false;
|
|
|
420 value = alnorm ( pn1, upper );
|
|
|
421 return value;
|
|
|
422 }
|
|
|
423 //
|
|
|
424 // If X is large set value = 1.
|
|
|
425 //
|
|
|
426 if ( xbig < x )
|
|
|
427 {
|
|
|
428 value = 1.0;
|
|
|
429 return value;
|
|
|
430 }
|
|
|
431 //
|
|
|
432 // Use Pearson's series expansion.
|
|
|
433 // (Note that P is not large enough to force overflow in ALOGAM).
|
|
|
434 // No need to test IFAULT on exit since P > 0.
|
|
|
435 //
|
|
|
436 if ( x <= 1.0 || x < p )
|
|
|
437 {
|
|
|
438 arg = p * log ( x ) - x - lgamma ( p + 1.0 );
|
|
|
439 c = 1.0;
|
|
|
440 value = 1.0;
|
|
|
441 a = p;
|
|
|
442
|
|
|
443 for ( ; ; )
|
|
|
444 {
|
|
|
445 a = a + 1.0;
|
|
|
446 c = c * x / a;
|
|
|
447 value = value + c;
|
|
|
448
|
|
|
449 if ( c <= tol )
|
|
|
450 {
|
|
|
451 break;
|
|
|
452 }
|
|
|
453 }
|
|
|
454
|
|
|
455 arg = arg + log ( value );
|
|
|
456
|
|
|
457 if ( elimit <= arg )
|
|
|
458 {
|
|
|
459 value = exp ( arg );
|
|
|
460 }
|
|
|
461 else
|
|
|
462 {
|
|
|
463 value = 0.0;
|
|
|
464 }
|
|
|
465 }
|
|
|
466 //
|
|
|
467 // Use a continued fraction expansion.
|
|
|
468 //
|
|
|
469 else
|
|
|
470 {
|
|
|
471 arg = p * log ( x ) - x - lgamma ( p );
|
|
|
472 a = 1.0 - p;
|
|
|
473 b = a + x + 1.0;
|
|
|
474 c = 0.0;
|
|
|
475 pn1 = 1.0;
|
|
|
476 pn2 = x;
|
|
|
477 pn3 = x + 1.0;
|
|
|
478 pn4 = x * b;
|
|
|
479 value = pn3 / pn4;
|
|
|
480
|
|
|
481 for ( ; ; )
|
|
|
482 {
|
|
|
483 a = a + 1.0;
|
|
|
484 b = b + 2.0;
|
|
|
485 c = c + 1.0;
|
|
|
486 an = a * c;
|
|
|
487 pn5 = b * pn3 - an * pn1;
|
|
|
488 pn6 = b * pn4 - an * pn2;
|
|
|
489
|
|
|
490 if ( pn6 != 0.0 )
|
|
|
491 {
|
|
|
492 rn = pn5 / pn6;
|
|
|
493
|
|
|
494 if ( fabs ( value - rn ) <= r8_min ( tol, tol * rn ) )
|
|
|
495 {
|
|
|
496 break;
|
|
|
497 }
|
|
|
498 value = rn;
|
|
|
499 }
|
|
|
500
|
|
|
501 pn1 = pn3;
|
|
|
502 pn2 = pn4;
|
|
|
503 pn3 = pn5;
|
|
|
504 pn4 = pn6;
|
|
|
505 //
|
|
|
506 // Re-scale terms in continued fraction if terms are large.
|
|
|
507 //
|
|
|
508 if ( oflo <= fabs ( pn5 ) )
|
|
|
509 {
|
|
|
510 pn1 = pn1 / oflo;
|
|
|
511 pn2 = pn2 / oflo;
|
|
|
512 pn3 = pn3 / oflo;
|
|
|
513 pn4 = pn4 / oflo;
|
|
|
514 }
|
|
|
515 }
|
|
|
516
|
|
|
517 arg = arg + log ( value );
|
|
|
518
|
|
|
519 if ( elimit <= arg )
|
|
|
520 {
|
|
|
521 value = 1.0 - exp ( arg );
|
|
|
522 }
|
|
|
523 else
|
|
|
524 {
|
|
|
525 value = 1.0;
|
|
|
526 }
|
|
|
527 }
|
|
|
528
|
|
|
529 return value;
|
|
|
530 }
|
|
|
531
|
|
|
532 double r8_min ( double x, double y )
|
|
|
533
|
|
|
534 //****************************************************************************80
|
|
|
535 //
|
|
|
536 // Purpose:
|
|
|
537 //
|
|
|
538 // R8_MIN returns the minimum of two R8's.
|
|
|
539 //
|
|
|
540 // Licensing:
|
|
|
541 //
|
|
|
542 // This code is distributed under the GNU LGPL license.
|
|
|
543 //
|
|
|
544 // Modified:
|
|
|
545 //
|
|
|
546 // 31 August 2004
|
|
|
547 //
|
|
|
548 // Author:
|
|
|
549 //
|
|
|
550 // John Burkardt
|
|
|
551 //
|
|
|
552 // Parameters:
|
|
|
553 //
|
|
|
554 // Input, double X, Y, the quantities to compare.
|
|
|
555 //
|
|
|
556 // Output, double R8_MIN, the minimum of X and Y.
|
|
|
557 //
|
|
|
558 {
|
|
|
559 double value;
|
|
|
560
|
|
|
561 if ( y < x )
|
|
|
562 {
|
|
|
563 value = y;
|
|
|
564 }
|
|
|
565 else
|
|
|
566 {
|
|
|
567 value = x;
|
|
|
568 }
|
|
|
569 return value;
|
|
|
570 }
|
|
|
571
|
|
|
572 // This function is implemented by Li Song.
|
|
|
573 // If Y follows (theta, k)=(1/alpha, k)-gamma distribution, then P_Y(Y<t)=gammad(alpha*t, k).
|
|
|
574 // Furthurmore, chi-square with d.f. k is gamma distribution (2, k/2)
|
|
|
575 double chicdf( double x, double df )
|
|
|
576 {
|
|
|
577 int ifault ;
|
|
|
578 return gammad( x / 2.0, df / 2.0, &ifault ) ;
|
|
|
579 }
|
|
|
580
|
