3 // PSPP - a program for statistical analysis.
4 // Copyright (C) 2005, 2006, 2009, 2010, 2011, 2012, 2015, 2016 Free Software Foundation, Inc.
6 // This program is free software: you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation, either version 3 of the License, or
9 // (at your option) any later version.
11 // This program is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
16 // You should have received a copy of the GNU General Public License
17 // along with this program. If not, see <http://www.gnu.org/licenses/>.
19 operator NEG (x) = -x;
21 operator ADD (a, b) = a + b;
22 operator SUB (a, b) = a - b;
24 absorb_miss operator MUL (a, b)
25 = (a == 0. || b == 0. ? 0.
26 : a == SYSMIS || b == SYSMIS ? SYSMIS
29 absorb_miss operator DIV (a, b)
31 : a == SYSMIS || b == SYSMIS ? SYSMIS
34 absorb_miss operator POW (a, b)
35 = (a == SYSMIS ? (b == 0. ? 1. : a)
36 : b == SYSMIS ? (a == 0. ? 0. : SYSMIS)
37 : a == 0. && b <= 0. ? SYSMIS
40 absorb_miss boolean operator AND (boolean a, boolean b)
43 : b == SYSMIS ? SYSMIS
46 absorb_miss boolean operator OR (boolean a, boolean b)
49 : b == SYSMIS ? SYSMIS
52 boolean operator NOT (boolean a)
57 // Numeric relational operators.
58 boolean operator EQ (a, b) = a == b;
59 boolean operator GE (a, b) = a >= b;
60 boolean operator GT (a, b) = a > b;
61 boolean operator LE (a, b) = a <= b;
62 boolean operator LT (a, b) = a < b;
63 boolean operator NE (a, b) = a != b;
65 // String relational operators.
66 boolean operator EQ_STRING (string a, string b) = compare_string_3way (&a, &b) == 0;
67 boolean operator GE_STRING (string a, string b) = compare_string_3way (&a, &b) >= 0;
68 boolean operator GT_STRING (string a, string b) = compare_string_3way (&a, &b) > 0;
69 boolean operator LE_STRING (string a, string b) = compare_string_3way (&a, &b) <= 0;
70 boolean operator LT_STRING (string a, string b) = compare_string_3way (&a, &b) < 0;
71 boolean operator NE_STRING (string a, string b) = compare_string_3way (&a, &b) != 0;
74 function ABS (x) = fabs (x);
75 extension function ACOS (x >= -1 && x <= 1) = acos (x);
76 function ASIN (x >= -1 && x <= 1) = asin (x);
77 function ATAN (x) = atan (x);
78 extension function ARCOS (x >= -1 && x <= 1) = acos (x);
79 function ARSIN (x >= -1 && x <= 1) = asin (x);
80 function ARTAN (x) = atan (x);
81 function COS (x) = cos (x);
82 function EXP (x) = check_errno (exp (x));
83 function LG10(x) = check_errno (log10 (x));
84 function LN (x) = check_errno (log (x));
85 function LNGAMMA (x >= 0) = gsl_sf_lngamma (x);
86 function MOD10 (x) = fmod (x, 10);
87 function RND (x) = round_nearest (x, 1, 0);
88 function RND (x, mult != 0) = round_nearest (x, mult, 0);
89 function RND (x, mult != 0, fuzzbits >= 0) = round_nearest (x, mult, fuzzbits);
90 function SIN (x) = sin (x);
91 function SQRT (x >= 0) = sqrt (x);
92 function TAN (x) = check_errno (tan (x));
93 function TRUNC (x) = x >= 0. ? floor (x) : -floor (-x);
95 absorb_miss function MOD (n, d)
98 return n != SYSMIS ? fmod (n, d) : SYSMIS;
100 return n != 0. ? SYSMIS : 0.;
103 // N-ary numeric functions.
104 absorb_miss boolean function ANY (x != SYSMIS, a[n])
109 for (i = 0; i < n; i++)
112 else if (a[i] == SYSMIS)
115 return sysmis ? SYSMIS : 0.;
118 boolean function ANY (string x, string a[n])
122 for (i = 0; i < n; i++)
123 if (!compare_string_3way (&x, &a[i]))
128 function CFVAR.2 (a[n])
130 double mean, variance;
132 moments_of_doubles (a, n, NULL, &mean, &variance, NULL, NULL);
134 if (mean == SYSMIS || mean == 0 || variance == SYSMIS)
137 return sqrt (variance) / mean;
140 function MAX.1 (a[n])
146 for (i = 0; i < n; i++)
147 if (a[i] != SYSMIS && a[i] > max)
152 string function MAX (string a[n])
154 struct substring *max;
158 for (i = 1; i < n; i++)
159 if (compare_string_3way (&a[i], max) > 0)
164 function MEAN.1 (a[n])
167 moments_of_doubles (a, n, NULL, &mean, NULL, NULL, NULL);
171 function MEDIAN.1 (a[n])
173 return median (a, n);
176 function MIN.1 (a[n])
182 for (i = 0; i < n; i++)
183 if (a[i] != SYSMIS && a[i] < min)
188 string function MIN (string a[n])
190 struct substring *min;
194 for (i = 1; i < n; i++)
195 if (compare_string_3way (&a[i], min) < 0)
200 absorb_miss function NMISS (a[n])
203 size_t missing_cnt = 0;
205 for (i = 0; i < n; i++)
206 missing_cnt += a[i] == SYSMIS;
210 absorb_miss function NVALID (a[n])
213 size_t valid_cnt = 0;
215 for (i = 0; i < n; i++)
216 valid_cnt += a[i] != SYSMIS;
220 absorb_miss boolean function RANGE (x != SYSMIS, a[n*2])
225 for (i = 0; i < n; i++)
228 double y = a[2 * i + 1];
229 if (w != SYSMIS && y != SYSMIS)
231 if (w <= x && x <= y)
237 return sysmis ? SYSMIS : 0.;
240 boolean function RANGE (string x, string a[n*2])
244 for (i = 0; i < n; i++)
246 struct substring *w = &a[2 * i];
247 struct substring *y = &a[2 * i + 1];
248 if (compare_string_3way (w, &x) <= 0 && compare_string_3way (&x, y) <= 0)
257 moments_of_doubles (a, n, NULL, NULL, &variance, NULL, NULL);
258 return sqrt (variance);
261 function SUM.1 (a[n])
267 for (i = 0; i < n; i++)
273 function VARIANCE.2 (a[n])
276 moments_of_doubles (a, n, NULL, NULL, &variance, NULL, NULL);
280 // Time construction & extraction functions.
281 function TIME.HMS (h, m, s)
283 if ((h > 0. || m > 0. || s > 0.) && (h < 0. || m < 0. || s < 0.))
285 msg (SW, _("TIME.HMS cannot mix positive and negative arguments."));
289 return H_S * h + MIN_S * m + s;
291 function TIME.DAYS (days) = days * DAY_S;
292 function CTIME.DAYS (time) = time / DAY_S;
293 function CTIME.HOURS (time) = time / H_S;
294 function CTIME.MINUTES (time) = time / MIN_S;
295 function CTIME.SECONDS (time) = time;
297 // Date construction functions.
298 function DATE.DMY (d, m, y) = expr_ymd_to_date (y, m, d);
299 function DATE.MDY (m, d, y) = expr_ymd_to_date (y, m, d);
300 function DATE.MOYR (m, y) = expr_ymd_to_date (y, m, 1);
301 function DATE.QYR (q, y)
303 if (q < 1.0 || q > 4.0 || q != (int) q)
305 msg (SW, _("The first argument to DATE.QYR must be 1, 2, 3, or 4."));
308 return expr_ymd_to_date (y, q * 3 - 2, 1);
310 function DATE.WKYR (w, y) = expr_wkyr_to_date (w, y);
311 function DATE.YRDAY (y, yday) = expr_yrday_to_date (y, yday);
312 function YRMODA (y, m, d) = expr_yrmoda (y, m, d);
314 // Date extraction functions.
315 function XDATE.TDAY (date) = floor (date / DAY_S);
316 function XDATE.HOUR (date) = fmod (floor (date / H_S), DAY_H);
317 function XDATE.MINUTE (date) = fmod (floor (date / H_MIN), H_MIN);
318 function XDATE.SECOND (date) = fmod (date, MIN_S);
319 function XDATE.DATE (date) = floor (date / DAY_S) * DAY_S;
320 function XDATE.TIME (date) = fmod (date, DAY_S);
322 function XDATE.JDAY (date >= DAY_S) = calendar_offset_to_yday (date / DAY_S);
323 function XDATE.MDAY (date >= DAY_S) = calendar_offset_to_mday (date / DAY_S);
324 function XDATE.MONTH (date >= DAY_S)
325 = calendar_offset_to_month (date / DAY_S);
326 function XDATE.QUARTER (date >= DAY_S)
327 = (calendar_offset_to_month (date / DAY_S) - 1) / 3 + 1;
328 function XDATE.WEEK (date >= DAY_S)
329 = (calendar_offset_to_yday (date / DAY_S) - 1) / 7 + 1;
330 function XDATE.WKDAY (date >= DAY_S) = calendar_offset_to_wday (date / DAY_S);
331 function XDATE.YEAR (date >= DAY_S) = calendar_offset_to_year (date / DAY_S);
333 // Date arithmetic functions.
334 no_abbrev function DATEDIFF (date2 >= DAY_S, date1 >= DAY_S, string unit)
335 = expr_date_difference (date1, date2, unit);
336 no_abbrev function DATESUM (date, quantity, string unit)
337 = expr_date_sum (date, quantity, unit, ss_cstr ("closest"));
338 no_abbrev function DATESUM (date, quantity, string unit, string method)
339 = expr_date_sum (date, quantity, unit, method);
343 string function CONCAT (string a[n])
346 struct substring dst;
349 dst = alloc_string (e, MAX_STRING);
351 for (i = 0; i < n; i++)
353 struct substring *src = &a[i];
356 copy_len = src->length;
357 if (dst.length + copy_len > MAX_STRING)
358 copy_len = MAX_STRING - dst.length;
359 memcpy (&dst.string[dst.length], src->string, copy_len);
360 dst.length += copy_len;
366 function INDEX (string haystack, string needle)
368 if (needle.length == 0)
372 int limit = haystack.length - needle.length + 1;
374 for (i = 1; i <= limit; i++)
375 if (!memcmp (&haystack.string[i - 1], needle.string, needle.length))
381 function INDEX (string haystack, string needles, needle_len_d)
383 if (needle_len_d <= INT_MIN || needle_len_d >= INT_MAX
384 || (int) needle_len_d != needle_len_d
385 || needles.length == 0)
389 int needle_len = needle_len_d;
390 if (needle_len < 0 || needle_len > needles.length
391 || needles.length % needle_len != 0)
395 int limit = haystack.length - needle_len + 1;
397 for (i = 1; i <= limit; i++)
398 for (j = 0; j < needles.length; j += needle_len)
399 if (!memcmp (&haystack.string[i - 1], &needles.string[j],
407 function RINDEX (string haystack, string needle)
409 if (needle.length == 0)
413 int limit = haystack.length - needle.length + 1;
415 for (i = limit; i >= 1; i--)
416 if (!memcmp (&haystack.string[i - 1], needle.string, needle.length))
422 function RINDEX (string haystack, string needles, needle_len_d)
424 if (needle_len_d <= 0 || needle_len_d >= INT_MAX
425 || (int) needle_len_d != needle_len_d
426 || needles.length == 0)
430 int needle_len = needle_len_d;
431 if (needle_len < 0 || needle_len > needles.length
432 || needles.length % needle_len != 0)
436 int limit = haystack.length - needle_len + 1;
438 for (i = limit; i >= 1; i--)
439 for (j = 0; j < needles.length; j += needle_len)
440 if (!memcmp (&haystack.string[i - 1],
441 &needles.string[j], needle_len))
448 function LENGTH (string s)
453 string function LOWER (string s)
457 for (i = 0; i < s.length; i++)
458 s.string[i] = tolower ((unsigned char) s.string[i]);
462 function MBLEN.BYTE (string s, idx)
464 if (idx < 0 || idx >= s.length || (int) idx != idx)
470 string function UPCASE (string s)
474 for (i = 0; i < s.length; i++)
475 s.string[i] = toupper ((unsigned char) s.string[i]);
479 absorb_miss string function LPAD (string s, n)
482 if (n < 0 || n > MAX_STRING || (int) n != n)
484 else if (s.length >= n)
488 struct substring t = alloc_string (e, n);
489 memset (t.string, ' ', n - s.length);
490 memcpy (&t.string[(int) n - s.length], s.string, s.length);
495 absorb_miss string function LPAD (string s, n, string c)
498 if (n < 0 || n > MAX_STRING || (int) n != n || c.length != 1)
500 else if (s.length >= n)
504 struct substring t = alloc_string (e, n);
505 memset (t.string, c.string[0], n - s.length);
506 memcpy (&t.string[(int) n - s.length], s.string, s.length);
511 string function REPLACE (string haystack, string needle, string replacement)
513 = replace_string (e, haystack, needle, replacement, DBL_MAX);
515 absorb_miss string function REPLACE (string haystack, string needle,
516 string replacement, n)
518 = replace_string (e, haystack, needle, replacement, n);
520 absorb_miss string function RPAD (string s, n)
523 if (n < 0 || n > MAX_STRING || (int) n != n)
525 else if (s.length >= n)
529 struct substring t = alloc_string (e, n);
530 memcpy (t.string, s.string, s.length);
531 memset (&t.string[s.length], ' ', n - s.length);
536 absorb_miss string function RPAD (string s, n, string c)
539 if (n < 0 || n > MAX_STRING || (int) n != n || c.length != 1)
541 else if (s.length >= n)
545 struct substring t = alloc_string (e, n);
546 memcpy (t.string, s.string, s.length);
547 memset (&t.string[s.length], c.string[0], n - s.length);
552 string function LTRIM (string s)
554 while (s.length > 0 && s.string[0] == ' ')
562 string function LTRIM (string s, string c)
566 while (s.length > 0 && s.string[0] == c.string[0])
577 string function RTRIM (string s)
579 while (s.length > 0 && s.string[s.length - 1] == ' ')
584 string function RTRIM (string s, string c)
588 while (s.length > 0 && s.string[s.length - 1] == c.string[0])
596 function NUMBER (string s, ni_format f)
603 error = data_in (s, C_ENCODING, f->type, &out, 0, NULL);
605 data_in_imply_decimals (s, C_ENCODING, f->type, f->d, &out);
608 msg (SE, "Cannot parse `%.*s' as format %s: %s",
609 (int) s.length, s.string, fmt_name (f->type), error);
615 absorb_miss string function STRING (x, no_format f)
619 struct substring dst;
624 assert (!fmt_is_string (f->type));
625 s = data_out (&v, C_ENCODING, f);
626 dst = alloc_string (e, strlen (s));
627 strcpy (dst.string, s);
632 absorb_miss string function STRUNC (string s, n)
634 if (n < 1 || n == SYSMIS)
639 while (s.length > 0 && s.string[s.length - 1] == ' ')
644 absorb_miss string function SUBSTR (string s, ofs)
647 if (ofs >= 1 && ofs <= s.length && (int) ofs == ofs)
648 return copy_string (e, &s.string[(int) ofs - 1], s.length - ofs + 1);
653 absorb_miss string function SUBSTR (string s, ofs, cnt)
656 if (ofs >= 1 && ofs <= s.length && (int) ofs == ofs
657 && cnt >= 1 && cnt <= INT_MAX && (int) cnt == cnt)
659 int cnt_max = s.length - (int) ofs + 1;
660 return copy_string (e, &s.string[(int) ofs - 1],
661 cnt <= cnt_max ? cnt : cnt_max);
667 absorb_miss no_opt no_abbrev string function VALUELABEL (var v)
671 const char *label = var_lookup_value_label (v, case_data (c, v));
673 return copy_string (e, label, strlen (label));
679 operator SQUARE (x) = x * x;
680 boolean operator NUM_TO_BOOLEAN (x, string op_name)
682 if (x == 0. || x == 1. || x == SYSMIS)
685 if (!ss_is_empty (op_name))
686 msg (SE, _("An operand of the %.*s operator was found to have a value "
687 "other than 0 (false), 1 (true), or the system-missing "
688 "value. The result was forced to 0."),
689 (int) op_name.length, op_name.string);
691 msg (SE, _("A logical expression was found to have a value other than 0 "
692 "(false), 1 (true), or the system-missing value. The result "
693 "was forced to 0."));
697 operator BOOLEAN_TO_NUM (boolean x) = x;
699 // Beta distribution.
700 function PDF.BETA (x >= 0 && x <= 1, a > 0, b > 0)
701 = gsl_ran_beta_pdf (x, a, b);
702 function CDF.BETA (x >= 0 && x <= 1, a > 0, b > 0) = gsl_cdf_beta_P (x, a, b);
703 function IDF.BETA (P >= 0 && P <= 1, a > 0, b > 0)
704 = gsl_cdf_beta_Pinv (P, a, b);
705 no_opt function RV.BETA (a > 0, b > 0) = gsl_ran_beta (get_rng (), a, b);
706 function NCDF.BETA (x >= 0, a > 0, b > 0, lambda > 0)
707 = ncdf_beta (x, a, b, lambda);
708 function NPDF.BETA (x >= 0, a > 0, b > 0, lambda > 0)
709 = npdf_beta (x, a, b, lambda);
711 // Bivariate normal distribution.
712 function CDF.BVNOR (x0, x1, r >= -1 && r <= 1) = cdf_bvnor (x0, x1, r);
713 function PDF.BVNOR (x0, x1, r >= -1 && r <= 1)
714 = gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
716 // Cauchy distribution.
717 function CDF.CAUCHY (x, a, b > 0) = gsl_cdf_cauchy_P ((x - a) / b, 1);
718 function IDF.CAUCHY (P > 0 && P < 1, a, b > 0)
719 = a + b * gsl_cdf_cauchy_Pinv (P, 1);
720 function PDF.CAUCHY (x, a, b > 0) = gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
721 no_opt function RV.CAUCHY (a, b > 0) = a + b * gsl_ran_cauchy (get_rng (), 1);
723 // Chi-square distribution.
724 function CDF.CHISQ (x >= 0, df > 0) = gsl_cdf_chisq_P (x, df);
725 function IDF.CHISQ (P >= 0 && P < 1, df > 0) = gsl_cdf_chisq_Pinv (P, df);
726 function PDF.CHISQ (x >= 0, df > 0) = gsl_ran_chisq_pdf (x, df);
727 no_opt function RV.CHISQ (df > 0) = gsl_ran_chisq (get_rng (), df);
728 function NCDF.CHISQ (x >= 0, df > 0, c) = unimplemented;
729 function NPDF.CHISQ (x >= 0, df > 0, c) = unimplemented;
730 function SIG.CHISQ (x >= 0, df > 0) = gsl_cdf_chisq_Q (x, df);
732 // Exponential distribution.
733 function CDF.EXP (x >= 0, a > 0) = gsl_cdf_exponential_P (x, 1. / a);
734 function IDF.EXP (P >= 0 && P < 1, a > 0)
735 = gsl_cdf_exponential_Pinv (P, 1. / a);
736 function PDF.EXP (x >= 0, a > 0) = gsl_ran_exponential_pdf (x, 1. / a);
737 no_opt function RV.EXP (a > 0) = gsl_ran_exponential (get_rng (), 1. / a);
739 // Exponential power distribution.
740 extension function PDF.XPOWER (x, a > 0, b >= 0)
741 = gsl_ran_exppow_pdf (x, a, b);
742 no_opt extension function RV.XPOWER (a > 0, b >= 0)
743 = gsl_ran_exppow (get_rng (), a, b);
746 function CDF.F (x >= 0, df1 > 0, df2 > 0) = gsl_cdf_fdist_P (x, df1, df2);
747 function IDF.F (P >= 0 && P < 1, df1 > 0, df2 > 0) = idf_fdist (P, df1, df2);
748 function PDF.F (x >= 0, df1 > 0, df2 > 0) = gsl_ran_fdist_pdf (x, df1, df2);
749 no_opt function RV.F (df1 > 0, df2 > 0) = gsl_ran_fdist (get_rng (), df1, df2);
750 function NCDF.F (x >= 0, df1 > 0, df2 > 0, lambda >= 0) = unimplemented;
751 function NPDF.F (x >= 0, df1 > 0, df2 > 0, lmabda >= 0) = unimplemented;
752 function SIG.F (x >= 0, df1 > 0, df2 > 0) = gsl_cdf_fdist_Q (x, df1, df2);
754 // Gamma distribution.
755 function CDF.GAMMA (x >= 0, a > 0, b > 0) = gsl_cdf_gamma_P (x, a, 1. / b);
756 function IDF.GAMMA (P >= 0 && P <= 1, a > 0, b > 0)
757 = gsl_cdf_gamma_Pinv (P, a, 1. / b);
758 function PDF.GAMMA (x >= 0, a > 0, b > 0) = gsl_ran_gamma_pdf (x, a, 1. / b);
759 no_opt function RV.GAMMA (a > 0, b > 0)
760 = gsl_ran_gamma (get_rng (), a, 1. / b);
762 // Half-normal distribution.
763 function CDF.HALFNRM (x, a, b > 0) = unimplemented;
764 function IDF.HALFNRM (P > 0 && P < 1, a, b > 0) = unimplemented;
765 function PDF.HALFNRM (x, a, b > 0) = unimplemented;
766 no_opt function RV.HALFNRM (a, b > 0) = unimplemented;
768 // Inverse Gaussian distribution.
769 function CDF.IGAUSS (x > 0, a > 0, b > 0) = unimplemented;
770 function IDF.IGAUSS (P >= 0 && P < 1, a > 0, b > 0) = unimplemented;
771 function PDF.IGAUSS (x > 0, a > 0, b > 0) = unimplemented;
772 no_opt function RV.IGAUSS (a > 0, b > 0) = unimplemented;
774 // Landau distribution.
775 extension function PDF.LANDAU (x) = gsl_ran_landau_pdf (x);
776 no_opt extension function RV.LANDAU () = gsl_ran_landau (get_rng ());
778 // Laplace distribution.
779 function CDF.LAPLACE (x, a, b > 0) = gsl_cdf_laplace_P ((x - a) / b, 1);
780 function IDF.LAPLACE (P > 0 && P < 1, a, b > 0)
781 = a + b * gsl_cdf_laplace_Pinv (P, 1);
782 function PDF.LAPLACE (x, a, b > 0) = gsl_ran_laplace_pdf ((x - a) / b, 1) / b;
783 no_opt function RV.LAPLACE (a, b > 0)
784 = a + b * gsl_ran_laplace (get_rng (), 1);
786 // Levy alpha-stable distribution.
787 no_opt extension function RV.LEVY (c, alpha > 0 && alpha <= 2)
788 = gsl_ran_levy (get_rng (), c, alpha);
790 // Levy skew alpha-stable distribution.
791 no_opt extension function RV.LVSKEW (c, alpha > 0 && alpha <= 2,
792 beta >= -1 && beta <= 1)
793 = gsl_ran_levy_skew (get_rng (), c, alpha, beta);
795 // Logistic distribution.
796 function CDF.LOGISTIC (x, a, b > 0) = gsl_cdf_logistic_P ((x - a) / b, 1);
797 function IDF.LOGISTIC (P > 0 && P < 1, a, b > 0)
798 = a + b * gsl_cdf_logistic_Pinv (P, 1);
799 function PDF.LOGISTIC (x, a, b > 0)
800 = gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
801 no_opt function RV.LOGISTIC (a, b > 0)
802 = a + b * gsl_ran_logistic (get_rng (), 1);
804 // Lognormal distribution.
805 function CDF.LNORMAL (x >= 0, m > 0, s > 0)
806 = gsl_cdf_lognormal_P (x, log (m), s);
807 function IDF.LNORMAL (P >= 0 && P < 1, m > 0, s > 0)
808 = gsl_cdf_lognormal_Pinv (P, log (m), s);
809 function PDF.LNORMAL (x >= 0, m > 0, s > 0)
810 = gsl_ran_lognormal_pdf (x, log (m), s);
811 no_opt function RV.LNORMAL (m > 0, s > 0)
812 = gsl_ran_lognormal (get_rng (), log (m), s);
814 // Normal distribution.
815 function CDF.NORMAL (x, u, s > 0) = gsl_cdf_gaussian_P (x - u, s);
816 function IDF.NORMAL (P > 0 && P < 1, u, s > 0)
817 = u + gsl_cdf_gaussian_Pinv (P, s);
818 function PDF.NORMAL (x, u, s > 0) = gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
819 no_opt function RV.NORMAL (u, s > 0) = u + gsl_ran_gaussian (get_rng (), s);
820 function CDFNORM (x) = gsl_cdf_ugaussian_P (x);
821 function PROBIT (P > 0 && P < 1) = gsl_cdf_ugaussian_Pinv (P);
822 no_opt function NORMAL (s > 0) = gsl_ran_gaussian (get_rng (), s);
824 // Normal tail distribution.
825 function PDF.NTAIL (x, a > 0, sigma > 0)
826 = gsl_ran_gaussian_tail_pdf (x, a, sigma);
827 no_opt function RV.NTAIL (a > 0, sigma > 0)
828 = gsl_ran_gaussian_tail (get_rng (), a, sigma);
830 // Pareto distribution.
831 function CDF.PARETO (x >= a, a > 0, b > 0) = gsl_cdf_pareto_P (x, b, a);
832 function IDF.PARETO (P >= 0 && P < 1, a > 0, b > 0)
833 = gsl_cdf_pareto_Pinv (P, b, a);
834 function PDF.PARETO (x >= a, a > 0, b > 0) = gsl_ran_pareto_pdf (x, b, a);
835 no_opt function RV.PARETO (a > 0, b > 0) = gsl_ran_pareto (get_rng (), b, a);
837 // Rayleigh distribution.
838 extension function CDF.RAYLEIGH (x, sigma > 0) = gsl_cdf_rayleigh_P (x, sigma);
839 extension function IDF.RAYLEIGH (P >= 0 && P <= 1, sigma > 0)
840 = gsl_cdf_rayleigh_Pinv (P, sigma);
841 extension function PDF.RAYLEIGH (x, sigma > 0)
842 = gsl_ran_rayleigh_pdf (x, sigma);
843 no_opt extension function RV.RAYLEIGH (sigma > 0)
844 = gsl_ran_rayleigh (get_rng (), sigma);
846 // Rayleigh tail distribution.
847 extension function PDF.RTAIL (x, a, sigma)
848 = gsl_ran_rayleigh_tail_pdf (x, a, sigma);
849 no_opt extension function RV.RTAIL (a, sigma)
850 = gsl_ran_rayleigh_tail (get_rng (), a, sigma);
852 // Studentized maximum modulus distribution.
853 function CDF.SMOD (x > 0, a >= 1, b >= 1) = unimplemented;
854 function IDF.SMOD (P >= 0 && P < 1, a >= 1, b >= 1) = unimplemented;
856 // Studentized range distribution.
857 function CDF.SRANGE (x > 0, a >= 1, b >= 1) = unimplemented;
858 function IDF.SRANGE (P >= 0 && P < 1, a >= 1, b >= 1) = unimplemented;
860 // Student t distribution.
861 function CDF.T (x, df > 0) = gsl_cdf_tdist_P (x, df);
862 function IDF.T (P > 0 && P < 1, df > 0) = gsl_cdf_tdist_Pinv (P, df);
863 function PDF.T (x, df > 0) = gsl_ran_tdist_pdf (x, df);
864 no_opt function RV.T (df > 0) = gsl_ran_tdist (get_rng (), df);
865 function NCDF.T (x, df > 0, nc) = unimplemented;
866 function NPDF.T (x, df > 0, nc) = unimplemented;
868 // Type-1 Gumbel distribution.
869 extension function CDF.T1G (x, a, b) = gsl_cdf_gumbel1_P (x, a, b);
870 extension function IDF.T1G (P >= 0 && P <= 1, a, b)
871 = gsl_cdf_gumbel1_P (P, a, b);
872 extension function PDF.T1G (x, a, b) = gsl_ran_gumbel1_pdf (x, a, b);
873 no_opt extension function RV.T1G (a, b) = gsl_ran_gumbel1 (get_rng (), a, b);
875 // Type-2 Gumbel distribution.
876 extension function CDF.T2G (x, a, b) = gsl_cdf_gumbel2_P (x, a, b);
877 extension function IDF.T2G (P >= 0 && P <= 1, a, b)
878 = gsl_cdf_gumbel2_P (P, a, b);
879 extension function PDF.T2G (x, a, b) = gsl_ran_gumbel2_pdf (x, a, b);
880 no_opt extension function RV.T2G (a, b) = gsl_ran_gumbel2 (get_rng (), a, b);
882 // Uniform distribution.
883 function CDF.UNIFORM (x <= b, a <= x, b) = gsl_cdf_flat_P (x, a, b);
884 function IDF.UNIFORM (P >= 0 && P <= 1, a <= b, b)
885 = gsl_cdf_flat_Pinv (P, a, b);
886 function PDF.UNIFORM (x <= b, a <= x, b) = gsl_ran_flat_pdf (x, a, b);
887 no_opt function RV.UNIFORM (a <= b, b) = gsl_ran_flat (get_rng (), a, b);
888 no_opt function UNIFORM (b >= 0) = gsl_ran_flat (get_rng (), 0, b);
890 // Weibull distribution.
891 function CDF.WEIBULL (x >= 0, a > 0, b > 0) = gsl_cdf_weibull_P (x, a, b);
892 function IDF.WEIBULL (P >= 0 && P < 1, a > 0, b > 0)
893 = gsl_cdf_weibull_Pinv (P, a, b);
894 function PDF.WEIBULL (x >= 0, a > 0, b > 0) = gsl_ran_weibull_pdf (x, a, b);
895 no_opt function RV.WEIBULL (a > 0, b > 0) = gsl_ran_weibull (get_rng (), a, b);
897 // Bernoulli distribution.
898 function CDF.BERNOULLI (k == 0 || k == 1, p >= 0 && p <= 1)
900 function PDF.BERNOULLI (k == 0 || k == 1, p >= 0 && p <= 1)
901 = gsl_ran_bernoulli_pdf (k, p);
902 no_opt function RV.BERNOULLI (p >= 0 && p <= 1)
903 = gsl_ran_bernoulli (get_rng (), p);
905 // Binomial distribution.
906 function CDF.BINOM (k, n > 0 && n == floor (n), p >= 0 && p <= 1)
907 = gsl_cdf_binomial_P (k, p, n);
908 function PDF.BINOM (k >= 0 && k == floor (k) && k <= n,
909 n > 0 && n == floor (n),
911 = gsl_ran_binomial_pdf (k, p, n);
912 no_opt function RV.BINOM (p > 0 && p == floor (p), n >= 0 && n <= 1)
913 = gsl_ran_binomial (get_rng (), p, n);
915 // Geometric distribution.
916 function CDF.GEOM (k >= 1 && k == floor (k), p >= 0 && p <= 1)
917 = gsl_cdf_geometric_P (k, p);
918 function PDF.GEOM (k >= 1 && k == floor (k),
920 = gsl_ran_geometric_pdf (k, p);
921 no_opt function RV.GEOM (p >= 0 && p <= 1) = gsl_ran_geometric (get_rng (), p);
923 // Hypergeometric distribution.
924 function CDF.HYPER (k >= 0 && k == floor (k) && k <= c,
925 a > 0 && a == floor (a),
926 b > 0 && b == floor (b) && b <= a,
927 c > 0 && c == floor (c) && c <= a)
928 = gsl_cdf_hypergeometric_P (k, c, a - c, b);
929 function PDF.HYPER (k >= 0 && k == floor (k) && k <= c,
930 a > 0 && a == floor (a),
931 b > 0 && b == floor (b) && b <= a,
932 c > 0 && c == floor (c) && c <= a)
933 = gsl_ran_hypergeometric_pdf (k, c, a - c, b);
934 no_opt function RV.HYPER (a > 0 && a == floor (a),
935 b > 0 && b == floor (b) && b <= a,
936 c > 0 && c == floor (c) && c <= a)
937 = gsl_ran_hypergeometric (get_rng (), c, a - c, b);
939 // Logarithmic distribution.
940 extension function PDF.LOG (k >= 1, p > 0 && p <= 1)
941 = gsl_ran_logarithmic_pdf (k, p);
942 no_opt extension function RV.LOG (p > 0 && p <= 1)
943 = gsl_ran_logarithmic (get_rng (), p);
945 // Negative binomial distribution.
946 function CDF.NEGBIN (k >= 1, n == floor (n), p > 0 && p <= 1)
947 = gsl_cdf_negative_binomial_P (k, p, n);
948 function PDF.NEGBIN (k >= 1, n == floor (n), p > 0 && p <= 1)
949 = gsl_ran_negative_binomial_pdf (k, p, n);
950 no_opt function RV.NEGBIN (n == floor (n), p > 0 && p <= 1)
951 = gsl_ran_negative_binomial (get_rng (), p, n);
953 // Poisson distribution.
954 function CDF.POISSON (k >= 0 && k == floor (k), mu > 0)
955 = gsl_cdf_poisson_P (k, mu);
956 function PDF.POISSON (k >= 0 && k == floor (k), mu > 0)
957 = gsl_ran_poisson_pdf (k, mu);
958 no_opt function RV.POISSON (mu > 0) = gsl_ran_poisson (get_rng (), mu);
961 absorb_miss boolean function MISSING (x) = x == SYSMIS || !finite (x);
962 absorb_miss boolean function SYSMIS (x) = x == SYSMIS || !finite (x);
963 no_opt boolean function SYSMIS (num_var v)
966 return case_num (c, v) == SYSMIS;
968 no_opt boolean function VALUE (num_var v)
971 return case_num (c, v);
974 no_opt operator VEC_ELEM_NUM (idx)
978 if (idx >= 1 && idx <= vector_get_var_cnt (v))
980 const struct variable *var = vector_get_var (v, (size_t) idx - 1);
981 double value = case_num (c, var);
982 return !var_is_num_missing (var, value, MV_USER) ? value : SYSMIS;
987 msg (SE, _("SYSMIS is not a valid index value for vector "
988 "%s. The result will be set to SYSMIS."),
989 vector_get_name (v));
991 msg (SE, _("%g is not a valid index value for vector %s. "
992 "The result will be set to SYSMIS."),
993 idx, vector_get_name (v));
998 absorb_miss no_opt string operator VEC_ELEM_STR (idx)
1003 if (idx >= 1 && idx <= vector_get_var_cnt (v))
1005 struct variable *var = vector_get_var (v, (size_t) idx - 1);
1006 return copy_string (e, CHAR_CAST_BUG (char *, case_str (c, var)),
1007 var_get_width (var));
1012 msg (SE, _("SYSMIS is not a valid index value for vector "
1013 "%s. The result will be set to the empty string."),
1014 vector_get_name (v));
1016 msg (SE, _("%g is not a valid index value for vector %s. "
1017 "The result will be set to the empty string."),
1018 idx, vector_get_name (v));
1019 return empty_string;
1025 no_opt operator NUM_VAR ()
1029 double d = case_num (c, v);
1030 return !var_is_num_missing (v, d, MV_USER) ? d : SYSMIS;
1033 no_opt string operator STR_VAR ()
1038 struct substring s = alloc_string (e, var_get_width (v));
1039 memcpy (s.string, case_str (c, v), var_get_width (v));
1043 no_opt perm_only function LAG (num_var v, pos_int n_before)
1046 const struct ccase *c = lagged_case (ds, n_before);
1049 double x = case_num (c, v);
1050 return !var_is_num_missing (v, x, MV_USER) ? x : SYSMIS;
1056 no_opt perm_only function LAG (num_var v)
1059 const struct ccase *c = lagged_case (ds, 1);
1062 double x = case_num (c, v);
1063 return !var_is_num_missing (v, x, MV_USER) ? x : SYSMIS;
1069 no_opt perm_only string function LAG (str_var v, pos_int n_before)
1073 const struct ccase *c = lagged_case (ds, n_before);
1075 return copy_string (e, CHAR_CAST_BUG (char *, case_str (c, v)),
1078 return empty_string;
1081 no_opt perm_only string function LAG (str_var v)
1085 const struct ccase *c = lagged_case (ds, 1);
1087 return copy_string (e, CHAR_CAST_BUG (char *, case_str (c, v)),
1090 return empty_string;
1093 no_opt operator NUM_SYS ()
1097 return case_num (c, v) == SYSMIS;
1100 no_opt operator NUM_VAL ()
1104 return case_num (c, v);
1107 no_opt operator CASENUM ()