1 /* -*- mode: C; mode: fold; -*- */
2
3 /*
4 Copyright (c) 2007 Massachusetts Institute of Technology
5
6 This software was developed by the MIT Center for Space Research
7 under contract SV1-61010 from the Smithsonian Institution.
8
9 Permission to use, copy, modify, distribute, and sell this software
10 and its documentation for any purpose is hereby granted without fee,
11 provided that the above copyright notice appear in all copies and
12 that both that copyright notice and this permission notice appear in
13 the supporting documentation, and that the name of the Massachusetts
14 Institute of Technology not be used in advertising or publicity
15 pertaining to distribution of the software without specific, written
16 prior permission. The Massachusetts Institute of Technology makes
17 no representations about the suitability of this software for any
18 purpose. It is provided "as is" without express or implied warranty.
19
20 THE MASSACHUSETTS INSTITUTE OF TECHNOLOGY DISCLAIMS ALL WARRANTIES
21 WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
22 MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE MASSACHUSETTS
23 INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY SPECIAL, INDIRECT OR
24 CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
25 OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
26 NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
27 WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
28 */
29
30 /* Author: John E. Davis <jed@jedsoft.org> */
31
32 #include "config.h"
33 #include <stdio.h>
34 #include <slang.h>
35 #include <math.h>
36 #include <string.h>
37
38 #include "stats-module.h"
39
40 #ifndef PI
41 #define PI 3.14159265358979323846264338327950288
42 #endif
43 /* <MODULE_INCLUDES> */
44
45 #ifdef __cplusplus
46 extern "C"
47 {
48 #endif
49 SLANG_MODULE(stats);
50 #ifdef __cplusplus
51 }
52 #endif
53
54 /* typedef unsigned long long uint64_t;
55 * typedef long long int64_t;
56 */
57
58 #define MODULE_VERSION_STRING "1.1.0"
59 #define MODULE_VERSION_NUMBER 10100
60
61 #if 0
62 double anderson_darling_statistic (double *f, unsigned int num)
63 {
64 double s;
65 double n;
66
67 n = (double)num;
68 for (i = 0; i < n; i++)
69 {
70 ((2*i+1)/n)* (log(f[i]) + log(f[n-i]));
71 }
72
73 }
74 #endif
75
76 /* Computes (n k) = n!/(k! (n-k)!) */
compute_binomial_coeff(unsigned int n,unsigned int k)77 static double compute_binomial_coeff (unsigned int n, unsigned int k) /*{{{*/
78 {
79 double f;
80 unsigned int i;
81
82 if (k > n)
83 return 0.0;
84 if ((k == 0) || (k == n))
85 return 1.0;
86
87 if (n-k < k)
88 k = n-k;
89
90 f = n;
91 for (i = 2; i <= k; i++)
92 {
93 n--;
94 f = (f/i)*n;
95 }
96 return f;
97 }
98
99 /*}}}*/
100
101 /* This function implements the _asymptotic_ form of the kolmogorov distribution
102 * for large N. Asymptotically, this is given by the series expansion:
103 *
104 * S(x) = 1 + 2\sum_{j=1}^{\infty} (-1)^j \exp(-2j^2 x^2)
105 *
106 * This is equivalent to
107 *
108 * S(x) = \sqrt{2\pi}/x \sum_{j=1}^\infty \exp[-(2j-1)^2\pi^2/(8x^2)]
109 *
110 * See eq. 1.4 of W. Feller, "On the Kolmogorov-Smirnov limit theorems for
111 * empirical distributions", Annals of Math. Stat, Vol 19 (1948), pp. 177-190.
112 *
113 * S(x) represents the prob of obtaining a (normalized) statistic <= x.
114 *
115 * I chose the value of 1.09 was chosen based upon some emperical testing. Of
116 * course this value will most likely vary with compiler.
117 */
smirnov_cdf(double x)118 static double smirnov_cdf (double x) /*{{{*/
119 {
120 int j;
121 unsigned int max_loops;
122 double a, b, sum;
123 max_loops = 5000;
124
125 if (x <= 0.15)
126 {
127 if (x < 0.0)
128 {
129 SLang_set_error (SL_INVALID_PARM);
130 return -1.0;
131 }
132 return 0.0;
133 }
134
135 if (x > 1.09)
136 {
137 int j3;
138
139 if (x > 19.4)
140 return 1.0;
141
142 x = 2*x*x;
143 sum = 0.0;
144 j = 1;
145 j3 = 3;
146 while (max_loops)
147 {
148 double dsum = exp(-j*j*x) * (1.0 - exp (-j3*x));
149 sum += dsum;
150 if (dsum == 0)
151 return 1-2*sum;
152
153 j += 2;
154 j3 += 4;
155 max_loops--;
156 }
157 return 1.0;
158 }
159
160 /* Use sqrt(X) = 0.5*log(X) */
161 b = log(sqrt(2*PI)/x);
162 a = (PI*PI/(8.0*x*x));
163
164 sum = 0.0;
165 j = 1;
166 while (max_loops)
167 {
168 double dsum = exp (-a*j*j + b);
169 sum += dsum;
170 if (dsum == 0.0)
171 return sum;
172 j += 2;
173 max_loops--;
174 }
175 return 0.0;
176 }
177
178 /*}}}*/
179
180 /* This gives the probability that an observed value of D_mn <= c/(mn).
181 */
kim_jennrich_cdf(unsigned int m,unsigned int n,unsigned int c)182 static double kim_jennrich_cdf (unsigned int m, unsigned int n, unsigned int c) /*{{{*/
183 {
184 unsigned int i, j;
185 double p;
186 double *u;
187
188 if (m > n)
189 {
190 unsigned int tmp = m;
191 m = n;
192 n = tmp;
193 }
194
195 u = (double *)SLmalloc (sizeof (double)*(n+1));
196 if (u == NULL)
197 return -1.0;
198
199 u[0] = 1.0;
200 for (j = 1; j <= n; j++)
201 {
202 u[j] = 1.0;
203 if (m*j > c)
204 u[j] = 0.0;
205 }
206
207 for (i = 1; i <= m; i++)
208 {
209 double w = i/(i + (double)n);
210 unsigned int ni = n*i;
211
212 u[0] = w*u[0];
213 if (ni > c) u[0] = 0;
214 for (j = 1; j <= n; j++)
215 {
216 unsigned int mj = m*j;
217 unsigned int diff;
218
219 diff = (ni >= mj) ? (ni-mj) : (mj-ni);
220 u[j] = (diff > c) ? 0.0 : u[j-1]+u[j]*w;
221 }
222 }
223 p = u[n];
224 if (p > 1.0) p = 1.0;
225 else if (p < 0.0) p = 0.0;
226
227 SLfree ((char *)u);
228 return p;
229 }
230
231 /*}}}*/
232
233 /* P(X<=s) */
mann_whitney_cdf(unsigned int m,unsigned int n,unsigned int s)234 static double mann_whitney_cdf (unsigned int m, unsigned int n, unsigned int s) /*{{{*/
235 {
236 unsigned int M, p;
237 unsigned int m_plus_n, t, u;
238 double *f, c, cf;
239
240 /* min allowed s: 1+2+...m = m(m+1)/2 = min
241 * max allowed s: (n+1)+(n+2)+ .. (n+m) = n*m + min
242 * avg = (min+max)/2 = nm/2 + m(m+1)/2 = (m/2)[n+ m+1]
243 */
244 M = m*(m+1)/2;
245 if (s < M)
246 return 0.0;
247 M += m*n;
248 if (M <= s)
249 return 1.0;
250
251 M = m*n/2;
252
253 f = (double *)SLmalloc ((M+1)*sizeof(double));
254 if (f == NULL)
255 return -1.0;
256 f[0] = 1.0;
257 for (p = 1; p <= M; p++)
258 f[p] = 0.0;
259
260 m_plus_n = m + n;
261 if (n + 1 < M)
262 {
263 p = m_plus_n;
264 if (p > M) p = M;
265 for (t = n + 1; t <= p; t++)
266 {
267 for (u = M; u >= t; u--)
268 f[u] -= f[u-t];
269 }
270 }
271 p = m;
272 if (M < p) p = M;
273 for (t = 1; t <= p; t++)
274 {
275 for (u = t; u <= M; u++)
276 f[u] += f[u-t];
277 }
278
279 c = compute_binomial_coeff (m+n, m);
280 cf = 0.0;
281 for (t = 0; t <= M; t++)
282 {
283 cf += f[t]/c;
284 f[t] = cf;
285 }
286
287 s -= m*(m+1)/2;
288 if (s > M)
289 {
290 s = m*n - s;
291 c = 1.0 - f[s];
292 }
293 else c = f[s];
294 SLfree ((char *) f);
295 return c;
296 }
297
298 /*}}}*/
299
300 /*{{{ Gamma Function */
301
302 #ifdef HAVE_LGAMMA
303 # define JDMlog_gamma lgamma
304 #else
305 /* The following code is adapted from my JDMath library */
306
307 /* This implementation is derived from
308 * Spouge JL. Computation of the gamma, digamma, and trigamma
309 * functions. SIAM J Numerical Anal 1994;31:931-44.
310 * as pointed out by Glynne Casteel <glynnec@ix.netcom.com> in the
311 * sci.math.num-analysis article <6potip$l56@sjx-ixn4.ix.netcom.com>.
312 */
313
314 #define NUM_COEFFS 18 /* make this even */
315 static double Param = NUM_COEFFS + 1; /* 0.840 */
316 /* determined empirically using the driver in main at end of the file. */
317
318 static double Coeffs[NUM_COEFFS+1];
319
320 static int Coeffs_Initialized = 0;
321
init_coefficients(void)322 static void init_coefficients (void)
323 {
324 int i;
325 double e = 2.718281828459045235360287471352;
326 double f = 2.506628274631000502415765284809; /* sqrt(2pi) */
327
328 Coeffs[0] = exp (-Param)*f;
329 Coeffs[1] = sqrt (Param-1)/e;
330 for (i = 1; i < NUM_COEFFS; i++)
331 {
332 register double x = Param - i;
333
334 Coeffs[i+1] = Coeffs[i]
335 * ((x - 1)*pow (1-1/x,i-0.5)/(i*e));
336 }
337 Coeffs_Initialized = 1;
338 }
339
JDMlog_gamma(double x)340 static double JDMlog_gamma (double x)
341 {
342 register double sum;
343 unsigned int i;
344
345 if (Coeffs_Initialized == 0)
346 init_coefficients ();
347
348 x -= 1.0;
349
350 sum = Coeffs[0];
351 i = 1;
352 while (i < NUM_COEFFS)
353 {
354 register double dsum;
355
356 dsum = Coeffs[i]/(x+i);
357 i++;
358 sum += (dsum - Coeffs[i]/(x+i));
359 i++;
360 }
361
362 return log(sum) + ((x+0.5)*log(x+Param) - x);
363 }
364 #endif /* HAVE_LGAMMA */
365
366 /* See A&S 6.5.4 for a definition. Here, the series expansion given by
367 * A&S 6.5.29 is used:
368 *
369 * gamma_star(a,x) = exp(-x) \sum_n x^n/Gamma(a+n+1)
370 *
371 * Here, Gamma(a+1)=a*Gamma(a), Gamma(a+2)=(a+1)*a*Gamma(a), ....
372 * Thus,
373 *
374 * gamma_star(a,x) = exp(-x)/Gamma(a) \sum_n x^n c_n
375 *
376 * where c_n = 1/(a(a+1)...(a+n)). For a > 0, c_n --> 0, hence based upon
377 * this c_n, the sum will terminate. However, x^n may overflow before this
378 * happens. However, as long as |x| < a, this should not happen. So avoid
379 * this function if |x| > a.
380 */
JDMlog_gamma_star(double a,double x)381 static double JDMlog_gamma_star (double a, double x)
382 {
383 double sum;
384 int n;
385 double cnxn;
386 double eps = 2.220446049250313e-16;
387
388 if (a == 0.0)
389 return 0.0;
390
391 if (a < 0.0)
392 {
393 /* FIXME!!! Handle a < 0 */
394 }
395
396 cnxn = 1.0 / a;
397 sum = cnxn;
398 n = 0;
399 while (n < 5000)
400 {
401 n++;
402 cnxn *= x/(a+n);
403 if (cnxn < sum*eps)
404 break;
405 sum += cnxn;
406 }
407
408 return log (sum) - x - JDMlog_gamma (a);
409 }
410
411 /* This incomplete gamma function has a continued fraction approximation
412 * given by A&S 6.5.31. This expansion may be used to compute the incomplete
413 * Gamma function (A&S 6.5.1) via A&S 6.5.2,3:
414 * Gamma(a,x) = Gamma(a) - P(a,x) Gamma(a)
415 * = Gamma(a) - gamma(a,x)
416 * or
417 * P(a,x) = 1 - Gamma(a,x)/Gamma(a).
418 *
419 * The continued fraction expansion (A&S 6.5.31) is given by:
420 *
421 * Gamma(a,x) = exp(-x)x^a (1/x+ (1-a)/1+ 1/x+ (2-a)/1+ 2/x+ ...)
422 *
423 * Use the recursion relations in theorem 2 of A&S 3.10.1 to evaluate this.
424 * That is, let f_n be defined as
425 *
426 * f_n = b_0 + a_1/b_1+ a_2/b_2+...a_n/b_n := A_n/B_n
427 *
428 * Then:
429 *
430 * A_n = b_n A_{n-1} + a_n A_{n-2}
431 * B_n = b_n B_{n-1} + a_n B_{n-2}
432 *
433 * where A_{-1}=1, A_0=b_0, B_{-1}=0, B_0=1
434 */
JDMlog_CapGamma(double a,double x)435 static double JDMlog_CapGamma (double a, double x)
436 {
437 double a1, a0, b0, b1;
438 double f;
439 double renorm_factor;
440 double eps = 2.220446049250313e-16;
441 int n;
442
443 if (x <= 0.0)
444 return log(x); /* NaN/-Inf */
445
446 /* Set up the recursion at the initial 1/x+ piece */
447 a1 = 1;
448 b1 = x;
449 a0 = 0;
450 b0 = 1;
451
452 renorm_factor = 1.0/b1;
453 f = a1 * renorm_factor;
454 n = 1;
455
456 if (renorm_factor != 0.0) while (n < 5000)
457 {
458 double f1, aa;
459
460 /* Note that the renormalization factor is 1/b1.
461 * So, replace renorm_factor*b1 combinations by 1 */
462
463 aa = n - a;
464 /* Now the (n-a)/(1+??) piece */
465 a0 = (a1 + aa * a0) * renorm_factor;
466 b0 = (b1 + aa * b0) * renorm_factor;
467
468 /* Since a0,b0 now have the renorm_factor applied, omit it below */
469 /* Handle {1,2,3...}/(x+??) piece */
470 a1 = x * a0 + n * a1 * renorm_factor;
471 /* b1 = x * b0 + n * b1 * renorm_factor;*/
472 b1 = x * b0 + n;
473
474 n++;
475
476 if (b1 == 0.0)
477 continue;
478
479 renorm_factor = 1.0 / b1;
480 f1 = a1 * renorm_factor;
481
482 if (fabs (f - f1) < eps * fabs(f1))
483 {
484 f = f1;
485 break;
486 }
487 f = f1;
488 }
489
490 return a * log(x) - x + log (f);
491 }
492
493 /* See A&S 6.5.1
494 * This is P(a,x) = gamma(a,x)/Gamma(a)
495 */
JDMincomplete_gamma(double a,double x)496 static double JDMincomplete_gamma (double a, double x)
497 {
498 if (a <= 0.0)
499 return log(a);
500
501 if (x >= a)
502 return 1.0 - exp (JDMlog_CapGamma (a,x) - JDMlog_gamma (a));
503
504 return exp (a * log(x) + JDMlog_gamma_star (a, x));
505 }
506
507 /*}}}*/
508
509 /*{{{ Incomplete Beta */
510
incbeta_cfe(double x,double a,double b,double * result)511 static int incbeta_cfe (double x, double a, double b, double *result)
512 {
513 unsigned int m, maxm = 1024;
514 double c = a+b;
515 double alpha_neg2, alpha_neg1, beta_neg2, beta_neg1;
516 double eps = 1e-14;
517 double last_f;
518 double factor;
519
520 /* factor = a*log(x)+b*log(1.0-x); */
521 factor = a*log(x)+b*log1p(-x);
522 factor += JDMlog_gamma(a+b)-JDMlog_gamma(a) - JDMlog_gamma(b);
523 factor = exp (factor)/a;
524
525 /* F = 1/(1+(d1/(1+(d2/(1+....)))))
526 * F0 = 1/1
527 * F1 = 1/1+ d1
528 * F2 = 1/1+ d1/1+ d2
529 * F_n = alpha_n/beta_n
530 * alpha_n = alpha_{n-1} + d_n*alpha_{n-2}
531 * beta_n = beta_{n-1} + d_n*beta_{n-2}
532 * (d0 = 1)
533 * alpha_{-2} = 1, alpha{-1} = 0, alpha_{0} = 1, alpha_1 = 1,
534 * beta_{-2} = 0, beta_{-1} = 1, beta_{0} = 1, beta_1 = (1+d1),
535 * alpha_2 = 1+d2
536 * beta_2 = (1+d1+d2)
537 * alpha_2/beta_2 = 1/(1+d1/(1+d2)) ok
538 */
539
540 /* This loop computes results for d2, d3, ... So initialize for d1 term */
541 alpha_neg2 = 1.0; alpha_neg1 = 1.0;
542 beta_neg2 = 1.0; beta_neg1 = 1.0-c/(a+1)*x;
543 last_f = alpha_neg1/beta_neg1;
544
545 m = 1;
546 while (m < maxm)
547 {
548 double d,g;
549 unsigned int m2 = 2*m;
550
551 g = a+m2;
552
553 /* even d_2m */
554 d = (m*(b-m)/((g-1.0)*g))*x;
555 alpha_neg2 = alpha_neg1 + d*alpha_neg2;
556 beta_neg2 = beta_neg1 + d*beta_neg2;
557
558 /* odd d_{2m+1} */
559 d = -(((a+m)*(c+m))/(g*(g+1.0)))*x;
560 alpha_neg1 = alpha_neg2 + d*alpha_neg1;
561 beta_neg1 = beta_neg2 + d*beta_neg1;
562
563 /* Renormalize to prevent large alpha/beta */
564 alpha_neg1 /= beta_neg1;
565 alpha_neg2 /= beta_neg1;
566 beta_neg2 /= beta_neg1;
567 beta_neg1 = 1.0;
568
569 if (fabs(alpha_neg1 - last_f) < eps * fabs(alpha_neg1))
570 {
571 *result = alpha_neg1*factor;
572 return 0;
573 }
574 last_f = alpha_neg1;
575 m++;
576 }
577
578 *result = last_f*factor;
579 return -1;
580 }
581
incbeta(double x,double a,double b,double * result)582 static int incbeta (double x, double a, double b, double *result)
583 {
584 double cfe;
585 int status;
586
587 if ((x < 0.0) || (x > 1.0))
588 {
589 SLang_verror (SL_INVALID_PARM, "Domain error for x in incbeta");
590 *result = -1;
591 return -1;
592 }
593
594 if ((x == 0.0) || (x == 1.0))
595 {
596 *result = x;
597 return 0;
598 }
599
600 /* Use 26.5.8 of A&S --- Note that this condition is different from that
601 * of A&S --- GSL uses this one and it seems to work better.
602 */
603 if (x * (a+b+2.0) < (a+1.0))
604 {
605 status = incbeta_cfe (x, a, b, &cfe);
606 }
607 else
608 {
609 /* Use I_x(a,b) = 1.0 - I_{1-x}(b,a) */
610 status = incbeta_cfe (1.0-x, b, a, &cfe);
611 cfe = 1.0 - cfe;
612 }
613
614 *result = cfe;
615 return status;
616 }
617
618 /*}}}*/
619
620 /*{{{ Intrinsics */
621
622 /* See A&S 6.5.13
623 *
624 * This computes: \sum_{0<=k<=n} lambda^k/k! exp(-lambda)
625 * From A&S 6.5.13, this is just 1-P(n+1,lambda)
626 */
poisson_cdf_intrin(double * lambdap,int * np)627 static double poisson_cdf_intrin (double *lambdap, int *np)
628 {
629 int n = *np;
630 double lambda = *lambdap;
631
632 if (n < 0)
633 return 0.0;
634
635 n = n+1;
636 if (lambda > 1000.0)
637 {
638 double dn = sqrt(n);
639 if (fabs(lambda - n) < dn)
640 {
641 /* Wilson and Hilferty, 1951 */
642 double c = pow (lambda/n, 1.0/3.0);
643 double mu = 1.0-1.0/(9.0*n);
644 double s = 1.0/(3*dn);
645
646 return 0.5 * (1 - erf((c-mu)/s/sqrt(2.0)));
647 }
648 /* drop */
649 }
650
651
652 return 1.0 - JDMincomplete_gamma ((double)n, lambda);
653 }
654
655 /* P(X^2 < chisqr)
656 * A&S 26.4.19
657 */
chisqr_cdf_intrin(int * dofp,double * chisqrp)658 static double chisqr_cdf_intrin (int *dofp, double *chisqrp)
659 {
660 int dof = *dofp;
661 double chisqr = *chisqrp;
662
663 if (dof <= 0)
664 {
665 SLang_verror (SL_INVALID_PARM, "The number of degrees of freedom should be positive");
666 return -1.0;
667 }
668 if (chisqr < 0)
669 {
670 SLang_verror (SL_INVALID_PARM, "Expecting a non-negative value for the chi-square statistic");
671 return -1.0;
672 }
673 return JDMincomplete_gamma (0.5*dof, 0.5*chisqr);
674 }
675
smirnov_cdf_intrin(double * x)676 static double smirnov_cdf_intrin (double *x)
677 {
678 return smirnov_cdf (*x);
679 }
680
kim_jennrich_cdf_intrin(unsigned int * m,unsigned int * n,unsigned int * c)681 static double kim_jennrich_cdf_intrin (unsigned int *m, unsigned int *n, unsigned int *c)
682 {
683 return kim_jennrich_cdf (*m, *n, *c);
684 }
685
mann_whitney_cdf_intrin(unsigned int * m,unsigned int * n,double * s)686 static double mann_whitney_cdf_intrin (unsigned int *m, unsigned int *n, double *s)
687 {
688 unsigned int sn = (unsigned int) (*s + 0.5);
689 return mann_whitney_cdf (*m, *n, sn);
690 }
691
692 /* df needs to be real valued for the welch test */
student_t_cdf_intrin(double * t,double * dfp)693 static double student_t_cdf_intrin (double *t, double *dfp)
694 {
695 double x = *t;
696 double n = *dfp;
697 double cdf;
698
699 (void) incbeta (n/(n+x*x), 0.5*n, 0.5, &cdf);
700 cdf *= 0.5;
701 if (x > 0)
702 cdf = 1.0 - cdf;
703 return cdf;
704 }
705
f_cdf_intrin(double * t,double * nu1p,double * nu2p)706 static double f_cdf_intrin (double *t, double *nu1p, double *nu2p)
707 {
708 double x = *t;
709 double cdf;
710 double nu1 = *nu1p, nu2 = *nu2p;
711
712 if (x < 0.0)
713 return 0.0;
714
715 (void) incbeta (nu2/(nu2+nu1*x), 0.5*nu2, 0.5*nu1, &cdf);
716 return 1.0-cdf;
717 }
718
719 /*}}}*/
720
721 /* Mean, Stdev, etc... */
722 #define GENERIC_TYPE double
723 #define MEAN_FUNCTION mean_doubles
724 #define MEAN_RESULT_TYPE double
725 #define STDDEV_FUNCTION stddev_doubles
726 #define STDDEV_RESULT_TYPE double
727 #define MEDIAN_FUNCTION median_doubles
728 #define NC_MEDIAN_FUNCTION nc_median_doubles
729 #include "stats-module.inc"
730
731 #define GENERIC_TYPE float
732 #define MEAN_FUNCTION mean_floats
733 #define MEAN_RESULT_TYPE float
734 #define STDDEV_FUNCTION stddev_floats
735 #define STDDEV_RESULT_TYPE float
736 #define MEDIAN_FUNCTION median_floats
737 #define NC_MEDIAN_FUNCTION nc_median_floats
738 #include "stats-module.inc"
739
740 #define GENERIC_TYPE int
741 #define MEAN_FUNCTION mean_ints
742 #define MEAN_RESULT_TYPE double
743 #define STDDEV_FUNCTION stddev_ints
744 #define STDDEV_RESULT_TYPE double
745 #define MEDIAN_FUNCTION median_ints
746 #define NC_MEDIAN_FUNCTION nc_median_ints
747 #include "stats-module.inc"
748
749 #define GENERIC_TYPE unsigned int
750 #define MEAN_FUNCTION mean_uints
751 #define MEAN_RESULT_TYPE double
752 #define STDDEV_FUNCTION stddev_uints
753 #define STDDEV_RESULT_TYPE double
754 #define MEDIAN_FUNCTION median_uints
755 #define NC_MEDIAN_FUNCTION nc_median_uints
756 #include "stats-module.inc"
757
758 #if SIZEOF_LONG != SIZEOF_INT
759 # define GENERIC_TYPE long
760 # define MEAN_FUNCTION mean_longs
761 # define MEAN_RESULT_TYPE double
762 # define STDDEV_FUNCTION stddev_longs
763 # define STDDEV_RESULT_TYPE double
764 # define MEDIAN_FUNCTION median_longs
765 # define NC_MEDIAN_FUNCTION nc_median_longs
766 # include "stats-module.inc"
767
768 # define GENERIC_TYPE unsigned long
769 # define MEAN_FUNCTION mean_ulongs
770 # define MEAN_RESULT_TYPE double
771 # define STDDEV_FUNCTION stddev_ulongs
772 # define STDDEV_RESULT_TYPE double
773 # define MEDIAN_FUNCTION median_ulongs
774 # define NC_MEDIAN_FUNCTION nc_median_ulongs
775 # include "stats-module.inc"
776 #else
777 # define mean_longs mean_ints
778 # define stddev_longs stddev_ints
779 # define median_longs median_ints
780 # define nc_median_longs nc_median_ints
781 # define mean_ulongs mean_uints
782 # define stddev_ulongs stddev_uints
783 # define median_ulongs median_uints
784 # define nc_median_ulongs nc_median_uints
785 #endif
786
787 #if SIZEOF_SHORT != SIZEOF_INT
788 # define GENERIC_TYPE short
789 # define MEAN_FUNCTION mean_shorts
790 # define MEAN_RESULT_TYPE float
791 # define STDDEV_FUNCTION stddev_shorts
792 # define STDDEV_RESULT_TYPE float
793 # define MEDIAN_FUNCTION median_shorts
794 # define NC_MEDIAN_FUNCTION nc_median_shorts
795 # include "stats-module.inc"
796
797 # define GENERIC_TYPE unsigned short
798 # define MEAN_FUNCTION mean_ushorts
799 # define MEAN_RESULT_TYPE float
800 # define STDDEV_FUNCTION stddev_ushorts
801 # define STDDEV_RESULT_TYPE float
802 # define MEDIAN_FUNCTION median_ushorts
803 # define NC_MEDIAN_FUNCTION nc_median_ushorts
804 # include "stats-module.inc"
805 #else
806 # define mean_shorts mean_ints
807 # define stddev_shorts stddev_ints
808 # define median_shorts median_ints
809 # define nc_median_shorts nc_median_ints
810 # define mean_ushorts mean_uints
811 # define stddev_ushorts stddev_uints
812 # define median_ushorts median_uints
813 # define nc_median_ushorts nc_median_uints
814 #endif
815
816 #define GENERIC_TYPE signed char
817 #define MEAN_FUNCTION mean_chars
818 #define MEAN_RESULT_TYPE float
819 #define STDDEV_FUNCTION stddev_chars
820 #define STDDEV_RESULT_TYPE float
821 #define MEDIAN_FUNCTION median_chars
822 #define NC_MEDIAN_FUNCTION nc_median_chars
823 #include "stats-module.inc"
824
825 #define GENERIC_TYPE unsigned char
826 #define MEAN_FUNCTION mean_uchars
827 #define MEAN_RESULT_TYPE float
828 #define STDDEV_FUNCTION stddev_uchars
829 #define STDDEV_RESULT_TYPE float
830 #define MEDIAN_FUNCTION median_uchars
831 #define NC_MEDIAN_FUNCTION nc_median_uchars
832 #include "stats-module.inc"
833
834 static SLarray_Contract_Type Stddev_Functions [] =
835 {
836 {SLANG_CHAR_TYPE, SLANG_CHAR_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) stddev_chars},
837 {SLANG_UCHAR_TYPE,SLANG_UCHAR_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) stddev_uchars},
838 {SLANG_SHORT_TYPE,SLANG_SHORT_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) stddev_shorts},
839 {SLANG_USHORT_TYPE, SLANG_USHORT_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) stddev_ushorts},
840 {SLANG_INT_TYPE, SLANG_INT_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) stddev_ints},
841 {SLANG_UINT_TYPE, SLANG_UINT_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) stddev_uints},
842 {SLANG_LONG_TYPE, SLANG_LONG_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) stddev_longs},
843 {SLANG_ULONG_TYPE, SLANG_ULONG_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) stddev_ulongs},
844 {SLANG_FLOAT_TYPE, SLANG_FLOAT_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) stddev_floats},
845 {SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) stddev_doubles},
846 {0, 0, 0, NULL}
847 };
848
stddev_intrin(void)849 static void stddev_intrin (void)
850 {
851 if (SLang_Num_Function_Args == 0)
852 {
853 SLang_verror (SL_Usage_Error, "x = stddev(X [,dim])");
854 return;
855 }
856 (void) SLarray_contract_array (Stddev_Functions);
857 }
858
859 static SLCONST SLarray_Contract_Type Mean_Functions [] =
860 {
861 {SLANG_CHAR_TYPE, SLANG_CHAR_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) mean_chars},
862 {SLANG_UCHAR_TYPE, SLANG_UCHAR_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) mean_uchars},
863 {SLANG_SHORT_TYPE, SLANG_SHORT_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) mean_shorts},
864 {SLANG_USHORT_TYPE, SLANG_USHORT_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) mean_ushorts},
865 {SLANG_UINT_TYPE, SLANG_UINT_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) mean_uints},
866 {SLANG_INT_TYPE, SLANG_INT_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) mean_ints},
867 {SLANG_LONG_TYPE, SLANG_LONG_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) mean_longs},
868 {SLANG_ULONG_TYPE, SLANG_ULONG_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) mean_ulongs},
869 {SLANG_FLOAT_TYPE, SLANG_FLOAT_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) mean_floats},
870 {SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) mean_doubles},
871 {0, 0, 0, NULL}
872 };
873
mean_intrin(void)874 static void mean_intrin (void)
875 {
876 if (SLang_Num_Function_Args == 0)
877 {
878 SLang_verror (SL_Usage_Error, "x = mean(X [,dim])");
879 return;
880 }
881 (void) SLarray_contract_array (Mean_Functions);
882 }
883
884 static SLarray_Contract_Type Median_Functions [] =
885 {
886 {SLANG_CHAR_TYPE, SLANG_CHAR_TYPE, SLANG_CHAR_TYPE, (SLarray_Contract_Fun_Type *) median_chars},
887 {SLANG_UCHAR_TYPE,SLANG_UCHAR_TYPE, SLANG_UCHAR_TYPE, (SLarray_Contract_Fun_Type *) median_uchars},
888 {SLANG_SHORT_TYPE,SLANG_SHORT_TYPE, SLANG_SHORT_TYPE, (SLarray_Contract_Fun_Type *) median_shorts},
889 {SLANG_USHORT_TYPE, SLANG_USHORT_TYPE, SLANG_USHORT_TYPE, (SLarray_Contract_Fun_Type *) median_ushorts},
890 {SLANG_INT_TYPE, SLANG_INT_TYPE, SLANG_INT_TYPE, (SLarray_Contract_Fun_Type *) median_ints},
891 {SLANG_UINT_TYPE, SLANG_UINT_TYPE, SLANG_UINT_TYPE, (SLarray_Contract_Fun_Type *) median_uints},
892 {SLANG_LONG_TYPE, SLANG_LONG_TYPE, SLANG_LONG_TYPE, (SLarray_Contract_Fun_Type *) median_longs},
893 {SLANG_ULONG_TYPE, SLANG_ULONG_TYPE, SLANG_ULONG_TYPE, (SLarray_Contract_Fun_Type *) median_ulongs},
894 {SLANG_FLOAT_TYPE, SLANG_FLOAT_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) median_floats},
895 {SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) median_doubles},
896 {0, 0, 0, NULL}
897 };
898
899 static SLarray_Contract_Type NC_Median_Functions [] =
900 {
901 {SLANG_CHAR_TYPE, SLANG_CHAR_TYPE, SLANG_CHAR_TYPE, (SLarray_Contract_Fun_Type *) nc_median_chars},
902 {SLANG_UCHAR_TYPE,SLANG_UCHAR_TYPE, SLANG_UCHAR_TYPE, (SLarray_Contract_Fun_Type *) nc_median_uchars},
903 {SLANG_SHORT_TYPE,SLANG_SHORT_TYPE, SLANG_SHORT_TYPE, (SLarray_Contract_Fun_Type *) nc_median_shorts},
904 {SLANG_USHORT_TYPE, SLANG_USHORT_TYPE, SLANG_USHORT_TYPE, (SLarray_Contract_Fun_Type *) nc_median_ushorts},
905 {SLANG_INT_TYPE, SLANG_INT_TYPE, SLANG_INT_TYPE, (SLarray_Contract_Fun_Type *) nc_median_ints},
906 {SLANG_UINT_TYPE, SLANG_UINT_TYPE, SLANG_UINT_TYPE, (SLarray_Contract_Fun_Type *) nc_median_uints},
907 {SLANG_LONG_TYPE, SLANG_LONG_TYPE, SLANG_LONG_TYPE, (SLarray_Contract_Fun_Type *) nc_median_longs},
908 {SLANG_ULONG_TYPE, SLANG_ULONG_TYPE, SLANG_ULONG_TYPE, (SLarray_Contract_Fun_Type *) nc_median_ulongs},
909 {SLANG_FLOAT_TYPE, SLANG_FLOAT_TYPE, SLANG_FLOAT_TYPE, (SLarray_Contract_Fun_Type *) nc_median_floats},
910 {SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, (SLarray_Contract_Fun_Type *) nc_median_doubles},
911 {0, 0, 0, NULL}
912 };
913
median_nc_intrin(void)914 static void median_nc_intrin (void)
915 {
916 if (SLang_Num_Function_Args == 0)
917 {
918 SLang_verror (SL_Usage_Error, "x = median_nc (X [,dim])");
919 return;
920 }
921 (void) SLarray_contract_array (NC_Median_Functions);
922 }
median_intrin(void)923 static void median_intrin (void)
924 {
925 if (SLang_Num_Function_Args == 0)
926 {
927 SLang_verror (SL_Usage_Error, "x = median (X [,dim])");
928 return;
929 }
930 (void) SLarray_contract_array (Median_Functions);
931 }
932
binomial_intrin(void)933 static void binomial_intrin (void)
934 {
935 unsigned int n, k, i;
936 double f;
937 SLindex_Type dims;
938 SLang_Array_Type *at;
939 double *data;
940
941 if (SLang_Num_Function_Args == 2)
942 {
943 if ((-1 == SLang_pop_uint (&k)) || (-1 == (SLang_pop_uint (&n))))
944 return;
945
946 (void) SLang_push_double (compute_binomial_coeff (n, k));
947 return;
948 }
949
950 if (-1 == SLang_pop_uint (&n))
951 return;
952
953 dims = (SLindex_Type) (n+1);
954 if (NULL == (at = SLang_create_array (SLANG_DOUBLE_TYPE, 0, NULL, &dims, 1)))
955 return;
956
957 data = (double *)at->data;
958 k = n;
959
960 data[0] = 1.0;
961 data[k] = 1.0;
962 f = 1.0;
963 for (i = 1; i <= k; i++)
964 {
965 f = (f/i)*k;
966 k--;
967 data[i] = data[k] = f;
968 }
969 (void) SLang_push_array (at, 1);
970 }
971
972 #if 0
973 static void erf_intrin (void)
974 {
975 int is_float;
976 SLang_Array_Type *at;
977
978 is_float == (SLang_peek_at_stack1 () == SLANG_FLOAT_TYPE);
979 if (SLang_peek_at_stack () != SLANG_ARRAY_TYPE)
980 {
981 double x;
982 if (is_float)
983 {
984 float xf;
985 if (-1 == SLang_pop_float (&xf))
986 return;
987 (void) SLang_push_float ((float) erf((double) xf));
988 }
989 if (0 == POP_DOUBLE (&x))
990 (void) SLang_push_double (x);
991 return;
992 }
993 if (is_float)
994 {
995 float *f,
996 if (-1 == SLang_pop_array_of_type (&at, SLANG_FLOAT_TYPE))
997 return;
998 }
999 xxxxxxxxxx
1000 }
1001 #endif
1002
normal_cdf_intrin(double * xp)1003 static double normal_cdf_intrin (double *xp)
1004 {
1005 double x = *xp;
1006 return 0.5 * (1.0 + erf (x/sqrt(2.0)));
1007 }
1008
kendall_tau_intrin(void)1009 static double kendall_tau_intrin (void)
1010 {
1011 double val, z;
1012 size_t n;
1013 SLang_Array_Type *at, *bt;
1014
1015 if (-1 == SLang_pop_array_of_type (&bt, SLANG_ARRAY_INDEX_TYPE))
1016 return -1.0;
1017 n = bt->num_elements;
1018
1019 if (-1 == SLang_pop_array_of_type (&at, SLANG_ARRAY_INDEX_TYPE))
1020 {
1021 SLang_free_array (bt);
1022 return -1.0;
1023 }
1024 if (at->num_elements != n)
1025 {
1026 SLang_verror (SL_TYPE_MISMATCH, "kendall_tau: arrays must have the same size");
1027 val = -1.0;
1028 goto free_and_return;
1029 }
1030
1031 val = _pSLstats_kendall_tau ((SLindex_Type *)at->data, (SLindex_Type *)bt->data, n, &z);
1032
1033 /* drop */
1034
1035 free_and_return:
1036
1037 SLang_free_array (at);
1038 SLang_free_array (bt);
1039 (void) SLang_push_double (z);
1040 return val;
1041 }
1042
1043 static SLang_Intrin_Fun_Type Module_Intrinsics [] =
1044 {
1045 MAKE_INTRINSIC_1("smirnov_cdf", smirnov_cdf_intrin, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE),
1046 MAKE_INTRINSIC_2("chisqr_cdf", chisqr_cdf_intrin, SLANG_DOUBLE_TYPE, SLANG_INT_TYPE, SLANG_DOUBLE_TYPE),
1047 MAKE_INTRINSIC_1("_normal_cdf", normal_cdf_intrin, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE),
1048 MAKE_INTRINSIC_2("student_t_cdf", student_t_cdf_intrin, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE),
1049 MAKE_INTRINSIC_2("_poisson_cdf", poisson_cdf_intrin, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, SLANG_INT_TYPE),
1050 MAKE_INTRINSIC_3("kim_jennrich_cdf", kim_jennrich_cdf_intrin, SLANG_DOUBLE_TYPE, SLANG_UINT_TYPE, SLANG_UINT_TYPE, SLANG_UINT_TYPE),
1051 MAKE_INTRINSIC_3("mann_whitney_cdf", mann_whitney_cdf_intrin, SLANG_DOUBLE_TYPE, SLANG_UINT_TYPE, SLANG_UINT_TYPE, SLANG_DOUBLE_TYPE),
1052 MAKE_INTRINSIC_3("f_cdf", f_cdf_intrin, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE, SLANG_DOUBLE_TYPE),
1053 MAKE_INTRINSIC_0("mean", mean_intrin, SLANG_VOID_TYPE),
1054 MAKE_INTRINSIC_0("stddev", stddev_intrin, SLANG_VOID_TYPE),
1055 MAKE_INTRINSIC_0("median", median_intrin, SLANG_VOID_TYPE),
1056 MAKE_INTRINSIC_0("median_nc", median_nc_intrin, SLANG_VOID_TYPE),
1057 MAKE_INTRINSIC_0("binomial", binomial_intrin, SLANG_VOID_TYPE),
1058 MAKE_INTRINSIC_0("_kendall_tau", kendall_tau_intrin, SLANG_DOUBLE_TYPE),
1059 SLANG_END_INTRIN_FUN_TABLE
1060 };
1061
1062 static SLFUTURE_CONST char *Module_Version_String = MODULE_VERSION_STRING;
1063 static SLang_Intrin_Var_Type Module_Variables [] =
1064 {
1065 MAKE_VARIABLE("_stats_module_version_string", &Module_Version_String, SLANG_STRING_TYPE, 1),
1066 SLANG_END_INTRIN_VAR_TABLE
1067 };
1068
1069 static SLang_IConstant_Type Module_IConstants [] =
1070 {
1071 MAKE_ICONSTANT("_stats_module_version", MODULE_VERSION_NUMBER),
1072 SLANG_END_ICONST_TABLE
1073 };
1074
1075 static SLang_DConstant_Type Module_DConstants [] =
1076 {
1077 SLANG_END_DCONST_TABLE
1078 };
1079
init_stats_module_ns(char * ns_name)1080 int init_stats_module_ns (char *ns_name)
1081 {
1082 SLang_NameSpace_Type *ns = SLns_create_namespace (ns_name);
1083 if (ns == NULL)
1084 return -1;
1085
1086 if (
1087 (-1 == SLns_add_intrin_var_table (ns, Module_Variables, NULL))
1088 || (-1 == SLns_add_intrin_fun_table (ns, Module_Intrinsics, NULL))
1089 || (-1 == SLns_add_iconstant_table (ns, Module_IConstants, NULL))
1090 || (-1 == SLns_add_dconstant_table (ns, Module_DConstants, NULL))
1091 )
1092 return -1;
1093
1094 return 0;
1095 }
1096
1097 /* This function is optional */
deinit_stats_module(void)1098 void deinit_stats_module (void)
1099 {
1100 }
1101