1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 /* pvalues.c - routines relating to computation of pvalues
21 of sample statistics
22 */
23
24 #include "libgretl.h"
25 #include "libset.h"
26 #include "../../cephes/libprob.h"
27
28 #include <errno.h>
29
30 #define NORM_CDF_MAX 0.9999999999999999
31
32 /**
33 * SECTION:pvalues
34 * @short_description: probability values for test statistics and
35 * related functionality
36 * @title: P-values
37 * @include: libgretl.h
38 *
39 * Libgretl uses the cephes library, developed by Stephen
40 * Moshier, as the basic engine for much of the functionality
41 * herein. We add some extra distributions, and wrap the cephes
42 * functions for ease of use with libgretl (e.g. on failure
43 * they return the libgretl missing value code, %NADBL).
44 */
45
46 static int nct_pdf_array (double df, double delta, double *x, int n);
47
48 /**
49 * gammafun:
50 * @x: argument.
51 *
52 * Returns: the gamma function of @x, or #NADBL on failure.
53 */
54
gammafun(double x)55 double gammafun (double x)
56 {
57 double ret = cephes_gamma(x);
58
59 if (get_cephes_errno()) {
60 ret = NADBL;
61 }
62
63 return ret;
64 }
65
66 /**
67 * lngamma:
68 * @x: argument.
69 *
70 * Returns: the log gamma function of @x, or #NADBL on failure.
71 */
72
lngamma(double x)73 double lngamma (double x)
74 {
75 double ret = cephes_lgamma(x);
76
77 if (get_cephes_errno()) {
78 ret = NADBL;
79 }
80
81 return ret;
82 }
83
84 /**
85 * digamma:
86 * @x: argument.
87 *
88 * Returns: the digamma (or Psi) function of @x, or #NADBL on failure.
89 */
90
digamma(double x)91 double digamma (double x)
92 {
93 double ret = psi(x);
94
95 if (get_cephes_errno()) {
96 ret = NADBL;
97 }
98
99 return ret;
100 }
101
102 /**
103 * trigamma:
104 * @x: argument.
105 *
106 * Returns: the trigamma function of @x, or #NADBL on failure. The code
107 * is adapted from
108 * https://people.sc.fsu.edu/~jburkardt/f_src/asa121/asa121.html
109 * See BE Schneider, Algorithm AS 121: Trigamma Function.
110 * Applied Statistics, Volume 27, Number 1, pages 97-99, 1978.
111 *
112 * The main modification with respect to the published version is the
113 * addition of three extra terms to the asymptotic expansion for x >= B.
114 */
115
trigamma(double x)116 double trigamma (double x)
117 {
118 double ret = 0;
119 double A = 0.000001; /* threshold for "small" argument */
120 double B = 5.0; /* threshold for "large" argument */
121
122 /* the Bernoulli numbers */
123 double b2 = 1.0/6;
124 double b4 = -1.0/30;
125 double b6 = 1.0/42;
126 double b8 = -1.0/30;
127 double b10 = 5.0/66;
128 double b12 = -691.0/2730;
129 double b14 = 7.0/6;
130
131 if (x <= 0) {
132 ret = NADBL;
133 } else if (x < A) {
134 ret = 1.0 / (x * x);
135 } else {
136 double y, a1, a2, z = x;
137
138 while (z < B) {
139 ret += 1.0 / (z * z);
140 z += 1.0;
141 }
142 /* Apply asymptotic formula for argument >= B */
143 y = 1.0 / (z * z);
144 a1 = 0.5 * y;
145 a2 = (1 + y *
146 (b2 + y *
147 (b4 + y *
148 (b6 + y *
149 (b8 + y *
150 (b10 + y *
151 (b12 + y * b14))))))) / z;
152 ret += a1 + a2;
153 }
154
155 return ret;
156 }
157
158 /**
159 * hypergeo:
160 * @a: argument.
161 * @b: argument.
162 * @c: argument.
163 * @x: absolute value must be less than 1.0.
164 *
165 * Returns: the Gauss hypergeometric function 2F1 of
166 * the given arguments.
167 */
168
hypergeo(double a,double b,double c,double x)169 double hypergeo (double a, double b, double c, double x)
170 {
171 double ret = hyp2f1(a, b, c, x);
172
173 if (get_cephes_errno()) {
174 ret = NADBL;
175 }
176
177 return ret;
178 }
179
180 /**
181 * binomial_cdf:
182 * @p: probability of success on each trial.
183 * @n: number of trials.
184 * @k: maximum number of successes.
185 *
186 * Returns: the probability of @k or less successes on
187 * @n trials given binomial probability @p, or
188 * #NADBL on failure.
189 */
190
binomial_cdf(double p,int n,int k)191 double binomial_cdf (double p, int n, int k)
192 {
193 double x = NADBL;
194
195 if (p >= 0 && n >= 0 && k >= 0) {
196 x = bdtr(k, n, p);
197 if (get_cephes_errno()) {
198 x = NADBL;
199 }
200 }
201
202 return x;
203 }
204
205 /**
206 * beta_cdf:
207 *
208 * Returns the probability that a B(a,b) random variable is
209 * between 0 and z, or #NADBL on failure.
210 */
211
beta_cdf(double a,double b,double z)212 double beta_cdf (double a, double b, double z)
213 {
214 double x = NADBL;
215
216 if (a > 0 && b > 0 && z >= 0 && z <= 1) {
217 x = btdtr(a, b, z);
218 if (get_cephes_errno()) {
219 x = NADBL;
220 }
221 }
222
223 return x;
224 }
225
226 /**
227 * binomial_cdf_comp:
228 * @p: probability of success on each trial.
229 * @n: number of trials.
230 * @k: maximum number of successes.
231 *
232 * Returns: the probability of @k + 1 or more successes on
233 * @n trials given binomial probability @p, or
234 * #NADBL on failure.
235 */
236
binomial_cdf_comp(double p,int n,int k)237 double binomial_cdf_comp (double p, int n, int k)
238 {
239 double x = NADBL;
240
241 if (p >= 0 && n >= 0 && k >= 0) {
242 x = bdtrc(k, n, p);
243 if (get_cephes_errno()) {
244 x = NADBL;
245 }
246 }
247
248 return x;
249 }
250
251 /* Binomial inverse CDF via binary search */
252
bininv_binary(double p,int n,double u)253 static int bininv_binary (double p, int n, double u)
254 {
255 int klo = 0;
256 int khi = n;
257 int kmid, ke, ke1;
258 double f, f1;
259 double fe, fe1;
260 double mu, pdmax;
261 int evals = 0;
262 int nk, emax;
263 int ret = -1;
264
265 if (u > 0.5) {
266 klo = floor(n*p) - 1;
267 } else if (u < 0.5) {
268 khi = floor(n*p) + 1;
269 }
270
271 kmid = (klo + khi) / 2;
272 nk = khi - klo + 1;
273 emax = 2 * log2(nk);
274
275 /* Find the max step-up of the CDF, which occurs
276 from expected value minus 1 to expected value.
277 */
278 mu = n*p;
279 ke = u > 0.5 ? ceil(mu) : floor(mu);
280 if (ke > 0) {
281 ke1 = ke - 1;
282 fe = binomial_cdf(p, n, ke);
283 fe1 = binomial_cdf(p, n, ke1);
284 pdmax = fe - fe1;
285 } else {
286 pdmax = 0;
287 ke = ke1 = -1;
288 fe = fe1 = 0;
289 }
290
291 while (evals < emax) {
292 if (kmid == ke || kmid == ke1) {
293 f = (kmid == ke)? fe : fe1;
294 } else {
295 f = binomial_cdf(p, n, kmid);
296 evals++;
297 }
298 if (f < u) {
299 /* we need to look on the right */
300 klo = kmid;
301 kmid = (klo + khi) / 2;
302 if (kmid == klo) {
303 ret = kmid;
304 break;
305 }
306 } else if (f == u) {
307 /* on the nail (unlikely) */
308 ret = kmid;
309 break;
310 } else {
311 /* f > u */
312 if (f - u <= pdmax && kmid > 0) {
313 /* check out k - 1 */
314 int k1 = kmid - 1;
315
316 if (k1 == ke || k1 == ke1) {
317 f1 = (k1 == ke)? fe : fe1;
318 } else {
319 f1 = binomial_cdf(p, n, k1);
320 evals++;
321 }
322 if (f1 < u) {
323 ret = kmid;
324 break;
325 }
326 }
327 /* we need to look on the left */
328 khi = kmid;
329 kmid = (klo + khi) / 2;
330 if (kmid == klo) {
331 ret = kmid;
332 break;
333 }
334 }
335 }
336
337 if (ret < 0) {
338 fprintf(stderr, "bininv_binary, bad: evals=%d, p=%g, n=%d, u=%g\n",
339 evals, p, n, u);
340 }
341
342 return ret;
343 }
344
345 /* Binomial inverse CDF via bottom-up or top-down summation */
346
bininv_sum(double p,int n,double u)347 static int bininv_sum (double p, int n, double u)
348 {
349 double f, a, s;
350 int k = 1;
351
352 f = binomial_cdf(p, n, n/2);
353
354 if (u <= f) {
355 /* bottom-up */
356 a = pow(1 - p, n);
357 s = a - u;
358 while (s < 0) {
359 a = (a*p/(1 - p)) * (n - k + 1)/k;
360 s += a;
361 k++;
362 }
363 return k - 1;
364 } else {
365 /* top-down */
366 a = pow(p, n);
367 s = 1 - u - a;
368 while (s >= 0) {
369 a = (a*(1 - p)/p) * (n - k + 1)/k;
370 s -= a;
371 k++;
372 }
373 return n - k + 1;
374 }
375
376 return -1;
377 }
378
379 /* Returns the smallest x such that the probability of obtaining
380 x successes on @n trials with success probability @p is at least
381 the given cumulative probability @u. FIXME efficiency?
382 */
383
binomial_cdf_inverse(double p,int n,double u)384 static double binomial_cdf_inverse (double p, int n, double u)
385 {
386 double x = NADBL;
387 int k;
388
389 if (u <= 0 || u > 1 || n <= 0 || p <= 0 || p >= 1) {
390 ;
391 } else if (u == 1.0) {
392 x = n;
393 } else if (n < 500) {
394 k = bininv_sum(p, n, u);
395 if (k >= 0) {
396 x = k;
397 }
398 } else {
399 k = bininv_binary(p, n, u);
400 if (k >= 0) {
401 x = k;
402 }
403 }
404
405 return x;
406 }
407
408 #if 0
409
410 /* Returns the event probability p such that the Binomial CDF
411 for @k successes on @n trials equals @a.
412 */
413
414 static double gretl_bdtri (int n, int k, double a)
415 {
416 double p = NADBL;
417
418 if (a > 0 && a < 1 && n >= 0 && k >= 0 && k <= n) {
419 p = bdtri(k, n, a);
420 if (get_cephes_errno()) {
421 p = NADBL;
422 }
423 }
424
425 return p;
426 }
427
428 #endif
429
430 /* The following is no doubt horribly inefficient */
431
binomial_critval(double p,int n,double a)432 static double binomial_critval (double p, int n, double a)
433 {
434 double pk, ac = 1 - a;
435 int k;
436
437 if (n <= 0 || p <= 0 || p >= 1 || a <= 0 || a >= 1) {
438 return NADBL;
439 }
440
441 for (k=n; k>0; k--) {
442 pk = binomial_cdf(p, n, k);
443 if (pk < ac) {
444 break;
445 }
446 }
447
448 return (double) (k + 1);
449 }
450
451 /**
452 * binomial_pmf:
453 * @p: success probability.
454 * @n: number of trials.
455 * @k: number of successes.
456 *
457 * Returns: the probability mass for @k successes in @n
458 * binomial trials with success probability @p.
459 */
460
binomial_pmf(double p,int n,int k)461 double binomial_pmf (double p, int n, int k)
462 {
463 double pm = binomial_cdf(p, n, k);
464
465 if (k > 0 && !na(pm)) {
466 pm -= binomial_cdf(p, n, k - 1);
467 }
468
469 return pm;
470 }
471
binomial_pmf_array(double p,int n,double * x,int T)472 static int binomial_pmf_array (double p, int n, double *x, int T)
473 {
474 double pm;
475 int i, k;
476
477 for (i=0; i<T; i++) {
478 k = x[i];
479 pm = binomial_cdf(p, n, k);
480 if (k > 0 && !na(pm)) {
481 pm -= binomial_cdf(p, n, k - 1);
482 }
483 x[i] = pm;
484 }
485
486 return 0;
487 }
488
489 /**
490 * x_factorial:
491 * @x: input value.
492 *
493 * Returns: the factorial of int(x), cast to a double, or
494 * NADBL on failure.
495 */
496
x_factorial(double x)497 double x_factorial (double x)
498 {
499 double fact;
500 int n = x;
501
502 if (x < 0) {
503 fact = NADBL;
504 } else if (x > 12) {
505 fact = cephes_gamma(1 + x);
506 if (get_cephes_errno()) {
507 fact = NADBL;
508 }
509 } else if (n == 0) {
510 fact = 1;
511 } else {
512 fact = n;
513 while (--n > 1) {
514 fact *= n;
515 }
516 }
517
518 return fact;
519 }
520
521 /**
522 * log_x_factorial:
523 * @x: input value.
524 *
525 * Returns: the log of the factorial of int(x), cast to a double, or
526 * NADBL on failure.
527 */
528
log_x_factorial(double x)529 double log_x_factorial (double x)
530 {
531 double lfact;
532 int n = x;
533
534 if (x < 0) {
535 lfact = NADBL;
536 } else if (x > 12) {
537 lfact = cephes_lgamma(1 + x);
538 if (get_cephes_errno()) {
539 lfact = NADBL;
540 }
541 } else if (n == 0) {
542 lfact = 0;
543 } else {
544 lfact = n;
545 while (--n > 1) {
546 lfact *= n;
547 }
548 lfact = log(lfact);
549 }
550
551 return lfact;
552 }
553
554 /**
555 * tcrit95:
556 * @df: degrees of freedom.
557 *
558 * Returns: the two-sided 95 percent critical value for the t
559 * distribution with @df degrees of freedom, or #NADBL on
560 * failure.
561 */
562
tcrit95(int df)563 double tcrit95 (int df)
564 {
565 double x = NADBL;
566
567 if (df > 0) {
568 x = stdtri(df, 0.975);
569 if (get_cephes_errno()) {
570 x = NADBL;
571 }
572 }
573
574 return x;
575 }
576
577 /**
578 * rhocrit95:
579 * @n: sample size.
580 *
581 * Computes the two-sided 5 percent critical value of the sample
582 * correlation coefficient for a sample of size @n. This is based
583 * on the inverse of the function which maps from the correlation
584 * coefficient, r, to a student t statistic, namely
585 *
586 * t = r / sqrt[(1—r^2) / (n-2)]
587 *
588 * The inverse is r = sqrt(t^2 / (t^2 + n - 2)).
589 *
590 * Returns: the critical value, or #NADBL on failure.
591 */
592
rhocrit95(int n)593 double rhocrit95 (int n)
594 {
595 double rc = NADBL;
596
597 if (n - 2 > 0) {
598 double tc = stdtri(n - 2, 0.975);
599
600 if (get_cephes_errno() == 0) {
601 double tc2 = tc * tc;
602
603 rc = sqrt(tc2 / (tc2 + n - 2));
604 }
605 }
606
607 return rc;
608 }
609
610 /**
611 * normal_pvalue_2:
612 * @x: double-precision value.
613 *
614 * Calculates the two-sided p-value for @x in relation to the
615 * standard normal distribution.
616 *
617 * Returns: 2 times (1 minus the value of the standard normal
618 * CDF evaluated at abs(@x)), or 0 on underflow.
619 */
620
normal_pvalue_2(double x)621 double normal_pvalue_2 (double x)
622 {
623 double p = (x < 0)? ndtr(x) : ndtr(-x);
624
625 return 2 * p;
626 }
627
628 /**
629 * normal_pvalue_1:
630 * @x: double-precision value.
631 *
632 * Calculates the one-sided p-value for @x in relation to the
633 * standard normal distribution (that is, the probability that a
634 * random variable distributed as N(0, 1) is greater than @x).
635 *
636 * Returns: 1 minus the value of the standard normal CDF
637 * evaluated at @x.
638 */
639
normal_pvalue_1(double x)640 double normal_pvalue_1 (double x)
641 {
642 double p = ndtr(x);
643
644 if (get_cephes_errno()) {
645 p = NADBL;
646 } else {
647 p = 1 - p;
648 }
649
650 return p;
651 }
652
653 /**
654 * student_cdf:
655 * @df: degrees of freedom.
656 * @x: the cutoff point in the distribution.
657 *
658 * Returns: the integral from minus infinity to @x of
659 * the Student's t distribution with @df degrees of freedom, or
660 * #NADBL on failure.
661 */
662
student_cdf(double df,double x)663 double student_cdf (double df, double x)
664 {
665 double p = NADBL;
666
667 if (df > 0) {
668 p = stdtr(df, x);
669 if (get_cephes_errno()) {
670 p = NADBL;
671 }
672 }
673
674 return p;
675 }
676
677 /**
678 * student_cdf_comp:
679 * @df: degrees of freedom.
680 * @x: the cutoff point in the distribution.
681 *
682 * Returns: the integral from @x to infinity of
683 * the t distribution with @df degrees of freedom, or
684 * #NADBL on failure.
685 */
686
student_cdf_comp(double df,double x)687 static double student_cdf_comp (double df, double x)
688 {
689 double p = NADBL;
690
691 if (df > 0) {
692 if (x > 0) {
693 p = stdtr(df, -x);
694 if (get_cephes_errno()) {
695 p = NADBL;
696 }
697 } else {
698 p = stdtr(df, x);
699 if (get_cephes_errno()) {
700 p = NADBL;
701 } else {
702 p = 1 - p;
703 }
704 }
705 }
706
707 return p;
708 }
709
710 /**
711 * normal_cdf_comp:
712 * @x: the cutoff point in the distribution.
713 *
714 * Returns: the integral from @x to infinity of the standard
715 * normal distribution, or #NADBL on failure.
716 */
717
normal_cdf_comp(double x)718 double normal_cdf_comp (double x)
719 {
720 double p;
721
722 if (x > 0) {
723 p = ndtr(-x);
724 if (get_cephes_errno()) {
725 p = NADBL;
726 }
727 } else {
728 p = ndtr(x);
729 if (get_cephes_errno()) {
730 p = NADBL;
731 } else {
732 p = 1 - p;
733 }
734 }
735
736 return p;
737 }
738
739 /**
740 * student_pvalue_1:
741 * @df: degrees of freedom.
742 * @x: the cutoff point in the distribution.
743 *
744 * Returns: the probability that t(@df) is greater than @x,
745 * or #NADBL on failure.
746 */
747
student_pvalue_1(double df,double x)748 double student_pvalue_1 (double df, double x)
749 {
750 double p = NADBL;
751
752 if (df > 0) {
753 p = stdtr(df, x);
754 if (get_cephes_errno()) {
755 p = NADBL;
756 } else {
757 p = (p >= 1.0)? 0.0 : 1 - p;
758 }
759 }
760
761 return p;
762 }
763
764 /**
765 * student_pvalue_2:
766 * @df: degrees of freedom.
767 * @x: the cutoff point in the distribution.
768 *
769 * Returns: the probability that t(@df) is greater than @x
770 * (two-sided, using the absolute value of @x), or
771 * #NADBL on failure.
772 */
773
student_pvalue_2(double df,double x)774 double student_pvalue_2 (double df, double x)
775 {
776 double p = NADBL;
777
778 if (df > 0) {
779 if (x < 0.0) {
780 p = stdtr(df, x);
781 } else {
782 p = stdtr(df, -x);
783 }
784 if (get_cephes_errno()) {
785 p = NADBL;
786 } else {
787 p *= 2;
788 }
789 }
790
791 return p;
792 }
793
794 /**
795 * student_critval:
796 * @df: degrees of freedom.
797 * @a: right-tail probability.
798 *
799 * Returns: the argument x such that the integral from x to
800 * infinity of the t(@df) density is equal to the given
801 * probability @a, or #NADBL on failure.
802 */
803
student_critval(double df,double a)804 double student_critval (double df, double a)
805 {
806 double x;
807
808 if (df < 0) {
809 return NADBL;
810 }
811
812 if (a > .10) {
813 x = stdtri(df, 1 - a);
814 } else {
815 x = -stdtri(df, a);
816 }
817
818 if (get_cephes_errno()) {
819 x = NADBL;
820 }
821
822 return x;
823 }
824
825 /**
826 * student_cdf_inverse:
827 * @df: degrees of freedom.
828 * @a: probability.
829 *
830 * Returns: the argument x such that the integral from
831 * minus infinity to @x of the t(@df) density is equal to
832 * the given probability @a, or #NADBL on failure.
833 */
834
student_cdf_inverse(double df,double a)835 double student_cdf_inverse (double df, double a)
836 {
837 double x;
838
839 if (df < 0) {
840 return NADBL;
841 }
842
843 x = stdtri(df, a);
844
845 if (get_cephes_errno()) {
846 x = NADBL;
847 }
848
849 return x;
850 }
851
852 /**
853 * chisq_cdf:
854 * @df: degrees of freedom.
855 * @x: the cutoff point in the distribution.
856 *
857 * Returns: the integral from 0 to @x of the chi-square
858 * distribution with @df degrees of freedom, or #NADBL
859 * on failure.
860 */
861
chisq_cdf(double df,double x)862 double chisq_cdf (double df, double x)
863 {
864 double p = chdtr(df, x);
865
866 /* 2016-07-22: should we return NA instead of zero
867 in case of CEPHES_UNDERFLOW? */
868 if (get_cephes_errno() == CEPHES_DOMAIN) {
869 p = NADBL;
870 }
871
872 return p;
873 }
874
875 /**
876 * chisq_cdf_comp:
877 * @df: degrees of freedom.
878 * @x: the cutoff point in the distribution.
879 *
880 * Returns: the integral from @x to infinity of the chi-square
881 * distribution with @df degrees of freedom, or #NADBL
882 * on failure.
883 */
884
chisq_cdf_comp(double df,double x)885 double chisq_cdf_comp (double df, double x)
886 {
887 double p = chdtrc(df, x);
888
889 if (get_cephes_errno() == CEPHES_DOMAIN) {
890 p = NADBL;
891 }
892
893 return p;
894 }
895
896 /**
897 * chisq_critval:
898 * @df: degrees of freedom.
899 * @a: right-tail probability.
900 *
901 * Returns: the Chi-square argument x such that the integral
902 * from x to infinity of the Chi-square density is equal
903 * to the given probability @a, or #NADBL on failure.
904 */
905
chisq_critval(int df,double a)906 static double chisq_critval (int df, double a)
907 {
908 double x = chdtri(df, a);
909
910 if (get_cephes_errno() == CEPHES_DOMAIN) {
911 x = NADBL;
912 }
913
914 return x;
915 }
916
chisq_cdf_inverse(int df,double a)917 static double chisq_cdf_inverse (int df, double a)
918 {
919 double x = NADBL;
920
921 if (df >= 1 && a >= 0.0) {
922 x = chdtri(df, 1 - a);
923 if (get_cephes_errno()) {
924 x = NADBL;
925 }
926 }
927
928 return x;
929 }
930
931 /**
932 * snedecor_cdf:
933 * @dfn: numerator degrees of freedom.
934 * @dfd: denominator degrees of freedom.
935 * @x: the cutoff point in the distribution.
936 *
937 * Returns: the integral of the F distribution with @dfn and
938 * @dfd degrees of freedom, from 0 to @x, or #NADBL on failure.
939 */
940
snedecor_cdf(double dfn,double dfd,double x)941 double snedecor_cdf (double dfn, double dfd, double x)
942 {
943 double p = NADBL;
944
945 if (dfn > 0 && dfd > 0 && x >= 0) {
946 p = fdtr(dfn, dfd, x);
947 if (get_cephes_errno()) {
948 p = NADBL;
949 }
950 }
951
952 return p;
953 }
954
955 /**
956 * snedecor_cdf_comp:
957 * @dfn: numerator degrees of freedom.
958 * @dfd: denominator degrees of freedom.
959 * @x: the cutoff point in the distribution.
960 *
961 * Returns: the integral of the F distribution with @dfn and
962 * @dfd degrees of freedom, from @x to infinity, or #NADBL
963 * on failure.
964 */
965
snedecor_cdf_comp(double dfn,double dfd,double x)966 double snedecor_cdf_comp (double dfn, double dfd, double x)
967 {
968 double p = NADBL;
969
970 if (dfn > 0 && dfd > 0 && x >= 0) {
971 p = fdtrc(dfn, dfd, x);
972 if (get_cephes_errno()) {
973 p = NADBL;
974 }
975 }
976
977 return p;
978 }
979
980 /**
981 * snedecor_critval:
982 * @dfn: numerator degrees of freedom.
983 * @dfd: denominator degrees of freedom.
984 * @a: right-tail probability.
985 *
986 * Returns: the F argument x such that the integral
987 * from x to infinity of the F density is equal
988 * to the given probability @a, or #NADBL on failure.
989 */
990
snedecor_critval(double dfn,double dfd,double a)991 double snedecor_critval (double dfn, double dfd, double a)
992 {
993 double x = NADBL;
994
995 if (dfn > 0 && dfd > 0 && a >= 0) {
996 x = fdtri(dfn, dfd, a);
997 if (get_cephes_errno()) {
998 x = NADBL;
999 }
1000 }
1001
1002 return x;
1003 }
1004
snedecor_cdf_inverse(double dfn,double dfd,double a)1005 static double snedecor_cdf_inverse (double dfn, double dfd, double a)
1006 {
1007 double x = NADBL;
1008
1009 if (dfn > 0 && dfd > 0 && a >= 0) {
1010 x = fdtri(dfn, dfd, 1 - a);
1011 if (get_cephes_errno()) {
1012 x = NADBL;
1013 }
1014 }
1015
1016 return x;
1017 }
1018
1019 /* PDFs */
1020
Binv(double p,double q)1021 static double Binv (double p, double q)
1022 {
1023 double f = NADBL;
1024
1025 errno = 0;
1026
1027 if (p > 0 && q > 0) {
1028 double x1 = lngamma(p + q);
1029 double x2 = lngamma(p);
1030 double x3 = lngamma(q);
1031
1032 if (!na(x1) && !na(x2) && !na(x3)) {
1033 f = exp(x1 - x2 - x3);
1034 if (errno) {
1035 f = NADBL;
1036 }
1037 }
1038 }
1039
1040 return f;
1041 }
1042
snedecor_pdf_array(double v1,double v2,double * x,int n)1043 static int snedecor_pdf_array (double v1, double v2, double *x, int n)
1044 {
1045 int i, err = 0;
1046
1047 errno = 0;
1048
1049 if (v1 > 0 && v2 > 0) {
1050 double xm = v1, xn = v2;
1051 double vr = xm / xn;
1052 double x1 = Binv(0.5*xm, 0.5*xn);
1053 double x2 = pow(vr, 0.5*xm);
1054 double x3, x4;
1055
1056 if (errno) {
1057 err = E_NAN;
1058 } else {
1059 for (i=0; i<n; i++) {
1060 if (!na(x[i]) && x[i] > 0) {
1061 errno = 0;
1062 x3 = pow(x[i], 0.5*xm - 1.0);
1063 x4 = pow(1.0 + vr * x[i], 0.5 * (xm+xn));
1064 x[i] = x1 * x2 * x3 / x4;
1065 if (errno) {
1066 x[i] = NADBL;
1067 }
1068 } else {
1069 x[i] = NADBL;
1070 }
1071 }
1072 }
1073 } else {
1074 err = E_DATA;
1075 }
1076
1077 if (err) {
1078 for (i=0; i<n; i++) {
1079 x[i] = NADBL;
1080 }
1081 }
1082
1083 errno = 0;
1084
1085 return err;
1086 }
1087
snedecor_pdf(double m,double n,double x)1088 static double snedecor_pdf (double m, double n, double x)
1089 {
1090 snedecor_pdf_array(m, n, &x, 1);
1091
1092 return x;
1093 }
1094
chisq_pdf_array(int m,double * x,int n)1095 static int chisq_pdf_array (int m, double *x, int n)
1096 {
1097 int i, err = 0;
1098
1099 errno = 0;
1100
1101 if (m > 0) {
1102 double m2 = m / 2.0;
1103 double x1 = pow(.5, m2);
1104 double x2 = gammafun(m2);
1105 double x3, x4;
1106
1107 if (errno) {
1108 err = E_NAN;
1109 } else {
1110 for (i=0; i<n; i++) {
1111 if (!na(x[i]) && x[i] > 0) {
1112 errno = 0;
1113 x3 = pow(x[i], m2 - 1.0);
1114 x4 = exp(-x[i] / 2.0);
1115 x[i] = (x1/x2) * x3 * x4;
1116 if (errno) {
1117 x[i] = NADBL;
1118 }
1119 } else {
1120 x[i] = NADBL;
1121 }
1122 }
1123 }
1124 } else {
1125 err = E_DATA;
1126 }
1127
1128 if (err) {
1129 for (i=0; i<n; i++) {
1130 x[i] = NADBL;
1131 }
1132 }
1133
1134 errno = 0;
1135
1136 return err;
1137 }
1138
chisq_pdf(int m,double x)1139 static double chisq_pdf (int m, double x)
1140 {
1141 chisq_pdf_array(m, &x, 1);
1142
1143 return x;
1144 }
1145
student_pdf_array(double m,double * x,int n)1146 static int student_pdf_array (double m, double *x, int n)
1147 {
1148 int i, err = 0;
1149
1150 errno = 0;
1151
1152 if (m > 0 && !na(m)) {
1153 double x1 = Binv(0.5 * m, 0.5) / sqrt(m);
1154 double x3 = 0.5 * (m + 1.0);
1155 double x2;
1156
1157 if (errno) {
1158 err = E_NAN;
1159 } else {
1160 for (i=0; i<n; i++) {
1161 if (!na(x[i])) {
1162 errno = 0;
1163 x2 = m / (m + x[i] * x[i]);
1164 x[i] = x1 * pow(x2, x3);
1165 if (errno) {
1166 x[i] = NADBL;
1167 }
1168 }
1169 }
1170 }
1171 } else {
1172 err = E_DATA;
1173 }
1174
1175 if (err) {
1176 for (i=0; i<n; i++) {
1177 x[i] = NADBL;
1178 }
1179 }
1180
1181 errno = 0;
1182
1183 return err;
1184 }
1185
student_pdf(double m,double x)1186 static double student_pdf (double m, double x)
1187 {
1188 student_pdf_array(m, &x, 1);
1189
1190 return x;
1191 }
1192
weibull_pdf_array(double k,double l,double * x,int n)1193 static int weibull_pdf_array (double k, double l,
1194 double *x, int n)
1195 {
1196 int i, err = 0;
1197
1198 errno = 0;
1199
1200 if (!na(k) && k > 0 && !na(l) && l > 0) {
1201 double x1 = k / l;
1202 double x2, x3, x4;
1203
1204 if (errno) {
1205 err = E_NAN;
1206 } else {
1207 for (i=0; i<n; i++) {
1208 if (!na(x[i]) && x[i] >= 0) {
1209 errno = 0;
1210 x3 = x[i] / l;
1211 x2 = pow(x3, k - 1.0);
1212 x3 *= x2;
1213 x4 = exp(-x3);
1214 x[i] = x1 * x2 * x4;
1215 if (errno) {
1216 x[i] = NADBL;
1217 }
1218 } else {
1219 x[i] = NADBL;
1220 }
1221 }
1222 }
1223 } else {
1224 err = E_DATA;
1225 }
1226
1227 if (err) {
1228 for (i=0; i<n; i++) {
1229 x[i] = NADBL;
1230 }
1231 }
1232
1233 errno = 0;
1234
1235 return err;
1236 }
1237
exponential_pdf_array(double mu,double * x,int n)1238 static int exponential_pdf_array (double mu, double *x, int n)
1239 {
1240 int i, err = 0;
1241
1242 if (na(mu) || mu <= 0) {
1243 err = E_INVARG;
1244 for (i=0; i<n; i++) {
1245 x[i] = NADBL;
1246 }
1247 } else {
1248 double b = 1.0 / mu;
1249
1250 for (i=0; i<n; i++) {
1251 if (x[i] < 0) {
1252 x[i] = 0;
1253 } else if (!na(x[i])) {
1254 x[i] = b * exp(-b * x[i]);
1255 } else {
1256 x[i] = NADBL;
1257 }
1258 }
1259 }
1260
1261 return err;
1262 }
1263
beta_pdf_array(double a,double b,double * x,int n)1264 static int beta_pdf_array (double a, double b, double *x, int n)
1265 {
1266 int i, err = 0;
1267
1268 if (na(a) || a <= 0 || na(b) || b <= 0) {
1269 err = E_INVARG;
1270 for (i=0; i<n; i++) {
1271 x[i] = NADBL;
1272 }
1273 } else {
1274 double k = exp(lngamma(a + b) - lngamma(a) - lngamma(b));
1275
1276 for (i=0; i<n; i++) {
1277 x[i] = k * pow(x[i], a-1) * pow(1.0 - x[i], b-1);
1278 }
1279 }
1280
1281 return err;
1282 }
1283
weibull_pdf(double k,double l,double x)1284 static double weibull_pdf (double k, double l, double x)
1285 {
1286 weibull_pdf_array(k, l, &x, 1);
1287
1288 return x;
1289 }
1290
exponential_pdf(double mu,double x)1291 static double exponential_pdf (double mu, double x)
1292 {
1293 exponential_pdf_array(mu, &x, 1);
1294
1295 return x;
1296 }
1297
beta_pdf(double a,double b,double x)1298 static double beta_pdf (double a, double b, double x)
1299 {
1300 beta_pdf_array(a, b, &x, 1);
1301
1302 return x;
1303 }
1304
1305 /**
1306 * normal_cdf:
1307 * @x: double-precision value.
1308 *
1309 * Returns: the value of the standard normal CDF evaluated
1310 * at @x, or #NADBL on failure.
1311 */
1312
normal_cdf(double x)1313 double normal_cdf (double x)
1314 {
1315 double y = ndtr(x);
1316
1317 if (get_cephes_errno()) {
1318 y = NADBL;
1319 } else if (y == 1.0) {
1320 /* do we want to do this? */
1321 y = NORM_CDF_MAX;
1322 }
1323
1324 return y;
1325 }
1326
1327 /**
1328 * normal_cdf_inverse:
1329 * @x: double-precision value.
1330 *
1331 * Returns: the argument, y, for which the area under the
1332 * Gaussian probability density function (integrated from
1333 * minus infinity to y) is equal to x, or #NADBL on failure.
1334 */
1335
normal_cdf_inverse(double x)1336 double normal_cdf_inverse (double x)
1337 {
1338 double y = ndtri(x);
1339
1340 if (get_cephes_errno()) {
1341 y = NADBL;
1342 }
1343
1344 return y;
1345 }
1346
1347 /**
1348 * normal_critval:
1349 * @a: right-tail probability.
1350 *
1351 * Returns: the argument z such that the integral from z to
1352 * infinity of the standard normal density is equal
1353 * to the given probability @a, or #NADBL on failure.
1354 */
1355
normal_critval(double a)1356 double normal_critval (double a)
1357 {
1358 double z;
1359
1360 if (a > 0.10) {
1361 z = ndtri(1.0 - a);
1362 } else {
1363 z = -ndtri(a);
1364 }
1365
1366 if (get_cephes_errno()) {
1367 z = NADBL;
1368 }
1369
1370 return z;
1371 }
1372
normal_pdf_array(double * x,int n)1373 static int normal_pdf_array (double *x, int n)
1374 {
1375 double s = 1.0 / SQRT_2_PI;
1376 int i;
1377
1378 for (i=0; i<n; i++) {
1379 if (!na(x[i])) {
1380 errno = 0;
1381 x[i] = s * exp(-0.5 * x[i] * x[i]);
1382 if (errno) {
1383 x[i] = NADBL;
1384 }
1385 }
1386 }
1387
1388 errno = 0;
1389
1390 return 0;
1391 }
1392
1393 /**
1394 * normal_pdf:
1395 * @x: double-precision value.
1396 *
1397 * Returns: the value of the standard normal PDF evaluated
1398 * at @x.
1399 */
1400
normal_pdf(double x)1401 double normal_pdf (double x)
1402 {
1403 return exp(-0.5 * x * x) / SQRT_2_PI;
1404 }
1405
1406 /**
1407 * log_normal_pdf:
1408 * @x: double-precision value.
1409 *
1410 * Returns: the value of the log-normal PDF evaluated
1411 * at @x.
1412 */
1413
log_normal_pdf(double x)1414 double log_normal_pdf (double x)
1415 {
1416 return -0.5 * x * x - LN_SQRT_2_PI;
1417 }
1418
1419 /**
1420 * gamma_cdf:
1421 * @s1: first parameter.
1422 * @s2: second parameter.
1423 * @x: reference value.
1424 * @control: see below.
1425 *
1426 * Calculates the value of the CDF of the gamma distribution
1427 * at @x. If @control equals 1, then it is assumed that the
1428 * parameters @s1 and @s2 represent the shape and scale,
1429 * respectively, otherwise it is assumed they give mean and
1430 * variance.
1431
1432 * Returns: the calculated probability, or #NADBL on failure.
1433 */
1434
gamma_cdf(double s1,double s2,double x,int control)1435 double gamma_cdf (double s1, double s2, double x, int control)
1436 {
1437 double scale, shape, p;
1438
1439 if (control == 1) {
1440 shape = s1;
1441 scale = s2;
1442 } else {
1443 scale = s2 / s1;
1444 shape = s1 / scale;
1445 }
1446
1447 /* for the cephes functions, the parameterization is
1448 inverse-scale (or "rate"), shape */
1449
1450 p = gdtr(1.0 / scale, shape, x);
1451 if (get_cephes_errno()) {
1452 p = NADBL;
1453 }
1454
1455 return p;
1456 }
1457
1458 /**
1459 * gamma_cdf_comp:
1460 * @s1: first parameter.
1461 * @s2: second parameter.
1462 * @x: reference value.
1463 * @control: see below.
1464 *
1465 * Calculates the complement of the CDF of the gamma distribution
1466 * at @x. If @control equals 1, then it is assumed that the
1467 * parameters @s1 and @s2 represent the shape and scale,
1468 * respectively, otherwise it is assumed they give mean and
1469 * variance.
1470
1471 * Returns: the calculated probability, or #NADBL on failure.
1472 */
1473
gamma_cdf_comp(double s1,double s2,double x,int control)1474 double gamma_cdf_comp (double s1, double s2, double x, int control)
1475 {
1476 double shape, scale, p;
1477
1478 if (control == 1) {
1479 shape = s1;
1480 scale = s2;
1481 } else {
1482 scale = s2 / s1; /* variance / mean */
1483 shape = s1 / scale; /* mean / scale */
1484 }
1485
1486 p = gdtrc(1.0 / scale, shape, x);
1487 if (get_cephes_errno()) {
1488 p = NADBL;
1489 }
1490
1491 return p;
1492 }
1493
1494 /**
1495 * gamma_cdf_inverse:
1496 * @shape: shape.
1497 * @scale: scale.
1498 * @p: probability.
1499 *
1500 * Returns: the argument x such that the integral from zero to @x of
1501 * the gamma density with given scale and shape parameters is equal to
1502 * the given probability @p, or #NADBL on failure. Note that the
1503 * alternate parametrization (mean, variance) is not supported.
1504 */
1505
gamma_cdf_inverse(double shape,double scale,double p)1506 double gamma_cdf_inverse (double shape, double scale, double p)
1507 {
1508 double x = NADBL;
1509
1510 if (p==0) {
1511 return 0;
1512 }
1513
1514 if (shape > 0 && scale > 0 && p > 0 && p < 1) {
1515 x = igami(shape, 1-p) * scale;
1516 if (get_cephes_errno()) {
1517 x = NADBL;
1518 }
1519 }
1520
1521 return x;
1522 }
1523
gamma_pdf_array(double shape,double scale,double * x,int n)1524 static int gamma_pdf_array (double shape, double scale,
1525 double *x, int n)
1526 {
1527 int i, err = 0;
1528
1529 errno = 0;
1530
1531 if (!na(shape) && shape > 0 && !na(scale) && scale > 0) {
1532 double x1, x2;
1533 double x3 = pow(scale, shape);
1534 double x4 = gammafun(shape);
1535
1536 if (errno || na(x4)) {
1537 err = E_NAN;
1538 } else {
1539 for (i=0; i<n; i++) {
1540 if (!na(x[i]) && x[i] > 0) {
1541 errno = 0;
1542 x1 = pow(x[i], shape - 1.0);
1543 x2 = exp(-x[i]/scale);
1544 x[i] = x1 * x2 / (x3 * x4);
1545 if (errno) {
1546 x[i] = NADBL;
1547 }
1548 } else {
1549 x[i] = NADBL;
1550 }
1551 }
1552 }
1553 } else {
1554 err = E_DATA;
1555 }
1556
1557 if (err) {
1558 for (i=0; i<n; i++) {
1559 x[i] = NADBL;
1560 }
1561 }
1562
1563 errno = 0;
1564
1565 return err;
1566 }
1567
1568 /**
1569 * gamma_pdf:
1570 * @shape: shape parameter.
1571 * @scale: scale parameter.
1572 * @x: reference value.
1573 *
1574 * Returns: the value of the gamma pdf at @x, or #NADBL on failure.
1575 */
1576
gamma_pdf(double shape,double scale,double x)1577 double gamma_pdf (double shape, double scale, double x)
1578 {
1579 gamma_pdf_array(shape, scale, &x, 1);
1580
1581 return x;
1582 }
1583
1584 /**
1585 * poisson_cdf:
1586 * @lambda: mean (also variance).
1587 * @k: test value.
1588 *
1589 * Returns: the probability of X <= @k, for X an r.v. that follows
1590 * the Poisson distribution with parameter @lambda.
1591 */
1592
poisson_cdf(double lambda,int k)1593 static double poisson_cdf (double lambda, int k)
1594 {
1595 double x = NADBL;
1596
1597 if (lambda >= 0 && k >= 0) {
1598 x = pdtr(k, lambda);
1599 if (get_cephes_errno()) {
1600 x = NADBL;
1601 }
1602 }
1603
1604 return x;
1605 }
1606
1607 /**
1608 * poisson_cdf_comp:
1609 * @lambda: mean (also variance).
1610 * @k: test value.
1611 *
1612 * Returns: the probability of X > @k, for X an r.v. that follows
1613 * the Poisson distribution with parameter @lambda.
1614 */
1615
poisson_cdf_comp(double lambda,int k)1616 static double poisson_cdf_comp (double lambda, int k)
1617 {
1618 double x = NADBL;
1619
1620 if (lambda >= 0 && k >= 0) {
1621 x = pdtrc(k, lambda);
1622 if (get_cephes_errno()) {
1623 x = NADBL;
1624 }
1625 }
1626
1627 return x;
1628 }
1629
poisson_pmf_array(double lambda,double * x,int n)1630 static int poisson_pmf_array (double lambda, double *x, int n)
1631 {
1632 double den, l0, p;
1633 int i, j, k;
1634
1635 if (lambda <= 0) {
1636 for (i=0; i<n; i++) {
1637 x[i] = NADBL;
1638 }
1639 return 0;
1640 }
1641
1642 l0 = exp(-lambda);
1643
1644 for (i=0; i<n; i++) {
1645 k = x[i];
1646 den = x_factorial((double) k);
1647 if (na(den)) {
1648 p = NADBL;
1649 } else {
1650 p = l0 * pow(lambda, (double) k) / den;
1651 }
1652 if (na(p)) {
1653 p = l0;
1654 for (j=1; j<=k; j++) {
1655 p *= lambda / j;
1656 }
1657 }
1658 x[i] = p;
1659 }
1660
1661 return 0;
1662 }
1663
1664 /**
1665 * poisson_pmf:
1666 * @lambda: mean (also variance).
1667 * @k: test value.
1668 *
1669 * Returns: the probability mass at @k, for an r.v. that follows
1670 * the Poisson distribution with parameter @lambda.
1671 */
1672
poisson_pmf(double lambda,int k)1673 double poisson_pmf (double lambda, int k)
1674 {
1675 double den, l0, p;
1676
1677 if (lambda <= 0 || k < 0) {
1678 return NADBL;
1679 }
1680
1681 den = x_factorial((double) k);
1682 l0 = exp(-lambda);
1683
1684 if (na(den)) {
1685 p = NADBL;
1686 } else {
1687 p = l0 * pow(lambda, (double) k) / den;
1688 }
1689
1690 if (na(p)) {
1691 int i;
1692
1693 p = l0;
1694 for (i=1; i<=k; i++) {
1695 p *= lambda / i;
1696 }
1697 }
1698
1699 return p;
1700 }
1701
1702 /* The following is probably horribly inefficient --
1703 investigate the possibility of using Chebyshev,
1704 or perhaps binary search.
1705 */
1706
poisson_critval(double mu,double a)1707 static double poisson_critval (double mu, double a)
1708 {
1709 double pk = 0.0;
1710 double ac = 1 - a;
1711 int k, k0 = 0;
1712
1713 if (mu <= 0 || a <= 0 || a >= 1) {
1714 return NADBL;
1715 }
1716
1717 if (mu >= 10 && a < 0.5) {
1718 k0 = mu - 1;
1719 pk = poisson_cdf(mu, k0++);
1720 }
1721
1722 for (k=k0; ; k++) {
1723 pk = poisson_cdf(mu, k);
1724 if (pk >= ac) {
1725 break;
1726 }
1727 }
1728
1729 return (double) k;
1730 }
1731
1732 /**
1733 * poisson_cdf_inverse:
1734 * @lambda: Poisson parameter (mean = variance).
1735 * @p: cumulative probability.
1736 *
1737 * Returns: the Poisson variable @x such that the integral
1738 * from 0 to @x of the Poisson density is equal to the
1739 * given probability @p.
1740 */
1741
1742 #include "poisson.h" /* use Mike Giles's code */
1743
poisson_cdf_inverse(double lambda,double p)1744 static double poisson_cdf_inverse (double lambda, double p)
1745 {
1746 return poissinv(p, lambda);
1747 }
1748
1749 #if 0
1750
1751 /* Returns the Poisson parameter lambda such that the
1752 Poisson CDF evaluated at @k equals @p.
1753 */
1754
1755 static double gretl_pdtri (int k, double p)
1756 {
1757 double lambda = NADBL;
1758
1759 if (k >= 0 && p >= 0 && p <= 1) {
1760 lambda = pdtri(k, p);
1761 if (get_cephes_errno()) {
1762 lambda = NADBL;
1763 }
1764 }
1765
1766 return lambda;
1767 }
1768
1769 #endif
1770
weibull_critval(double shape,double scale,double rtail)1771 static double weibull_critval (double shape, double scale,
1772 double rtail)
1773 {
1774 double ret = NADBL;
1775
1776 if (shape > 0 && scale > 0 && rtail > 0 && rtail < 1) {
1777 ret = scale * pow(-log(rtail), 1.0 / shape);
1778 }
1779
1780 return ret;
1781 }
1782
1783 /**
1784 * weibull_cdf:
1785 * @shape: shape parameter > 0.
1786 * @scale: scale parameter > 0.
1787 * @x: test value.
1788 *
1789 * Returns: the probability of X <= @x, for X an r.v. that follows
1790 * the Weibull distribution with parameters @shape and @scale.
1791 */
1792
weibull_cdf(double shape,double scale,double x)1793 double weibull_cdf (double shape, double scale, double x)
1794 {
1795 double ret = NADBL;
1796
1797 if (shape > 0 && scale > 0 && !na(x)) {
1798 if (x == 0.0) {
1799 ret = 0.0;
1800 } else if (x > 0.0) {
1801 ret = 1.0 - exp(-pow(x/scale, shape));
1802 }
1803 }
1804
1805 return ret;
1806 }
1807
1808 /**
1809 * exponential_cdf:
1810 * @mu: scale parameter > 0.
1811 * @x: test value.
1812 *
1813 * Returns: the probability of X <= @x, for X an r.v. that follows
1814 * the exponential distribution with scale parameter @mu.
1815 */
1816
exponential_cdf(double mu,double x)1817 double exponential_cdf (double mu, double x)
1818 {
1819 double ret = NADBL;
1820
1821 if (mu > 0 && !na(x)) {
1822 if (x < 0.0) {
1823 ret = 0.0;
1824 } else {
1825 ret = 1.0 - exp(-x / mu);
1826 }
1827 }
1828
1829 return ret;
1830 }
1831
weibull_cdf_comp(double shape,double scale,double x)1832 static double weibull_cdf_comp (double shape, double scale, double x)
1833 {
1834 double ret = NADBL;
1835
1836 if (shape > 0 && scale > 0 && !na(x)) {
1837 if (x == 0.0) {
1838 ret = 1.0;
1839 } else if (x > 0.0) {
1840 ret = exp(-pow(x/scale, shape));
1841 }
1842 }
1843
1844 return ret;
1845 }
1846
1847 /**
1848 * GED_pdf:
1849 * @nu: shape parameter.
1850 * @x: reference value.
1851 *
1852 * Returns: the density function of the Generalized Error distribution
1853 * with shape parameter @nu at @x, or #NADBL on failure.
1854 */
1855
GED_pdf(double nu,double x)1856 double GED_pdf (double nu, double x)
1857 {
1858 if (nu > 0) {
1859 double lg1 = lngamma(1/nu);
1860 double lg3 = lngamma(3/nu);
1861 double lC = 0.5*(lg3 - 3*lg1);
1862 double k = pow(0.5, 1/nu) * exp(0.5*(lg1 - lg3));
1863 double znu = pow(fabs(x/k), nu);
1864
1865 return (0.5 * nu) * exp(lC - 0.5 * znu);
1866 } else {
1867 return NADBL;
1868 }
1869 }
1870
GED_pdf_array(double nu,double * x,int n)1871 static int GED_pdf_array (double nu, double *x, int n)
1872 {
1873 int i, err = 0;
1874
1875 if (nu > 0) {
1876 double lg1 = lngamma(1/nu);
1877 double lg3 = lngamma(3/nu);
1878 double lC = 0.5*(lg3 - 3*lg1);
1879 double k = pow(0.5, 1/nu) * exp(0.5*(lg1 - lg3));
1880 double znu;
1881
1882 for (i=0; i<n; i++) {
1883 if (!na(x[i])) {
1884 znu = pow(fabs(x[i]/k), nu);
1885 x[i] = (0.5 * nu) * exp(lC - 0.5 * znu);
1886 } else {
1887 x[i] = NADBL;
1888 }
1889 }
1890 } else {
1891 err = E_DATA;
1892 }
1893
1894 if (err) {
1895 for (i=0; i<n; i++) {
1896 x[i] = NADBL;
1897 }
1898 }
1899
1900 return err;
1901 }
1902
1903 /**
1904 * GED_cdf:
1905 * @nu: shape parameter.
1906 * @x: reference value.
1907 *
1908 * Calculates the value of the CDF of the Generalized Error
1909 * distribution with parameter @nu at @x. We exploit the fact that
1910 * if x ~ GED(n), then |x/k|^n is a Gamma rv.
1911 *
1912 * Returns: the calculated probability, or #NADBL on failure.
1913 */
1914
GED_cdf(double nu,double x)1915 double GED_cdf (double nu, double x)
1916 {
1917 if (nu > 0) {
1918 int sgn = (x > 0)? 1 : -1;
1919 double p = 1/nu;
1920 double lg1 = lngamma(p);
1921 double lg3 = lngamma(3*p);
1922 double k = pow(0.5, p) * exp(0.5*(lg1 - lg3));
1923 double znu = pow(fabs(x/k), nu);
1924 double P = gamma_cdf(p, 2, znu, 1);
1925
1926 P = 0.5 * (1 + sgn*P);
1927 return P;
1928 } else {
1929 return NADBL;
1930 }
1931 }
1932
1933 /**
1934 * GED_cdf_comp:
1935 * @nu: shape parameter.
1936 * @x: reference value.
1937 *
1938 * Calculates the complement of the CDF of the Generalized Error
1939 * distribution with parameter @nu at @x. We exploit the fact that
1940 * if x ~ GED(n), then |x/k|^n is a Gamma rv.
1941 *
1942 * Returns: the calculated probability, or #NADBL on failure.
1943 */
1944
GED_cdf_comp(double nu,double x)1945 double GED_cdf_comp (double nu, double x)
1946 {
1947 if (nu > 0) {
1948 int sgn = (x > 0)? 1 : -1;
1949 double p = 1/nu;
1950 double lg1 = lngamma(p);
1951 double lg3 = lngamma(3*p);
1952 double k = pow(0.5, p) * exp(0.5*(lg1 - lg3));
1953 double znu = pow(fabs(x/k), nu);
1954 double P = gamma_cdf_comp(p, 2, znu, 1);
1955
1956 P = (sgn == 1) ? 0.5 * P : 1 - 0.5 * P;
1957 return P;
1958 } else {
1959 return NADBL;
1960 }
1961 }
1962
1963 /**
1964 * GED_cdf_inverse:
1965 * @nu: shape parameter.
1966 * @a: probability.
1967 *
1968 * Returns: the argument x such that the integral from minus infinity
1969 * to @x of the standardized GED density with shape parameter @nu is
1970 * equal to the given probability @a, or #NADBL on failure. We exploit
1971 * the well-known relationship between the standardized GED and the
1972 * Gamma variates.
1973 */
1974
GED_cdf_inverse(double nu,double a)1975 double GED_cdf_inverse (double nu, double a)
1976 {
1977 if (nu > 0 && a < 1 && a > 0) {
1978 double a2, p, lg1, lg3, sd, x;
1979 int sgn;
1980
1981 if (a > 0.5) {
1982 a2 = 2*a - 1;
1983 sgn = 1;
1984 } else {
1985 a2 = 1 - 2*a;
1986 sgn = -1;
1987 }
1988
1989 p = 1/nu;
1990 lg1 = lngamma(p);
1991 lg3 = lngamma(3*p);
1992 sd = pow(2.0, p) * exp(0.5*(lg3 - lg1));
1993 x = gamma_cdf_inverse(p, 2.0, a2);
1994
1995 return sgn * pow(x, p) / sd;
1996 } else {
1997 return NADBL;
1998 }
1999 }
2000
2001 /**
2002 * laplace_pdf:
2003 * @mu: mean.
2004 * @b: scale (greater than 0).
2005 * @x: reference value.
2006 *
2007 * Returns: the density function of the Laplace distribution
2008 * with mean @mu and scale @b evaluated at @x, or #NADBL on failure.
2009 */
2010
laplace_pdf(double mu,double b,double x)2011 double laplace_pdf (double mu, double b, double x)
2012 {
2013 if (b > 0) {
2014 return exp(-fabs(x - mu)/b) / (2*b);
2015 } else {
2016 return NADBL;
2017 }
2018 }
2019
laplace_pdf_array(double mu,double b,double * x,int n)2020 static int laplace_pdf_array (double mu, double b,
2021 double *x, int n)
2022 {
2023 int i, err = 0;
2024
2025 if (b > 0) {
2026 for (i=0; i<n; i++) {
2027 if (!na(x[i])) {
2028 x[i] = exp(-fabs(x[i] - mu)/b) / (2*b);
2029 } else {
2030 x[i] = NADBL;
2031 }
2032 }
2033 } else {
2034 err = E_DATA;
2035 for (i=0; i<n; i++) {
2036 x[i] = NADBL;
2037 }
2038 }
2039
2040 return err;
2041 }
2042
2043 /**
2044 * laplace_cdf:
2045 * @mu: mean.
2046 * @b: scale (greater than 0).
2047 * @x: reference value.
2048 *
2049 * Returns: the CDF of the Laplace distribution
2050 * with mean @mu and scale @b evaluated at @x, or #NADBL on failure.
2051 */
2052
laplace_cdf(double mu,double b,double x)2053 double laplace_cdf (double mu, double b, double x)
2054 {
2055 if (b > 0) {
2056 if (x < mu) {
2057 return 0.5 * exp((x-mu)/b);
2058 } else {
2059 return 1 - 0.5 * exp(-(x-mu)/b);
2060 }
2061 } else {
2062 return NADBL;
2063 }
2064 }
2065
2066 /**
2067 * laplace_cdf_comp:
2068 * @mu: mean.
2069 * @b: scale (greater than 0).
2070 * @x: reference value.
2071 *
2072 * Returns: the complement of the CDF of the Laplace distribution
2073 * with mean @mu and scale @b evaluated at @x, or #NADBL on failure.
2074 */
2075
laplace_cdf_comp(double mu,double b,double x)2076 double laplace_cdf_comp (double mu, double b, double x)
2077 {
2078 if (b > 0) {
2079 if (x < mu) {
2080 return 1 - 0.5 * exp((x-mu)/b);
2081 } else {
2082 return 0.5 * exp(-(x-mu)/b);
2083 }
2084 } else {
2085 return NADBL;
2086 }
2087 }
2088
2089 /**
2090 * laplace_cdf_inverse:
2091 * @mu: mean.
2092 * @b: scale (greater than 0).
2093 * @a: probability.
2094 *
2095 * Returns: the argument x such that the integral from minus infinity
2096 * to @x of the Laplace density with mean @mu and scale @b is
2097 * equal to the given probability @a, or #NADBL on failure.
2098 */
2099
laplace_cdf_inverse(double mu,double b,double a)2100 double laplace_cdf_inverse (double mu, double b, double a)
2101 {
2102 if (b <= 0 || a < 0 || a > 1) {
2103 return NADBL;
2104 } else if (a == 0.0) {
2105 return -1.0 / 0.0;
2106 } else if (a == 1.0) {
2107 return 1.0 / 0.0;
2108 } else {
2109 int sgn = a - 0.5 < 0 ? -1 : 1;
2110
2111 return mu - b*sgn * log(1.0 - 2*fabs(a - 0.5));
2112 }
2113 }
2114
2115 /**
2116 * johansen_trace_pval:
2117 * @N: the number of potentially cointegrated variables
2118 * minus the cointegration rank under H0.
2119 * @det: index of deterministic setup of model (0 to 4).
2120 * @T: the sample size, or 0 for an asymptotic result.
2121 * @tr: the trace statistic.
2122 *
2123 * Returns: the p-value of the trace statistic, computed
2124 * via Doornik's gamma approximation, or #NADBL on failure.
2125 */
2126
johansen_trace_pval(int N,int det,int T,double tr)2127 double johansen_trace_pval (int N, int det, int T, double tr)
2128 {
2129 double (*pvfunc) (double, int, int, int);
2130 double pv = NADBL;
2131
2132 pvfunc = get_plugin_function("trace_pvalue");
2133
2134 if (pvfunc != NULL) {
2135 pv = (*pvfunc) (tr, N, det, T);
2136 }
2137
2138 return pv;
2139 }
2140
2141 /* below: non-central distributions: chi-square, F and Student's t
2142 */
2143
2144 #define qsmall(sum,x) (sum < 1.0e-30 || (x) < 1.0e-8*sum)
2145
2146 /**
2147 * nc_chisq_cdf:
2148 * @df: degrees of freedom.
2149 * @delta: noncentrality parameter.
2150 * @x: reference value.
2151 *
2152 * Calculates the value at @x of the CDF of the noncentral chi^2
2153 * distribution with @df dof and noncentrality parameter equal to
2154 * @delta.
2155 *
2156 * This is a version of cumchn() from dcdflib, de-spaghettized by
2157 * Jack Lucchetti (2015-06-21). The original algorithm uses formula
2158 * 26.4.25 from Abramowitz and Stegun, Handbook of Mathematical
2159 * Functions, US NBS (1966).
2160 *
2161 * Returns: the calculated probability, or #NADBL on failure.
2162 */
2163
nc_chisq_cdf(double df,double delta,double x)2164 double nc_chisq_cdf (double df, double delta, double x)
2165 {
2166 double adj, centaj, centwt, chid2, dfd2, lcntaj, lcntwt;
2167 double lfact, pcent, pterm, sum, sumadj, term, wt, xnonc;
2168 double T1, T2, T3;
2169 int i, icent, iterb, iterf;
2170 int itermax = 1000;
2171
2172 if (x < 0.0) {
2173 return 1.0;
2174 }
2175
2176 if (df <= 0.0) {
2177 return NADBL;
2178 }
2179
2180 if (delta <= 1.0e-10) {
2181 /*
2182 When non-centrality parameter is (essentially) zero,
2183 use ordinary chi-square distribution
2184 */
2185 return chisq_cdf(df, x);
2186 }
2187
2188 xnonc = delta / 2.0;
2189
2190 /*
2191 The following code calculates the weight, chi-square, and
2192 adjustment term for the central term in the infinite series.
2193 The central term is the one in which the poisson weight is
2194 greatest. The adjustment term is the amount that must
2195 be subtracted from the chi-square to move up two degrees
2196 of freedom.
2197 */
2198 icent = (xnonc < 1.0) ? 1 : (int) trunc(xnonc);
2199 chid2 = x / 2.0;
2200
2201 /*
2202 Calculate central weight term
2203 */
2204
2205 T1 = (double) (icent + 1);
2206 lfact = lngamma(T1);
2207 lcntwt = -xnonc + (double)icent*log(xnonc) - lfact;
2208 centwt = exp(lcntwt);
2209
2210 /*
2211 Calculate central chi-square
2212 */
2213 T2 = df + 2.0 * (double) icent;
2214 pcent = chisq_cdf(T2, x);
2215
2216 /*
2217 Calculate central adjustment term
2218 */
2219
2220 dfd2 = df / 2.0 + icent;
2221 T3 = dfd2 + 1.0;
2222 lfact = lngamma(T3);
2223 lcntaj = dfd2 * log(chid2) - chid2 - lfact;
2224 centaj = exp(lcntaj);
2225 sum = centwt * pcent;
2226
2227 /*
2228 Sum backwards from the central term towards zero.
2229 Quit whenever either
2230 (1) the zero term is reached, or
2231 (2) the term gets small relative to the sum, or
2232 (3) More than NTIRED terms are totaled.
2233 */
2234
2235 iterb = 0;
2236 sumadj = 0.0;
2237 adj = centaj;
2238 wt = centwt;
2239 i = icent;
2240
2241 do {
2242 dfd2 = df/2.0 + i;
2243 /*
2244 Adjust chi-square for two fewer degrees of freedom.
2245 The adjusted value ends up in PTERM.
2246 */
2247 adj = adj * dfd2/chid2;
2248 sumadj += adj;
2249 pterm = pcent + sumadj;
2250 /*
2251 Adjust poisson weight for J decreased by one
2252 */
2253 wt *= ((double)i/xnonc);
2254 term = wt * pterm;
2255 sum += term;
2256 i--;
2257 iterb++;
2258 } while (iterb <= itermax && !qsmall(sum, term) && i > 0);
2259
2260 /*
2261 Now sum forward from the central term towards infinity.
2262 Quit when either
2263 (1) the term gets small relative to the sum, or
2264 (2) More than NTIRED terms are totaled.
2265 */
2266 iterf = 0;
2267 sumadj = adj = centaj;
2268 wt = centwt;
2269 i = icent;
2270
2271 do {
2272 /*
2273 Update weights for next higher J
2274 */
2275 wt *= (xnonc/(double)(i+1));
2276 /*
2277 Calculate PTERM and add term to sum
2278 */
2279 pterm = pcent - sumadj;
2280 term = wt * pterm;
2281 sum += term;
2282 /*
2283 Update adjustment term for DF for next iteration
2284 */
2285 i++;
2286 dfd2 = df/2.0 + i;
2287 adj = adj * chid2/dfd2;
2288 sumadj += adj;
2289 iterf++;
2290 } while (iterf <= itermax && !qsmall(sum, term));
2291
2292 return sum;
2293 }
2294
2295 /**
2296 * nc_chisq_pdf_array:
2297 * @p: degrees of freedom.
2298 * @c: noncentrality parameter.
2299 * @x: array of arguments (overwritten on exit).
2300 * @n: no. of elements in x.
2301 *
2302 * Calculates the value at @x of the CDF of the noncentral chi^2
2303 * distribution with @p dof and noncentrality parameter equal
2304 * to @c.
2305 *
2306 * Returns: an error code, as appropriate.
2307 */
2308
nc_chisq_pdf_array(double p,double c,double * x,int n)2309 static int nc_chisq_pdf_array (double p, double c, double *x, int n)
2310 {
2311 int i, err = 0;
2312 double k, a, b;
2313
2314 if (fabs(c) < 1.0e-10) {
2315 return chisq_pdf_array((int) floor(p), x, n);
2316 }
2317
2318 if (p <= 0 || c < 0) {
2319 return E_DATA;
2320 }
2321
2322 k = p/2.0 - 1;
2323
2324 for (i=0; i<n; i++) {
2325 if (na(x[i]) || x[i] < 0) {
2326 x[i] = NADBL;
2327 } else {
2328 a = exp(-0.5*(x[i]+c) + k/2.0 * log(x[i]/c)) / 2.0;
2329 b = gretl_bessel('I', k, sqrt(x[i]*c), &err);
2330 if (err) {
2331 break;
2332 }
2333 x[i] = a*b;
2334 }
2335 }
2336
2337 return err;
2338 }
2339
nc_chisq_pdf(double p,double c,double x)2340 double nc_chisq_pdf (double p, double c, double x)
2341 {
2342 nc_chisq_pdf_array(p, c, &x, 1);
2343 return x;
2344 }
2345
2346 /**
2347 * nc_chisq_cdf_inverse:
2348 * @p: degrees of freedom
2349 * @c: noncentrality parameter.
2350 * @q: probability.
2351 *
2352 * Calculates the @q-th quantile of the noncentral chi^2
2353 * distribution with @n1, @n2 dof and noncentrality parameter equal to
2354 * @c via a rough and not particularly clever root-finding
2355 * algorithm. Maybe this can be more efficient by using logs. Some
2356 * experimentation needed.
2357 *
2358 * Returns: the calculated quantile, or #NADBL on failure.
2359 */
2360
nc_chisq_cdf_inverse(double p,double c,double q)2361 static double nc_chisq_cdf_inverse (double p, double c, double q)
2362 {
2363 double x, d0, d1;
2364 int iter, subiter, retry;
2365 double F, f, dir;
2366
2367 if (p < 0 || c < 0 || q <= 0 || q >= 1) {
2368 return NADBL;
2369 }
2370
2371 if (fabs(c) < 1.0e-10) {
2372 /* don't bother for infinitesimal c */
2373 return chisq_cdf_inverse(p, q);
2374 }
2375
2376 /* start from the mean (safe bet) */
2377 x = p + c;
2378 d0 = 1.0e7;
2379 iter = 0;
2380
2381 while (fabs(d0) > 1.0e-10 && iter < 1000) {
2382 F = nc_chisq_cdf(p, c, x);
2383 f = nc_chisq_pdf(p, c, x);
2384 d0 = F - q;
2385 dir = d0/f;
2386 d1 = 1.0e7;
2387 retry = 1;
2388 subiter = 0;
2389
2390 while (retry && subiter < 100) {
2391 if ((x-dir) > 0) {
2392 d1 = F - nc_chisq_cdf(p, c, x - dir);
2393 }
2394 dir /= 2.0;
2395 retry = (x-dir) < 0 || fabs(d1) > fabs(d0);
2396 subiter++;
2397 }
2398
2399 if (subiter >= 100) {
2400 x = NADBL;
2401 break;
2402 } else {
2403 x -= dir*2;
2404 d0 = d1;
2405 iter++;
2406 }
2407 }
2408
2409 if (iter >= 1000) {
2410 x = NADBL;
2411 }
2412
2413 return x;
2414 }
2415
2416 /**
2417 * nc_snedecor_cdf:
2418 * @dfn: degrees of freedom (num).
2419 * @dfd: degrees of freedom (den).
2420 * @delta: noncentrality parameter.
2421 * @x: reference value.
2422 *
2423 * Calculates the value at @x of the CDF of the noncentral F
2424 * distribution with @dfn, @dfd dof and noncentrality parameter equal
2425 * to @delta.
2426 *
2427 * This is a version of cumfnc() from dcdflib, de-spaghettized by
2428 * Jack Lucchetti (2015-06-21). The original algorithm uses formula
2429 * 26.6.18 from Abramowitz and Stegun, Handbook of Mathematical
2430 * Functions, US NBS (1966).
2431 *
2432 * Returns: the calculated probability, or #NADBL on failure.
2433 */
2434
2435
nc_snedecor_cdf(double dfn,double dfd,double delta,double x)2436 double nc_snedecor_cdf (double dfn, double dfd, double delta, double x)
2437 {
2438 double dsum, prod, xx, yy, adn, aup, b, betdn, betup;
2439 double centwt, dnterm, sum, upterm, xmult, xnonc;
2440 double T1, T2, T3, T4, T5, T6;
2441 int i, icent;
2442
2443 if (x < 0.0) {
2444 return 1.0;
2445 }
2446
2447 if (dfn <= 0.0 || dfd <= 0.0) {
2448 return NADBL;
2449 }
2450
2451 if (delta <= 1.0e-10) {
2452 /*
2453 When non-centrality parameter is (essentially) zero, use
2454 ordinary F distribution
2455 */
2456 return snedecor_cdf(dfn, dfd, x);
2457 } else {
2458 xnonc = delta / 2.0;
2459 }
2460
2461 /*
2462 Calculate the central term of the poisson weighting factor.
2463 */
2464 icent = (xnonc < 1.0) ? 1 : (int) trunc(xnonc);
2465
2466 /*
2467 Compute central weight term
2468 */
2469 T1 = (double) (icent + 1);
2470 centwt = exp(-xnonc + (double) icent * log(xnonc) - lngamma(T1));
2471
2472 /*
2473 Compute central incomplete beta term
2474 Assure that minimum of arg to beta and 1 - arg is computed
2475 accurately.
2476 */
2477 prod = dfn * x;
2478 dsum = dfd + prod;
2479 yy = dfd / dsum;
2480 if (yy > 0.5) {
2481 xx = prod / dsum;
2482 yy = 1.0 - xx;
2483 } else {
2484 xx = 1.0 - yy;
2485 }
2486
2487 T2 = dfn / 2.0 + (double) icent;
2488 T3 = dfd / 2.0;
2489 betdn = incbet(T2, T3, xx);
2490 adn = dfn / 2.0 + (double) icent;
2491 aup = adn;
2492 b = dfd / 2.0;
2493 betup = betdn;
2494 sum = centwt * betdn;
2495
2496 /*
2497 Now sum terms backward from icent until convergence or all done
2498 */
2499
2500 xmult = centwt;
2501 i = icent;
2502 T4 = adn + b;
2503 T5 = adn + 1.0;
2504 dnterm = exp(lngamma(T4)-lngamma(T5)-lngamma(b) + adn*log(xx) + b*log(yy));
2505
2506 while (!qsmall(sum, xmult*betdn) && i > 0) {
2507 xmult *= (double) i /xnonc;
2508 i--;
2509 adn -= 1.0;
2510 dnterm *= (adn + 1.0)/((adn+b)*xx);
2511 betdn += dnterm;
2512 sum += xmult * betdn;
2513 }
2514
2515 i = icent+1;
2516
2517 /*
2518 Now sum forwards until convergence
2519 */
2520 xmult = centwt;
2521 if (aup-1.0+b == 0) {
2522 upterm = exp(-lngamma(aup) - lngamma(b) +
2523 (aup-1.0)*log(xx) + b*log(yy));
2524 } else {
2525 T6 = aup - 1.0 + b;
2526 upterm = exp(lngamma(T6) - lngamma(aup) - lngamma(b) +
2527 (aup-1.0)*log(xx) + b*log(yy));
2528 }
2529
2530 while (!qsmall(sum, xmult*betup)) {
2531 xmult *= (xnonc/(double)i);
2532 i++;
2533 aup += 1.0;
2534 upterm *= (aup + b - 2.0) * xx/(aup-1.0);
2535 betup -= upterm;
2536 sum += (xmult*betup);
2537 }
2538
2539 return sum;
2540 }
2541 #undef qsmall
2542
2543 /**
2544 * nc_snedecor_pdf_array:
2545 * @dfn: degrees of freedom (num).
2546 * @dfd: degrees of freedom (den).
2547 * @c: noncentrality parameter.
2548 * @x: array of arguments (overwritten on exit).
2549 * @n: no. of elements in x.
2550 *
2551 * Calculates the value at @x of the CDF of the noncentral F
2552 * distribution with @dfn, @dfd dof and noncentrality parameter equal
2553 * to @c.
2554 *
2555 * Source: S. Kay, Fundamentals of Statistical Signal Processing:
2556 * Detection Theory, (New Jersey: Prentice Hall, 1998),
2557 * <TeX>
2558 * p(x) = \sum\limits_{k=0}^\infty
2559 * \frac{e^{-c/2}(c/2)^k}{k!} % Poisson weights
2560 * \frac{1}{B\left(\frac{\nu_2}{2},\frac{\nu_1}{2}+k\right)} % Beta function
2561 * \left(\frac{\nu_1}{\nu_2}\right)^{\frac{\nu_1}{2}+k}
2562 * \left(\frac{\nu_2}{\nu_2+\nu_1f}\right)^{\frac{\nu_1+\nu_2}{2}+k}
2563 * x^{\nu_1/2-1+k}
2564 * </TeX>
2565 * coded in C by Jack Lucchetti (2015-06-27).
2566 *
2567 * Returns: an error code, as appropriate.
2568 */
2569
ncf_pdf_array(double dfn,double dfd,double c,double * x,int n)2570 static int ncf_pdf_array (double dfn, double dfd, double c,
2571 double *x, int n)
2572 {
2573 double ch, k1, k2, k;
2574 double a, b, l, pw, beta;
2575 double pwi, betai;
2576 double errtol = 1.0e-16;
2577 double maxit = 256;
2578 double *vx, *vz;
2579 int i, t, start, iter;
2580 int err = 0;
2581
2582 if (dfd <= 0.0 || dfn <= 0.0 || c < 0.0) {
2583 return E_DATA;
2584 }
2585
2586 if (fabs(c) <= 1.0e-10) {
2587 /*
2588 When non-centrality parameter is (essentially) zero, use
2589 ordinary F distribution
2590 */
2591 return snedecor_pdf_array(dfn, dfd, x, n);
2592 }
2593
2594 vx = malloc(n * sizeof *vx);
2595 vz = malloc(n * sizeof *vz);
2596
2597 if (vx == NULL || vz == NULL) {
2598 err = E_ALLOC;
2599 goto bailout;
2600 }
2601
2602 /* fill up auxiliary vectors */
2603 ch = c/2.0;
2604 k1 = dfn/2;
2605 k2 = (dfn + dfd)/2;
2606 k = log(dfn) - log(dfd);
2607
2608 for(t=0; t<n; t++) {
2609 if (na(x[t]) || x[t] < 0) {
2610 vx[t] = vz[t] = NADBL;
2611 } else {
2612 vx[t] = log(x[t]);
2613 vz[t] = log(dfd) - log(dfd + dfn * x[t]);
2614 }
2615 }
2616
2617 /* start from central Poisson weight */
2618
2619 start = (int) floor(ch);
2620
2621 /* Poisson weight */
2622 pw = exp(-ch - lngamma(start+1) + start * log(ch));
2623 /* Beta */
2624 beta = exp(lngamma(k1+start) + lngamma(dfd/2.0) - lngamma(start + k2));
2625 a = pw / beta;
2626 b = (k1 + start) * k;
2627
2628 for (t=0; t<n; t++) {
2629 if (na(x[t]) || x[t] < 0) {
2630 x[t] = NADBL;
2631 } else {
2632 l = b + (start + k1 - 1) * vx[t] + (start + k2) * vz[t];
2633 x[t] = a * exp(l);
2634 }
2635 }
2636
2637 /*
2638 First, go back from start to 0
2639 */
2640
2641 pwi = pw;
2642 betai = beta;
2643 iter = 0;
2644
2645 for (i = start-1; i>=0 && pwi>errtol && iter < maxit; i--) {
2646 iter++;
2647 pwi *= (i + 1.0)/ch;
2648 betai *= (k2 + i)/(k1 + i);
2649 a = pwi / betai;
2650 b = (k1 + i) * k;
2651 for (t=0; t<n; t++) {
2652 if (!na(x[t])) {
2653 l = b + (i + k1 - 1) * vx[t] + (i + k2) * vz[t];
2654 x[t] += a * exp(l);
2655 }
2656 }
2657 }
2658
2659 /*
2660 Then, go from start all the way up as necessary
2661 */
2662
2663 iter = 0;
2664 pwi = pw;
2665 betai = beta;
2666
2667 for (i = start+1; pwi>errtol && iter<maxit; i++) {
2668 iter++;
2669 pwi *= ch/i;
2670 betai *= (k1 + i - 1.0)/(k2 + i - 1.0);
2671 a = pwi / betai;
2672 b = (k1 + i) * k;
2673 for (t=0; t<n; t++) {
2674 if (!na(x[t])) {
2675 l = b + (i + k1 - 1) * vx[t] + (i + k2) * vz[t];
2676 x[t] += a * exp(l);
2677 }
2678 }
2679 }
2680
2681 bailout:
2682
2683 free(vx);
2684 free(vz);
2685
2686 return err;
2687 }
2688
ncf_pdf(double dfn,double dfd,double c,double x)2689 double ncf_pdf (double dfn, double dfd, double c, double x)
2690 {
2691 ncf_pdf_array(dfn, dfd, c, &x, 1);
2692
2693 return x;
2694 }
2695
2696 /**
2697 * ncf_cdf_inverse:
2698 * @n1: degrees of freedom (numerator).
2699 * @n2: degrees of freedom (denominator).
2700 * @c: noncentrality parameter.
2701 * @q: probability.
2702 *
2703 * Calculates the @q-th quantile of the noncentral F distribution with
2704 * @n1, @n2 dof and noncentrality parameter equal to @c via a rough
2705 * and not particularly clever root-finding algorithm. Maybe this can
2706 * be more efficient by using logs. Some experimentation needed.
2707 *
2708 * Returns: the calculated quantile, or #NADBL on failure.
2709 */
2710
ncf_cdf_inverse(double n1,double n2,double c,double q)2711 static double ncf_cdf_inverse (double n1, double n2, double c, double q)
2712 {
2713 double x, d0, d1;
2714 int iter, subiter;
2715 double F, f, dir;
2716
2717 if (n2 < 1 || n1 < 1 || c < 0 || q <= 0 || q >= 1) {
2718 return NADBL;
2719 }
2720
2721 x = 0.5;
2722 d0 = 1.0e7;
2723 iter = 0;
2724
2725 while (fabs(d0) > 1.0e-10 && iter < 1000) {
2726 F = nc_snedecor_cdf(n1, n2, c, x);
2727 f = ncf_pdf(n1, n2, c, x);
2728 d0 = F - q;
2729 dir = d0/f;
2730 d1 = 1.0e7;
2731 subiter = 0;
2732
2733 while (fabs(d1) > fabs(d0) && subiter < 100) {
2734 d1 = F - nc_snedecor_cdf(n1, n2, c, x - dir);
2735 dir /= 2.0;
2736 subiter++;
2737 }
2738
2739 if (subiter >= 100) {
2740 x = NADBL;
2741 break;
2742 } else {
2743 x -= dir*2;
2744 d0 = d1;
2745 iter++;
2746 }
2747 }
2748
2749 if (iter >= 1000) {
2750 x = NADBL;
2751 }
2752
2753 return x;
2754 }
2755
2756 #ifndef ISQRT_2
2757 #define ISQRT_2 .707106781186547524401
2758 #endif
2759
2760 /**
2761 * nc_student_cdf:
2762 * @df: degrees of freedom.
2763 * @delta: noncentrality parameter.
2764 * @x: reference value.
2765 *
2766 * Calculates the value at @x of the CDF of the noncentral Student t
2767 * distribution with @df dof and noncentrality parameter equal to
2768 * @delta. The algorithm is by Benson-Krishnamoorthy (2003) CSDA 43,
2769 * with minimal changes.
2770 *
2771 * Returns: the calculated probability, or #NADBL on failure.
2772 */
2773
nc_student_cdf(double df,double delta,double x)2774 double nc_student_cdf (double df, double delta, double x)
2775 {
2776 double errtol = 1.0e-16;
2777 double maxit = 512;
2778 double ax, y, del, dels, k, a, b, c;
2779 double pkf, pkb, qkf, qkb, pbetaf, pbetab, qbetaf, qbetab;
2780 double pgamf, pgamb, qgamf, qgamb, tmp, ret;
2781 double rempois, sum, ptermf, qtermf, ptermb, qtermb;
2782 int i;
2783
2784 if (df <= 0.0) {
2785 return NADBL;
2786 }
2787
2788 if (fabs(delta) <= 1.0e-10) {
2789 /*
2790 When non-centrality parameter is (essentially) zero, use
2791 ordinary t distribution
2792 */
2793 return student_cdf(df, x);
2794 }
2795
2796 ax = fabs(x);
2797
2798 del = x > 0 ? delta : -delta;
2799 ret = normal_cdf(-del);
2800
2801 if (ax < 1.0e-12) {
2802 return 1.0 - ret;
2803 }
2804
2805 /* settings */
2806
2807 y = ax*ax / (df + ax*ax);
2808 dels = del * del / 2.0;
2809 k = (int) floor(dels);
2810 a = k + 0.5;
2811 c = k + 1;
2812 b = df * 0.5;
2813
2814 /*
2815 Initialization to compute the P_k and Q_k terms
2816 and the respective incomplete beta functions
2817 */
2818
2819 tmp = -dels + k * log(dels);
2820 pkf = pkb = exp(tmp - lngamma(k + 1));
2821 qkf = qkb = exp(tmp - lngamma(k + 1.5));
2822 pbetaf = pbetab = incbet(a, b, y);
2823 qbetaf = qbetab = incbet(c, b, y);
2824
2825 /*
2826 Initialization to compute the incomplete beta functions
2827 associated with the P_i and the Q_i recursively:
2828 */
2829
2830 tmp = b * log(1-y) - lngamma(b);
2831 pgamf = exp(lngamma(a+b-1) - lngamma(a) + (a-1) * log(y) + tmp);
2832 pgamb = pgamf * y * (a + b - 1)/a;
2833
2834 qgamf = exp(lngamma(c+b-1) - lngamma(c) + (c-1) * log(y) + tmp);
2835 qgamb = qgamf * y * (c + b - 1)/c;
2836
2837 /*
2838 Compute the remainder of the Poisson weights
2839 */
2840
2841 rempois = 1.0 - pkf;
2842 sum = pkf * pbetaf + del * qkf * qbetaf * ISQRT_2;
2843
2844 for (i = 1; i<=k && rempois>errtol; i++) {
2845 /* first block --- backwards */
2846 pgamb *= (a-i+1)/(y * (a+b-i));
2847 pbetab += pgamb;
2848 pkb *= (k-i+1)/dels;
2849 ptermb = pkb * pbetab;
2850
2851 /* second block --- backwards */
2852 qgamb *= (c-i+1)/(y * (c+b-i));
2853 qbetab += qgamb;
2854 qkb *= (k-i+1.5)/dels;
2855 qtermb = qkb * qbetab;
2856
2857 /* accumulate */
2858 sum += ptermb + del * qtermb * ISQRT_2;
2859 rempois -= pkb;
2860 }
2861
2862 for (i = 1; i<maxit && rempois > errtol; i++) {
2863 /* first block --- forwards */
2864 pgamf *= y * (a+b-2+i)/(a+i-1);
2865 pbetaf -= pgamf;
2866 pkf *= dels/(k+i);
2867 ptermf = pkf * pbetaf;
2868
2869 /* second block --- forwards */
2870 qgamf *= y * (c+b-2+i)/(c+i-1);
2871 qbetaf -= qgamf;
2872 qkf *= dels/(k+i+0.5);
2873 qtermf = qkf * qbetaf;
2874
2875 /* accumulate */
2876 sum += ptermf + del * qtermf * ISQRT_2;
2877 rempois -= pkf;
2878 }
2879
2880 ret += sum/2.0;
2881
2882 return x < 0 ? (1.0 - ret) : ret;
2883 }
2884
2885 /**
2886 * nc_student_pdf:
2887 * @df: degrees of freedom.
2888 * @delta: noncentrality parameter.
2889 * @x: reference value.
2890 *
2891 * Calculates the value at @x of the PDF of the noncentral Student t
2892 * distribution with @df dof and noncentrality parameter equal to
2893 * @delta. The algorithm is from Wikipedia, apparently used in R too.
2894 *
2895 * Returns: the calculated density, or #NADBL on failure.
2896 */
2897
nc_student_pdf(double df,double delta,double x)2898 double nc_student_pdf (double df, double delta, double x)
2899 {
2900 double ret, tmp;
2901
2902 if (df <= 0.0) {
2903 return NADBL;
2904 }
2905
2906 if (fabs(delta) <= 1.0e-10) {
2907 /*
2908 When non-centrality parameter is (essentially) zero, use
2909 ordinary t distribution
2910 */
2911 return student_pdf(df, x);
2912 }
2913
2914 if (fabs(x) < 1.0e-12) {
2915 tmp = lngamma((df+1)/2) - lngamma(df/2);
2916 ret = exp(tmp - 0.5 * delta*delta) / (sqrt(M_PI * df));
2917 } else {
2918 tmp = nc_student_cdf(df+2, delta, x * sqrt(1 + 2.0/df)) -
2919 nc_student_cdf(df, delta, x);
2920 ret = tmp * (df / x);
2921 }
2922
2923 return ret;
2924 }
2925
2926
nct_pdf_array(double df,double delta,double * x,int n)2927 static int nct_pdf_array (double df, double delta, double *x, int n)
2928 {
2929 int i, err = 0;
2930
2931 if (df > 0) {
2932 for (i=0; i<n; i++) {
2933 if (!na(x[i])) {
2934 x[i] = nc_student_pdf(df, delta, x[i]);
2935 } else {
2936 x[i] = NADBL;
2937 }
2938 }
2939 } else {
2940 err = E_DATA;
2941 }
2942
2943 if (err) {
2944 for (i=0; i<n; i++) {
2945 x[i] = NADBL;
2946 }
2947 }
2948
2949 return err;
2950 }
2951
2952 /**
2953 * nct_cdf_inverse:
2954 * @p: degrees of freedom.
2955 * @c: noncentrality parameter.
2956 * @q: probability.
2957 *
2958 * Calculates the @q-th quantile of the noncentral Student t
2959 * distribution with @c dof and noncentrality parameter equal to @c
2960 * via a rough and not particularly clever root-finding
2961 * algorithm. Maybe this can be more efficient by using logs. Some
2962 * experimentation needed.
2963 *
2964 * Returns: the calculated quantile, or #NADBL on failure.
2965 */
2966
nct_cdf_inverse(double p,double c,double q)2967 static double nct_cdf_inverse (double p, double c, double q)
2968 {
2969 double x, d0, d1;
2970 int iter, subiter;
2971 double F, f, dir;
2972
2973 if (p < 1 || c < 0 || q <= 0 || q >= 1) {
2974 return NADBL;
2975 }
2976
2977 if (fabs(c) < 1.0e-10) {
2978 /* don't bother for infinitesimal c */
2979 return student_cdf_inverse(p, q);
2980 }
2981
2982 x = c + student_cdf_inverse(p, q) / sqrt(p - 0.5);
2983 d0 = 1.0e7;
2984 iter = 0;
2985
2986 while (fabs(d0) > 1.0e-10 && iter < 1000) {
2987 F = nc_student_cdf(p, c, x);
2988 f = nc_student_pdf(p, c, x);
2989 d0 = F - q;
2990 dir = d0/f;
2991 d1 = 1.0e7;
2992 subiter = 0;
2993
2994 while (fabs(d1) > fabs(d0) && subiter < 100) {
2995 d1 = F - nc_student_cdf(p, c, x - dir);
2996 dir /= 2.0;
2997 subiter++;
2998 }
2999
3000 if (subiter >= 100) {
3001 x = NADBL;
3002 break;
3003 } else {
3004 x -= dir*2;
3005 d0 = d1;
3006 iter++;
3007 }
3008 }
3009
3010 if (iter >= 1000) {
3011 x = NADBL;
3012 }
3013
3014 return x;
3015 }
3016
3017 struct distmap {
3018 int code;
3019 char *s;
3020 };
3021
dist_code_from_string(const char * s)3022 int dist_code_from_string (const char *s)
3023 {
3024 struct distmap dmap[] = {
3025 { D_UNIFORM, "u" },
3026 { D_UDISCRT, "i" },
3027 { D_NORMAL, "z" },
3028 { D_STUDENT, "t" },
3029 { D_CHISQ, "x" },
3030 { D_SNEDECOR, "f" },
3031 { D_BINOMIAL, "b" },
3032 { D_POISSON, "p" },
3033 { D_EXPON, "exp" },
3034 { D_WEIBULL, "w" },
3035 { D_GAMMA, "g" },
3036 { D_GED, "e" },
3037 { D_LAPLACE, "l" },
3038 { D_BETA, "beta" },
3039 { D_DW, "d" },
3040 { D_BINORM, "D" },
3041 { D_JOHANSEN, "J" },
3042 { D_BETABIN, "bb" },
3043 { D_NC_CHISQ, "ncx" },
3044 { D_NC_F, "ncf" },
3045 { D_NC_T, "nct" },
3046 { D_LOGISTIC, "s" },
3047 { D_NONE, NULL }
3048 };
3049 char test[8];
3050 int i;
3051
3052 if (!strcmp(s, "D")) {
3053 /* special: case counts for bivariate normal */
3054 return D_BINORM;
3055 }
3056
3057 /* otherwise we'll ignore case */
3058 for (i=0; i<8 && s[i]; i++) {
3059 test[i] = tolower(s[i]);
3060 }
3061 test[i] = '\0';
3062
3063 for (i=0; dmap[i].code; i++) {
3064 if (!strcmp(test, dmap[i].s)) {
3065 return dmap[i].code;
3066 }
3067 }
3068
3069 /* backward compatibility */
3070 if (!strcmp(test, "n")) {
3071 return D_NORMAL;
3072 } else if (!strcmp(test, "c")) {
3073 return D_CHISQ;
3074 } else if (!strcmp(test, "lgt")) {
3075 return D_LOGISTIC;
3076 }
3077
3078 return D_NONE;
3079 }
3080
3081 /**
3082 * print_critval:
3083 * @dist: distribution code.
3084 * @parm: array holding 0 to 2 parameter values.
3085 * @a: alpha.
3086 * @c: the critical value.
3087 * @prn: gretl printer.
3088 *
3089 * Prints the critical value information in a consistent manner.
3090 */
3091
print_critval(int dist,const double * parm,double a,double c,PRN * prn)3092 void print_critval (int dist, const double *parm, double a, double c, PRN *prn)
3093 {
3094 switch (dist) {
3095 case D_NORMAL:
3096 pprintf(prn, "%s", _("Standard normal distribution"));
3097 break;
3098 case D_STUDENT:
3099 pprintf(prn, "t(%g)", parm[0]);
3100 break;
3101 case D_CHISQ:
3102 pprintf(prn, "%s(%g)", _("Chi-square"), parm[0]);
3103 break;
3104 case D_SNEDECOR:
3105 pprintf(prn, "F(%g, %g)", parm[0], parm[1]);
3106 break;
3107 case D_BINOMIAL:
3108 pprintf(prn, "Binomial (P = %g, %g trials)", parm[0], parm[1]);
3109 break;
3110 case D_POISSON:
3111 pprintf(prn, "Poisson (mean = %g)", parm[0]);
3112 break;
3113 case D_EXPON:
3114 pprintf(prn, "Exponential (scale = %g)", parm[0]);
3115 break;
3116 case D_WEIBULL:
3117 pprintf(prn, "Weibull (shape = %g, scale = %g)", parm[0], parm[1]);
3118 break;
3119 }
3120
3121 pputs(prn, "\n ");
3122 pprintf(prn, _("right-tail probability = %g"), a);
3123 pputs(prn, "\n ");
3124 pprintf(prn, _("complementary probability = %g"), 1.0 - a);
3125 if (a < 0.5 && (dist == D_NORMAL || dist == D_STUDENT)) {
3126 pputs(prn, "\n ");
3127 pprintf(prn, _("two-tailed probability = %g"), 2.0 * a);
3128 }
3129 pputs(prn, "\n\n ");
3130 pprintf(prn, _("Critical value = %g"), c);
3131 pputc(prn, '\n');
3132 }
3133
3134 /* This apparatus is for use with the "batch p-value"
3135 routine: it remembers the parameters (p) and
3136 argument (x) from the last internal p-value
3137 assessment.
3138 */
3139
3140 static double pvargs[3];
3141
remember_pvalue_args(const double * p,double x)3142 static void remember_pvalue_args (const double *p, double x)
3143 {
3144 pvargs[0] = p[0];
3145 pvargs[1] = p[1];
3146 pvargs[2] = x;
3147 }
3148
3149 /* end remember parameters */
3150
pdist_check_input(int dist,const double * parm,double x)3151 static int pdist_check_input (int dist, const double *parm,
3152 double x)
3153 {
3154 int i, np = 1; /* default */
3155
3156 if (na(x)) {
3157 return E_MISSDATA;
3158 }
3159
3160 if (dist == D_NORMAL) {
3161 np = 0;
3162 } else if (dist == D_SNEDECOR || dist == D_GAMMA ||
3163 dist == D_BINOMIAL || dist == D_WEIBULL ||
3164 dist == D_NC_CHISQ || dist == D_NC_T ||
3165 dist == D_LAPLACE || dist == D_BETA) {
3166 np = 2;
3167 } else if (dist == D_JOHANSEN || dist == D_BETABIN ||
3168 dist == D_NC_F) {
3169 np = 3;
3170 }
3171
3172 for (i=0; i<np; i++) {
3173 if (na(parm[i])) {
3174 return E_MISSDATA;
3175 }
3176 }
3177
3178 return 0;
3179 }
3180
3181
3182 /**
3183 * gretl_get_cdf_inverse:
3184 * @dist: distribution code.
3185 * @parm: array holding from zero to two parameter values,
3186 * depending on the distribution.
3187 * @a: probability value.
3188 *
3189 * Returns: the argument, y, for which the area under the PDF
3190 * specified by @dist and @parm, integrated from its minimum to y,
3191 * is equal to @a, or #NADBL on failure.
3192 */
3193
gretl_get_cdf_inverse(int dist,const double * parm,double a)3194 double gretl_get_cdf_inverse (int dist, const double *parm,
3195 double a)
3196 {
3197 double y = NADBL;
3198
3199 if (pdist_check_input(dist, parm, a) == E_MISSDATA) {
3200 return y;
3201 }
3202
3203 if (dist == D_NORMAL) {
3204 y = normal_cdf_inverse(a);
3205 } else if (dist == D_STUDENT) {
3206 y = student_cdf_inverse(parm[0], a);
3207 } else if (dist == D_CHISQ) {
3208 y = chisq_cdf_inverse((int) parm[0], a);
3209 } else if (dist == D_GAMMA) {
3210 y = gamma_cdf_inverse(parm[0], parm[1], a);
3211 } else if (dist == D_SNEDECOR) {
3212 y = snedecor_cdf_inverse(parm[0], parm[1], a);
3213 } else if (dist == D_BINOMIAL) {
3214 y = binomial_cdf_inverse(parm[0], (int) parm[1], a);
3215 } else if (dist == D_POISSON) {
3216 y = poisson_cdf_inverse(parm[0], a);
3217 } else if (dist == D_GED) {
3218 y = GED_cdf_inverse(parm[0], a);
3219 } else if (dist == D_LAPLACE) {
3220 y = laplace_cdf_inverse(parm[0], parm[1], a);
3221 } else if (dist == D_NC_F) {
3222 y = ncf_cdf_inverse(parm[0], parm[1], parm[2], a);
3223 } else if (dist == D_NC_CHISQ) {
3224 y = nc_chisq_cdf_inverse(parm[0], parm[1], a);
3225 } else if (dist == D_NC_T) {
3226 y = nct_cdf_inverse(parm[0], parm[1], a);
3227 }
3228
3229 return y;
3230 }
3231
3232
3233 /**
3234 * gretl_get_critval:
3235 * @dist: distribution code.
3236 * @parm: array holding from zero to two parameter values,
3237 * depending on the distribution.
3238 * @a: right-tail probability.
3239 *
3240 * Returns: the abcsissa value x for the distribution specified
3241 * by @dist and @parm, such that P(X >= x) = @a, or #NADBL on
3242 * failure.
3243 */
3244
gretl_get_critval(int dist,const double * parm,double a)3245 double gretl_get_critval (int dist, const double *parm, double a)
3246 {
3247 double x = NADBL;
3248
3249 if (pdist_check_input(dist, parm, a) == E_MISSDATA) {
3250 return x;
3251 }
3252
3253 if (dist == D_NORMAL) {
3254 x = normal_critval(a);
3255 } else if (dist == D_STUDENT) {
3256 x = student_critval(parm[0], a);
3257 } else if (dist == D_CHISQ) {
3258 x = chisq_critval((int) parm[0], a);
3259 } else if (dist == D_SNEDECOR) {
3260 x = snedecor_critval((int) parm[0], (int) parm[1], a);
3261 } else if (dist == D_BINOMIAL) {
3262 x = binomial_critval(parm[0], (int) parm[1], a);
3263 } else if (dist == D_POISSON) {
3264 x = poisson_critval(parm[0], a);
3265 } else if (dist == D_WEIBULL) {
3266 x = weibull_critval(parm[0], parm[1], a);
3267 } else if (dist == D_EXPON) {
3268 /* special case of Weibull */
3269 x = weibull_critval(1.0, parm[0], a);
3270 } else if (dist == D_GED) {
3271 x = GED_cdf_inverse(parm[0], 1-a);
3272 } else if (dist == D_LAPLACE) {
3273 x = laplace_cdf_inverse(parm[0], parm[1], 1-a);
3274 }
3275
3276 return x;
3277 }
3278
3279 /**
3280 * gretl_get_cdf:
3281 * @dist: distribution code.
3282 * @parm: array holding from zero to two parameter values,
3283 * depending on the distribution.
3284 * @x: abscissa value.
3285 *
3286 * Evaluates the CDF for the distribution specified by
3287 * @dist and @parm applicable at @x.
3288 *
3289 * Returns: the CDF value, or #NADBL on error.
3290 */
3291
gretl_get_cdf(int dist,const double * parm,double x)3292 double gretl_get_cdf (int dist, const double *parm, double x)
3293 {
3294 double y = NADBL;
3295
3296 if (pdist_check_input(dist, parm, x) == E_MISSDATA) {
3297 return y;
3298 }
3299
3300 if (dist == D_NORMAL) {
3301 y = normal_cdf(x);
3302 } else if (dist == D_STUDENT) {
3303 y = student_cdf(parm[0], x);
3304 } else if (dist == D_CHISQ) {
3305 y = chisq_cdf((int) parm[0], x);
3306 } else if (dist == D_SNEDECOR) {
3307 y = snedecor_cdf((int) parm[0], (int) parm[1], x);
3308 } else if (dist == D_GAMMA) {
3309 y = gamma_cdf(parm[0], parm[1], x, 1);
3310 } else if (dist == D_BINOMIAL) {
3311 y = binomial_cdf(parm[0], (int) parm[1], (int) x);
3312 } else if (dist == D_POISSON) {
3313 y = poisson_cdf(parm[0], (int) x);
3314 } else if (dist == D_EXPON) {
3315 y = exponential_cdf(parm[0], x);
3316 } else if (dist == D_WEIBULL) {
3317 y = weibull_cdf(parm[0], parm[1], x);
3318 } else if (dist == D_GED) {
3319 y = GED_cdf(parm[0], x);
3320 } else if (dist == D_LAPLACE) {
3321 y = laplace_cdf(parm[0], parm[1], x);
3322 } else if (dist == D_NC_CHISQ) {
3323 y = nc_chisq_cdf(parm[0], parm[1], x);
3324 } else if (dist == D_NC_F) {
3325 y = nc_snedecor_cdf(parm[0], parm[1], parm[2], x);
3326 } else if (dist == D_NC_T) {
3327 y = nc_student_cdf(parm[0], parm[1], x);
3328 } else if (dist == D_LOGISTIC) {
3329 y = logistic_cdf(x);
3330 } else if (dist == D_BETA) {
3331 y = beta_cdf(parm[0], parm[1], x);
3332 }
3333
3334 return y;
3335 }
3336
3337 /**
3338 * gretl_get_pdf:
3339 * @dist: distribution code.
3340 * @parm: array holding from zero to two parameter values,
3341 * depending on the distribution.
3342 * @x: abscissa value.
3343 *
3344 * Evaluates the PDF for the distribution specified by
3345 * @dist and @parm at @x.
3346 *
3347 * Returns: the PDF value, or #NADBL on error.
3348 */
3349
gretl_get_pdf(int dist,const double * parm,double x)3350 double gretl_get_pdf (int dist, const double *parm, double x)
3351 {
3352 double y = NADBL;
3353
3354 if (pdist_check_input(dist, parm, x) == E_MISSDATA) {
3355 return y;
3356 }
3357
3358 if (dist == D_NORMAL) {
3359 y = normal_pdf(x);
3360 } else if (dist == D_STUDENT) {
3361 y = student_pdf(parm[0], x);
3362 } else if (dist == D_CHISQ) {
3363 y = chisq_pdf((int) parm[0], x);
3364 } else if (dist == D_SNEDECOR) {
3365 y = snedecor_pdf((int) parm[0], (int) parm[1], x);
3366 } else if (dist == D_GAMMA) {
3367 y = gamma_pdf(parm[0], parm[1], x);
3368 } else if (dist == D_BINOMIAL) {
3369 y = binomial_pmf(parm[0], parm[1], x);
3370 } else if (dist == D_POISSON) {
3371 y = poisson_pmf(parm[0], x);
3372 } else if (dist == D_EXPON) {
3373 y = exponential_pdf(parm[0], x);
3374 } else if (dist == D_WEIBULL) {
3375 y = weibull_pdf(parm[0], parm[1], x);
3376 } else if (dist == D_GED) {
3377 y = GED_pdf(parm[0], x);
3378 } else if (dist == D_LAPLACE) {
3379 y = laplace_pdf(parm[0], parm[1], x);
3380 } else if (dist == D_NC_F) {
3381 y = ncf_pdf(parm[0], parm[1], parm[2], x);
3382 } else if (dist == D_NC_T) {
3383 y = nc_student_pdf(parm[0], parm[1], x);
3384 } else if (dist == D_NC_CHISQ) {
3385 y = nc_chisq_pdf(parm[0], parm[1], x);
3386 } else if (dist == D_BETA) {
3387 y = beta_pdf(parm[0], parm[1], x);
3388 }
3389
3390 return y;
3391 }
3392
3393 /**
3394 * gretl_fill_pdf_array:
3395 * @dist: distribution code.
3396 * @parm: array holding from zero to two parameter values,
3397 * depending on the distribution.
3398 * @x: see below.
3399 * @n: number of elements in @x.
3400 *
3401 * On input, @x contains an array of abscissae at which the
3402 * PDF specified by @dist and @parm should be evaluated. On
3403 * output it contains the corresponding PDF values.
3404 *
3405 * Returns: 0 on success, non-zero on error.
3406 */
3407
gretl_fill_pdf_array(int dist,const double * parm,double * x,int n)3408 int gretl_fill_pdf_array (int dist, const double *parm,
3409 double *x, int n)
3410 {
3411 int err = E_DATA;
3412
3413 if (pdist_check_input(dist, parm, 0) == E_MISSDATA) {
3414 return E_MISSDATA;
3415 }
3416
3417 if (dist == D_NORMAL) {
3418 err = normal_pdf_array(x, n);
3419 } else if (dist == D_STUDENT) {
3420 err = student_pdf_array(parm[0], x, n);
3421 } else if (dist == D_CHISQ) {
3422 err = chisq_pdf_array(parm[0], x, n);
3423 } else if (dist == D_SNEDECOR) {
3424 err = snedecor_pdf_array((int) parm[0], (int) parm[1], x, n);
3425 } else if (dist == D_GAMMA) {
3426 err = gamma_pdf_array(parm[0], parm[1], x, n);
3427 } else if (dist == D_BINOMIAL) {
3428 err = binomial_pmf_array(parm[0], parm[1], x, n);
3429 } else if (dist == D_POISSON) {
3430 err = poisson_pmf_array(parm[0], x, n);
3431 } else if (dist == D_EXPON) {
3432 err = exponential_pdf_array(parm[0], x, n);
3433 } else if (dist == D_WEIBULL) {
3434 err = weibull_pdf_array(parm[0], parm[1], x, n);
3435 } else if (dist == D_GED) {
3436 err = GED_pdf_array(parm[0], x, n);
3437 } else if (dist == D_LAPLACE) {
3438 err = laplace_pdf_array(parm[0], parm[1], x, n);
3439 } else if (dist == D_NC_F) {
3440 err = ncf_pdf_array(parm[0], parm[1], parm[2], x, n);
3441 } else if (dist == D_NC_T) {
3442 err = nct_pdf_array(parm[0], parm[1], x, n);
3443 } else if (dist == D_NC_CHISQ) {
3444 err = nc_chisq_pdf_array(parm[0], parm[1], x, n);
3445 } else if (dist == D_BETA) {
3446 err = beta_pdf_array(parm[0], parm[1], x, n);
3447 }
3448
3449 return err;
3450 }
3451
3452 /**
3453 * gretl_get_pvalue:
3454 * @dist: distribution code.
3455 * @parm: array holding from zero to two parameter values,
3456 * depending on the distribution.
3457 * @x: abscissa value.
3458 *
3459 * Returns: the integral of the PDF specified by @dist and
3460 * @parm from @x to infinity, or #NADBL on error.
3461 */
3462
gretl_get_pvalue(int dist,const double * parm,double x)3463 double gretl_get_pvalue (int dist, const double *parm, double x)
3464 {
3465 double y = NADBL;
3466
3467 if (pdist_check_input(dist, parm, x) == E_MISSDATA) {
3468 return y;
3469 }
3470
3471 if (dist == D_NORMAL) {
3472 y = normal_cdf_comp(x);
3473 } else if (dist == D_STUDENT) {
3474 y = student_cdf_comp(parm[0], x);
3475 } else if (dist == D_CHISQ) {
3476 y = chisq_cdf_comp((int) parm[0], x);
3477 } else if (dist == D_SNEDECOR) {
3478 y = snedecor_cdf_comp(parm[0], parm[1], x);
3479 } else if (dist == D_GAMMA) {
3480 y = gamma_cdf_comp(parm[0], parm[1], x, 1);
3481 } else if (dist == D_BINOMIAL) {
3482 y = binomial_cdf_comp(parm[0], (int) parm[1], x);
3483 } else if (dist == D_POISSON) {
3484 y = poisson_cdf_comp(parm[0], x);
3485 } else if (dist == D_EXPON) {
3486 y = weibull_cdf_comp(1.0, parm[0], x);
3487 } else if (dist == D_WEIBULL) {
3488 y = weibull_cdf_comp(parm[0], parm[1], x);
3489 } else if (dist == D_GED) {
3490 y = GED_cdf_comp(parm[0], x);
3491 } else if (dist == D_LAPLACE) {
3492 y = laplace_cdf_comp(parm[0], parm[1], x);
3493 } else if (dist == D_JOHANSEN) {
3494 y = johansen_trace_pval((int) parm[0], (int) parm[1],
3495 (int) parm[2], x);
3496 }
3497
3498 remember_pvalue_args(parm, x);
3499
3500 return y;
3501 }
3502
gretl_fill_random_array(double * x,int t1,int t2,int dist,const double * parm,const double * vecp1,const double * vecp2)3503 static int gretl_fill_random_array (double *x, int t1, int t2,
3504 int dist, const double *parm,
3505 const double *vecp1,
3506 const double *vecp2)
3507 {
3508 int t, err = 0;
3509
3510 if (dist == D_UNIFORM) {
3511 /* uniform, continuous */
3512 double min = parm[0], max = parm[1];
3513
3514 if (vecp1 != NULL || vecp2 != NULL) {
3515 for (t=t1; t<=t2 && !err; t++) {
3516 if (vecp1 != NULL) min = vecp1[t];
3517 if (vecp2 != NULL) max = vecp2[t];
3518 err = gretl_rand_uniform_minmax(x, t, t, min, max);
3519 }
3520 } else {
3521 err = gretl_rand_uniform_minmax(x, t1, t2, min, max);
3522 }
3523 } else if (dist == D_UDISCRT) {
3524 /* uniform, discrete */
3525 int min = (int) parm[0], max = (int) parm[1];
3526
3527 if (vecp1 != NULL || vecp2 != NULL) {
3528 for (t=t1; t<=t2 && !err; t++) {
3529 if (vecp1 != NULL) min = (int) vecp1[t];
3530 if (vecp2 != NULL) max = (int) vecp2[t];
3531 err = gretl_rand_uniform_int_minmax(x, t, t, min, max,
3532 OPT_NONE);
3533 }
3534 } else {
3535 err = gretl_rand_uniform_int_minmax(x, t1, t2, min, max,
3536 OPT_NONE);
3537 }
3538 } else if (dist == D_NORMAL) {
3539 double mu = parm[0], sd = parm[1];
3540
3541 if (vecp1 != NULL || vecp2 != NULL) {
3542 for (t=t1; t<=t2 && !err; t++) {
3543 if (vecp1 != NULL) mu = vecp1[t];
3544 if (vecp2 != NULL) sd = vecp2[t];
3545 err = gretl_rand_normal_full(x, t, t, mu, sd);
3546 }
3547 } else {
3548 err = gretl_rand_normal_full(x, t1, t2, mu, sd);
3549 }
3550 } else if (dist == D_STUDENT) {
3551 /* Student's t */
3552 double v = parm[0];
3553
3554 if (vecp1 != NULL) {
3555 for (t=t1; t<=t2 && !err; t++) {
3556 v = vecp1[t];
3557 err = gretl_rand_student(x, t, t, v);
3558 }
3559 } else {
3560 err = gretl_rand_student(x, t1, t2, v);
3561 }
3562 } else if (dist == D_CHISQ) {
3563 /* chi-square */
3564 int v = parm[0];
3565
3566 if (vecp1 != NULL) {
3567 for (t=t1; t<=t2 && !err; t++) {
3568 v = vecp1[t];
3569 err = gretl_rand_chisq(x, t, t, v);
3570 }
3571 } else {
3572 err = gretl_rand_chisq(x, t1, t2, v);
3573 }
3574 } else if (dist == D_SNEDECOR) {
3575 int v1 = parm[0], v2 = parm[1];
3576
3577 if (vecp1 != NULL || vecp2 != NULL) {
3578 for (t=t1; t<=t2 && !err; t++) {
3579 if (vecp1 != NULL) v1 = vecp1[t];
3580 if (vecp2 != NULL) v2 = vecp2[t];
3581 err = gretl_rand_F(x, t, t, v1, v2);
3582 }
3583 } else {
3584 err = gretl_rand_F(x, t1, t2, v1, v2);
3585 }
3586 } else if (dist == D_GAMMA) {
3587 double shape = parm[0], scale = parm[1];
3588
3589 if (vecp1 != NULL || vecp2 != NULL) {
3590 for (t=t1; t<=t2 && !err; t++) {
3591 if (vecp1 != NULL) shape = vecp1[t];
3592 if (vecp2 != NULL) scale = vecp2[t];
3593 err = gretl_rand_gamma(x, t, t, shape, scale);
3594 }
3595 } else {
3596 err = gretl_rand_gamma(x, t1, t2, shape, scale);
3597 }
3598 } else if (dist == D_BINOMIAL) {
3599 double pr = parm[0];
3600 int n = parm[1];
3601
3602 if (vecp1 != NULL || vecp2 != NULL) {
3603 for (t=t1; t<=t2 && !err; t++) {
3604 if (vecp1 != NULL) pr = vecp1[t];
3605 if (vecp2 != NULL) n = vecp2[1];
3606 err = gretl_rand_binomial(x, t, t, n, pr);
3607 }
3608 } else {
3609 err = gretl_rand_binomial(x, t1, t2, n, pr);
3610 }
3611 } else if (dist == D_POISSON) {
3612 double m = parm[0];
3613
3614 if (vecp1 != NULL) {
3615 for (t=t1; t<=t2 && !err; t++) {
3616 m = vecp1[t];
3617 err = gretl_rand_poisson(x, t, t, &m, 0);
3618 }
3619 } else {
3620 err = gretl_rand_poisson(x, t1, t2, &m, 0);
3621 }
3622 } else if (dist == D_EXPON) {
3623 double scale = parm[0];
3624
3625 if (vecp1 != NULL) {
3626 for (t=t1; t<=t2 && !err; t++) {
3627 scale = vecp1[t];
3628 err = gretl_rand_exponential(x, t, t, scale);
3629 }
3630 } else {
3631 err = gretl_rand_exponential(x, t1, t2, scale);
3632 }
3633 } else if (dist == D_WEIBULL) {
3634 double shape = parm[0], scale = parm[1];
3635
3636 if (vecp1 != NULL || vecp2 != NULL) {
3637 for (t=t1; t<=t2 && !err; t++) {
3638 if (vecp1 != NULL) shape = vecp1[t];
3639 if (vecp2 != NULL) scale = vecp2[t];
3640 err = gretl_rand_weibull(x, t, t, shape, scale);
3641 }
3642 } else {
3643 err = gretl_rand_weibull(x, t1, t2, shape, scale);
3644 }
3645 } else if (dist == D_GED) {
3646 double nu = parm[0];
3647
3648 if (vecp1 != NULL) {
3649 for (t=t1; t<=t2 && !err; t++) {
3650 if (vecp1 != NULL) nu = vecp1[t];
3651 err = gretl_rand_GED(x, t, t, nu);
3652 }
3653 } else {
3654 err = gretl_rand_GED(x, t1, t2, nu);
3655 }
3656 } else if (dist == D_LAPLACE) {
3657 double mu = parm[0], b = parm[1];
3658
3659 if (vecp1 != NULL || vecp2 != NULL) {
3660 for (t=t1; t<=t2 && !err; t++) {
3661 if (vecp1 != NULL) mu = vecp1[t];
3662 if (vecp2 != NULL) b = vecp1[t];
3663 err = gretl_rand_laplace(x, t, t, mu, b);
3664 }
3665 } else {
3666 err = gretl_rand_laplace(x, t1, t2, mu, b);
3667 }
3668 } else if (dist == D_BETA) {
3669 double shape1 = parm[0], shape2 = parm[1];
3670
3671 if (vecp1 != NULL || vecp2 != NULL) {
3672 for (t=t1; t<=t2 && !err; t++) {
3673 if (vecp1 != NULL) shape1 = vecp1[t];
3674 if (vecp2 != NULL) shape2 = vecp2[t];
3675 err = gretl_rand_beta(x, t, t, shape1, shape2);
3676 }
3677 } else {
3678 err = gretl_rand_beta(x, t1, t2, shape1, shape2);
3679 }
3680 } else if (dist == D_BETABIN) {
3681 double shape1 = parm[1], shape2 = parm[2];
3682 int n = parm[0];
3683
3684 err = gretl_rand_beta_binomial(x, t1, t2, n, shape1, shape2);
3685 } else if (dist == D_LOGISTIC) {
3686 double loc = parm[0], shape = parm[1];
3687
3688 err = gretl_rand_logistic(x, t1, t2, loc, shape);
3689 }
3690
3691 return err;
3692 }
3693
3694 /**
3695 * gretl_fill_random_series:
3696 * @x: series to fill (must be of length dset->n).
3697 * @dist: distribution code.
3698 * @parm: array holding either one or two scalar
3699 * parameter values, depending on the distribution.
3700 * @vecp1: series containing values for first param,
3701 * or %NULL.
3702 * @vecp2: series containing values for second param,
3703 * or %NULL.
3704 * @dset: dataset information.
3705 *
3706 * Fills @x with random values conforming to the distribution
3707 * given by @dist, which may require specification of either
3708 * one or two parameters. These parameters are either
3709 * given as (scalar) elements of @parm, or they may vary
3710 * by observation, in which case they are given as the
3711 * elements of @vecp1 or @vecp2.
3712 *
3713 * Returns: 0 on success, non-zero code on error.
3714 */
3715
gretl_fill_random_series(double * x,int dist,const double * parm,const double * vecp1,const double * vecp2,const DATASET * dset)3716 int gretl_fill_random_series (double *x, int dist,
3717 const double *parm,
3718 const double *vecp1,
3719 const double *vecp2,
3720 const DATASET *dset)
3721 {
3722 return gretl_fill_random_array(x, dset->t1, dset->t2,
3723 dist, parm, vecp1, vecp2);
3724 }
3725
gretl_get_random_matrix(int dist,const double * parm,const double * vecp1,const double * vecp2,int rows,int cols,int * err)3726 gretl_matrix *gretl_get_random_matrix (int dist,
3727 const double *parm,
3728 const double *vecp1,
3729 const double *vecp2,
3730 int rows, int cols,
3731 int *err)
3732 {
3733 gretl_matrix *m = NULL;
3734 int n = rows * cols;
3735
3736 if (n <= 0) {
3737 *err = E_INVARG;
3738 } else {
3739 m = gretl_matrix_alloc(rows, cols);
3740 if (m == NULL) {
3741 *err = E_ALLOC;
3742 return NULL;
3743 } else {
3744 *err = gretl_fill_random_array(m->val, 0, n-1, dist,
3745 parm, vecp1, vecp2);
3746 }
3747 }
3748
3749 return m;
3750 }
3751
gretl_get_random_scalar(int dist,const double * parm,int * err)3752 double gretl_get_random_scalar (int dist, const double *parm,
3753 int *err)
3754 {
3755 double x;
3756
3757 *err = gretl_fill_random_array(&x, 0, 0, dist,
3758 parm, NULL, NULL);
3759
3760 return x;
3761 }
3762
3763 static int
print_pv_string(double x,double p,PRN * prn)3764 print_pv_string (double x, double p, PRN *prn)
3765 {
3766 char numstr[32];
3767
3768 if (na(p)) {
3769 pprintf(prn, _("area to the right of %g: NA\n"), x);
3770 return 1;
3771 }
3772
3773 sprintf(numstr, "%g", p);
3774
3775 if (!strcmp(numstr, "1") || !strcmp(numstr, "0")) {
3776 pprintf(prn, _("area to the right of %g =~ %g\n"), x, p);
3777 } else {
3778 pprintf(prn, _("area to the right of %g = %g\n"), x, p);
3779 }
3780
3781 return 0;
3782 }
3783
3784 /**
3785 * print_pvalue:
3786 * @dist: distribution code.
3787 * @parm: array holding 1 or 2 parameter values.
3788 * @x: the value in the distribution.
3789 * @pv: the p-value.
3790 * @prn: gretl printer.
3791 *
3792 * Prints the p-value information in a consistent manner.
3793 */
3794
print_pvalue(int dist,const double * parm,double x,double pv,PRN * prn)3795 void print_pvalue (int dist, const double *parm, double x,
3796 double pv, PRN *prn)
3797 {
3798 double pc;
3799 int err;
3800
3801 switch (dist) {
3802
3803 case D_NORMAL:
3804 pprintf(prn, "%s: ", _("Standard normal"));
3805 err = print_pv_string(x, pv, prn);
3806 if (err) return;
3807 if (pv < 0.5) {
3808 pprintf(prn, _("(two-tailed value = %g; complement = %g)\n"),
3809 2 * pv, 1 - 2 * pv);
3810 } else {
3811 pc = normal_cdf(x);
3812 pprintf(prn, _("(to the left: %g)\n"), pc);
3813 pprintf(prn, _("(two-tailed value = %g; complement = %g)\n"),
3814 2 * pc, 1 - 2 * pc);
3815 }
3816 break;
3817
3818 case D_STUDENT:
3819 pprintf(prn, "t(%d): ", (int) parm[0]);
3820 err = print_pv_string(x, pv, prn);
3821 if (err) return;
3822 if (pv < 0.5) {
3823 pprintf(prn, _("(two-tailed value = %g; complement = %g)\n"),
3824 2 * pv, 1 - 2 * pv);
3825 } else {
3826 pc = student_cdf(parm[0], x);
3827 pprintf(prn, _("(to the left: %g)\n"), pc);
3828 pprintf(prn, _("(two-tailed value = %g; complement = %g)\n"),
3829 2 * pc, 1 - 2 * pc);
3830 }
3831 break;
3832
3833 case D_CHISQ:
3834 pprintf(prn, "%s(%d): ", _("Chi-square"), (int) parm[0]);
3835 err = print_pv_string(x, pv, prn);
3836 if (err) return;
3837 pc = chisq_cdf(parm[0], x);
3838 pprintf(prn, _("(to the left: %g)\n"), pc);
3839 break;
3840
3841 case D_SNEDECOR:
3842 pprintf(prn, "F(%d, %d): ", (int) parm[0], (int) parm[1]);
3843 err = print_pv_string(x, pv, prn);
3844 if (err) return;
3845 pc = snedecor_cdf((int) parm[0], (int) parm[1], x);
3846 pprintf(prn, _("(to the left: %g)\n"), pc);
3847 break;
3848
3849 case D_GAMMA:
3850 pprintf(prn, _("Gamma (shape %g, scale %g, mean %g, variance %g):"
3851 "\n area to the right of %g = %g\n"),
3852 parm[0], parm[1], parm[0] * parm[1],
3853 parm[0] * parm[1] * parm[1],
3854 x, pv);
3855 break;
3856
3857 case D_BINOMIAL:
3858 pprintf(prn, _("Binomial (p = %g, n = %d):"
3859 "\n Prob(x > %d) = %g\n"),
3860 parm[0], (int) parm[1], (int) x, pv);
3861 pc = binomial_cdf(parm[0], parm[1], x);
3862 if (x > 0) {
3863 pprintf(prn, _(" Prob(x <= %d) = %g\n"), (int) x, pc);
3864 pprintf(prn, _(" Prob(x = %d) = %g\n"), (int) x,
3865 pc - binomial_cdf(parm[0], parm[1], x - 1));
3866 } else {
3867 pprintf(prn, _(" Prob(x = %d) = %g\n"), (int) x, pc);
3868 }
3869 break;
3870
3871 case D_POISSON:
3872 pprintf(prn, _("Poisson (mean = %g): "), parm[0]);
3873 err = print_pv_string(x, pv, prn);
3874 if (err) return;
3875 pc = poisson_cdf(parm[0], (int) x);
3876 if (x > 0) {
3877 pprintf(prn, _(" Prob(x <= %d) = %g\n"), (int) x, pc);
3878 pprintf(prn, _(" Prob(x = %d) = %g\n"), (int) x,
3879 poisson_pmf(parm[0], x));
3880 } else {
3881 pprintf(prn, _(" Prob(x = %d) = %g\n"), (int) x, pc);
3882 }
3883 break;
3884
3885 case D_EXPON:
3886 pprintf(prn, _("Exponential (scale = %g): "), parm[0]);
3887 err = print_pv_string(x, pv, prn);
3888 if (err) return;
3889 pc = exponential_cdf(parm[0], x);
3890 pprintf(prn, _("(to the left: %g)\n"), pc);
3891 break;
3892
3893 case D_WEIBULL:
3894 pprintf(prn, _("Weibull (shape = %g, scale = %g): "),
3895 parm[0], parm[1]);
3896 err = print_pv_string(x, pv, prn);
3897 if (err) return;
3898 pc = weibull_cdf(parm[0], parm[1], x);
3899 pprintf(prn, _("(to the left: %g)\n"), pc);
3900 break;
3901
3902 case D_GED:
3903 pprintf(prn, _("GED (shape = %g): "), parm[0]);
3904 err = print_pv_string(x, pv, prn);
3905 if (err) return;
3906 pc = GED_cdf(parm[0], x);
3907 pprintf(prn, _("(to the left: %g)\n"), pc);
3908 break;
3909
3910 case D_LAPLACE:
3911 pprintf(prn, _("Laplace (mean = %g, scale = %g): "), parm[0], parm[1]);
3912 err = print_pv_string(x, pv, prn);
3913 if (err) return;
3914 pc = laplace_cdf(parm[0], parm[1], x);
3915 pprintf(prn, _("(to the left: %g)\n"), pc);
3916 break;
3917
3918 default:
3919 break;
3920 }
3921 }
3922
3923 /**
3924 * batch_pvalue:
3925 * @str: the command line, which should be of one of the following forms:
3926 * pvalue z x (Normal distribution);
3927 * pvalue t df x (t-distribution);
3928 * pvalue X df x (Chi-square);
3929 * pvalue F dfn dfd x (F-distribution); or
3930 * pvalue G mean variance x (Gamma distribution).
3931 * pvalue B prob n x (Binomial distribution).
3932 * pvalue P mean k (Poisson distribution).
3933 * pvalue W shape scale x (Weibull distribution).
3934 * @dset: dataset struct.
3935 * @prn: gretl printing struct.
3936 *
3937 * Calculates and prints the probability that a random variable
3938 * distributed as specified in the command line @str exceeds the
3939 * value indicated in @str.
3940 *
3941 * Returns: 0 on success, non-zero code on error.
3942 */
3943
batch_pvalue(const char * str,DATASET * dset,PRN * prn)3944 int batch_pvalue (const char *str, DATASET *dset, PRN *prn)
3945 {
3946 double pv = NADBL;
3947 char line[MAXLEN];
3948 char **S;
3949 int dist;
3950 int i, n, m;
3951 int err = 0;
3952
3953 if (str == NULL || *str == '\0') {
3954 return E_ARGS;
3955 }
3956
3957 if (!strncmp(str, "pvalue ", 7)) {
3958 str += 7;
3959 }
3960
3961 S = gretl_string_split(str, &n, NULL);
3962 if (S == NULL) {
3963 return E_ALLOC;
3964 }
3965
3966 dist = dist_code_from_string(S[0]);
3967
3968 if (dist == D_NONE) {
3969 err = E_INVARG;
3970 } else {
3971 strcpy(line, "pvalue(");
3972 m = 8;
3973 for (i=0; i<n && !err; i++) {
3974 m += strlen(S[i]) + 1;
3975 if (m > MAXLEN) {
3976 err = E_DATA;
3977 } else {
3978 strcat(line, S[i]);
3979 strcat(line, (i == n - 1)? ")" : ",");
3980 }
3981 }
3982 }
3983
3984 strings_array_free(S, n);
3985
3986 if (!err) {
3987 pv = generate_scalar(line, dset, &err);
3988 }
3989
3990 if (!err) {
3991 print_pvalue(dist, pvargs, pvargs[2], pv, prn);
3992 }
3993
3994 return err;
3995 }
3996
3997 /**
3998 * gretl_get_DW:
3999 * @n: number of observations
4000 * @k: number of regressors excluding the constant.
4001 * @err: location to receive error code.
4002 *
4003 * Consults a table of Durbin-Watson critical values and
4004 * returns the results in a gretl_matrix.
4005 *
4006 * Returns: on success, a 4-vector containing the lower
4007 * and upper Durbin-Watson values, dl and du, along with
4008 * the values actually used for @n and @k (which may differ
4009 * from those given on input if the exact values are not
4010 * found in the table and have to be approximated).
4011 * On error, returns %NULL.
4012 */
4013
gretl_get_DW(int n,int k,int * err)4014 gretl_matrix *gretl_get_DW (int n, int k, int *err)
4015 {
4016 int (*dw_lookup) (int, int, gretl_matrix **);
4017 gretl_matrix *m = NULL;
4018
4019 dw_lookup = get_plugin_function("dw_lookup");
4020
4021 if (dw_lookup == NULL) {
4022 *err = E_FOPEN;
4023 return NULL;
4024 }
4025
4026 *err = (*dw_lookup) (n, k, &m);
4027
4028 return m;
4029 }
4030