Mercurial > repos > lsong10 > psiclass

comparison PsiCLASS-1.0.2/stats.cpp @ 0:903fc43d6227 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author lsong10
date Fri, 26 Mar 2021 16:52:45 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:903fc43d6227
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