3 // PSPP - computes sample statistics.
4 // Copyright (C) 2005, 2006 Free Software Foundation, Inc.
5 // Written by Ben Pfaff <blp@gnu.org>.
7 // This program is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU General Public License as
9 // published by the Free Software Foundation; either version 2 of the
10 // License, or (at your option) any later version.
12 // This program is distributed in the hope that it will be useful, but
13 // WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 // General Public License for more details.
17 // You should have received a copy of the GNU General Public License
18 // along with this program; if not, write to the Free Software
19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 // 02110-1301, USA. */
22 operator NEG (x) = -x;
24 operator ADD (a, b) = a + b;
25 operator SUB (a, b) = a - b;
27 absorb_miss operator MUL (a, b)
28 = (a == 0. || b == 0. ? 0.
29 : a == SYSMIS || b == SYSMIS ? SYSMIS
32 absorb_miss operator DIV (a, b)
34 : a == SYSMIS || b == SYSMIS ? SYSMIS
37 absorb_miss operator POW (a, b)
38 = (a == SYSMIS ? (b == 0. ? 1. : a)
39 : b == SYSMIS ? (a == 0. ? 0. : SYSMIS)
40 : a == 0. && b <= 0. ? SYSMIS
43 absorb_miss boolean operator AND (boolean a, boolean b)
46 : b == SYSMIS ? SYSMIS
49 absorb_miss boolean operator OR (boolean a, boolean b)
52 : b == SYSMIS ? SYSMIS
55 boolean operator NOT (boolean a)
60 // Numeric relational operators.
61 boolean operator EQ (a, b) = a == b;
62 boolean operator GE (a, b) = a >= b;
63 boolean operator GT (a, b) = a > b;
64 boolean operator LE (a, b) = a <= b;
65 boolean operator LT (a, b) = a < b;
66 boolean operator NE (a, b) = a != b;
68 // String relational operators.
69 boolean operator EQ_STRING (string a, string b) = compare_string (&a, &b) == 0;
70 boolean operator GE_STRING (string a, string b) = compare_string (&a, &b) >= 0;
71 boolean operator GT_STRING (string a, string b) = compare_string (&a, &b) > 0;
72 boolean operator LE_STRING (string a, string b) = compare_string (&a, &b) <= 0;
73 boolean operator LT_STRING (string a, string b) = compare_string (&a, &b) < 0;
74 boolean operator NE_STRING (string a, string b) = compare_string (&a, &b) != 0;
77 function ABS (x) = fabs (x);
78 extension function ACOS (x >= -1 && x <= 1) = acos (x);
79 function ASIN (x >= -1 && x <= 1) = asin (x);
80 function ATAN (x) = atan (x);
81 extension function ARCOS (x >= -1 && x <= 1) = acos (x);
82 function ARSIN (x >= -1 && x <= 1) = asin (x);
83 function ARTAN (x) = atan (x);
84 function COS (x) = cos (x);
85 function EXP (x) = check_errno (exp (x));
86 function LG10(x) = check_errno (log10 (x));
87 function LN (x) = check_errno (log (x));
88 function LNGAMMA (x >= 0) = gsl_sf_lngamma (x);
89 function MOD10 (x) = fmod (x, 10);
90 function RND (x) = x >= 0. ? floor (x + .5) : -floor (-x + .5);
91 function SIN (x) = sin (x);
92 function SQRT (x >= 0) = sqrt (x);
93 function TAN (x) = check_errno (tan (x));
94 function TRUNC (x) = x >= 0. ? floor (x) : -floor (-x);
96 absorb_miss function MOD (n, d)
99 return n != SYSMIS ? fmod (n, d) : SYSMIS;
101 return n != 0. ? SYSMIS : 0.;
104 // N-ary numeric functions.
105 absorb_miss boolean function ANY (x != SYSMIS, a[n])
110 for (i = 0; i < n; i++)
113 else if (a[i] == SYSMIS)
116 return sysmis ? SYSMIS : 0.;
119 boolean function ANY (string x, string a[n])
123 for (i = 0; i < n; i++)
124 if (!compare_string (&x, &a[i]))
129 function CFVAR.2 (a[n])
131 double mean, variance;
133 moments_of_doubles (a, n, NULL, &mean, &variance, NULL, NULL);
135 if (mean == SYSMIS || mean == 0 || variance == SYSMIS)
138 return sqrt (variance) / mean;
141 function MAX.1 (a[n])
147 for (i = 0; i < n; i++)
148 if (a[i] != SYSMIS && a[i] > max)
153 string function MAX (string a[n])
155 struct substring *max;
159 for (i = 1; i < n; i++)
160 if (compare_string (&a[i], max) > 0)
165 function MEAN.1 (a[n])
168 moments_of_doubles (a, n, NULL, &mean, NULL, NULL, NULL);
172 function MIN.1 (a[n])
178 for (i = 0; i < n; i++)
179 if (a[i] != SYSMIS && a[i] < min)
184 string function MIN (string a[n])
186 struct substring *min;
190 for (i = 1; i < n; i++)
191 if (compare_string (&a[i], min) < 0)
196 absorb_miss function NMISS (a[n])
199 size_t missing_cnt = 0;
201 for (i = 0; i < n; i++)
202 missing_cnt += a[i] == SYSMIS;
206 absorb_miss function NVALID (a[n])
209 size_t valid_cnt = 0;
211 for (i = 0; i < n; i++)
212 valid_cnt += a[i] != SYSMIS;
216 absorb_miss boolean function RANGE (x != SYSMIS, a[n*2])
221 for (i = 0; i < n; i++)
224 double y = a[2 * i + 1];
225 if (w != SYSMIS && y != SYSMIS)
227 if (w <= x && x <= y)
233 return sysmis ? SYSMIS : 0.;
236 boolean function RANGE (string x, string a[n*2])
240 for (i = 0; i < n; i++)
242 struct substring *w = &a[2 * i];
243 struct substring *y = &a[2 * i + 1];
244 if (compare_string (w, &x) <= 0 && compare_string (&x, y) <= 0)
253 moments_of_doubles (a, n, NULL, NULL, &variance, NULL, NULL);
254 return sqrt (variance);
257 function SUM.1 (a[n])
263 for (i = 0; i < n; i++)
269 function VARIANCE.2 (a[n])
272 moments_of_doubles (a, n, NULL, NULL, &variance, NULL, NULL);
276 // Time construction & extraction functions.
277 function TIME.HMS (h, m, s)
279 if ((h > 0. || m > 0. || s > 0.) && (h < 0. || m < 0. || s < 0.))
281 msg (SW, _("TIME.HMS cannot mix positive and negative arguments."));
285 return H_S * h + MIN_S * m + s;
287 function TIME.DAYS (days) = days * DAY_S;
288 function CTIME.DAYS (time) = time / DAY_S;
289 function CTIME.HOURS (time) = time / H_S;
290 function CTIME.MINUTES (time) = time / MIN_S;
291 function CTIME.SECONDS (time) = time;
293 // Date construction functions.
294 function DATE.DMY (d, m, y) = expr_ymd_to_date (y, m, d);
295 function DATE.MDY (m, d, y) = expr_ymd_to_date (y, m, d);
296 function DATE.MOYR (m, y) = expr_ymd_to_date (y, m, 1);
297 function DATE.QYR (q, y) = expr_ymd_to_date (y, q * 3 - 2, 1);
298 function DATE.WKYR (w, y) = expr_wkyr_to_date (w, y);
299 function DATE.YRDAY (y, yday) = expr_yrday_to_date (y, yday);
300 function YRMODA (y, m, d) = expr_yrmoda (y, m, d);
302 // Date extraction functions.
303 function XDATE.TDAY (date) = floor (date / DAY_S);
304 function XDATE.HOUR (date) = fmod (floor (date / H_S), DAY_H);
305 function XDATE.MINUTE (date) = fmod (floor (date / H_MIN), H_MIN);
306 function XDATE.SECOND (date) = fmod (date, MIN_S);
307 function XDATE.DATE (date) = floor (date / DAY_S) * DAY_S;
308 function XDATE.TIME (date) = fmod (date, DAY_S);
310 function XDATE.JDAY (date >= DAY_S) = calendar_offset_to_yday (date / DAY_S);
311 function XDATE.MDAY (date >= DAY_S) = calendar_offset_to_mday (date / DAY_S);
312 function XDATE.MONTH (date >= DAY_S)
313 = calendar_offset_to_month (date / DAY_S);
314 function XDATE.QUARTER (date >= DAY_S)
315 = (calendar_offset_to_month (date / DAY_S) - 1) / 3 + 1;
316 function XDATE.WEEK (date >= DAY_S)
317 = (calendar_offset_to_yday (date / DAY_S) - 1) / 7 + 1;
318 function XDATE.WKDAY (date >= DAY_S) = calendar_offset_to_wday (date / DAY_S);
319 function XDATE.YEAR (date >= DAY_S) = calendar_offset_to_year (date / DAY_S);
321 // Date arithmetic functions.
322 no_abbrev function DATEDIFF (date1 >= DAY_S, date2 >= DAY_S, string unit)
323 = expr_date_difference (date1, date2, unit);
324 no_abbrev function DATESUM (date, quantity, string unit)
325 = expr_date_sum (date, quantity, unit, ss_cstr ("closest"));
326 no_abbrev function DATESUM (date, quantity, string unit, string method)
327 = expr_date_sum (date, quantity, unit, method);
331 string function CONCAT (string a[n])
334 struct substring dst;
337 dst = alloc_string (e, MAX_STRING);
339 for (i = 0; i < n; i++)
341 struct substring *src = &a[i];
344 copy_len = src->length;
345 if (dst.length + copy_len > MAX_STRING)
346 copy_len = MAX_STRING - dst.length;
347 memcpy (&dst.string[dst.length], src->string, copy_len);
348 dst.length += copy_len;
354 function INDEX (string haystack, string needle)
356 if (needle.length == 0)
360 int limit = haystack.length - needle.length + 1;
362 for (i = 1; i <= limit; i++)
363 if (!memcmp (&haystack.string[i - 1], needle.string, needle.length))
369 function INDEX (string haystack, string needles, needle_len_d)
371 if (needle_len_d <= INT_MIN || needle_len_d >= INT_MAX
372 || (int) needle_len_d != needle_len_d
373 || needles.length == 0)
377 int needle_len = needle_len_d;
378 if (needle_len < 0 || needle_len > needles.length
379 || needles.length % needle_len != 0)
383 int limit = haystack.length - needle_len + 1;
385 for (i = 1; i <= limit; i++)
386 for (j = 0; j < needles.length; j += needle_len)
387 if (!memcmp (&haystack.string[i - 1], &needles.string[j],
396 function RINDEX (string haystack, string needle)
398 if (needle.length == 0)
402 int limit = haystack.length - needle.length + 1;
404 for (i = limit; i >= 1; i--)
405 if (!memcmp (&haystack.string[i - 1], needle.string, needle.length))
411 function RINDEX (string haystack, string needles, needle_len_d)
413 if (needle_len_d <= INT_MIN || needle_len_d >= INT_MAX
414 || (int) needle_len_d != needle_len_d
415 || needles.length == 0)
419 int needle_len = needle_len_d;
420 if (needle_len < 0 || needle_len > needles.length
421 || needles.length % needle_len != 0)
425 int limit = haystack.length - needle_len + 1;
427 for (i = limit; i >= 1; i--)
428 for (j = 0; j < needles.length; j += needle_len)
429 if (!memcmp (&haystack.string[i - 1],
430 &needles.string[j], needle_len))
437 function LENGTH (string s)
442 string function LOWER (string s)
446 for (i = 0; i < s.length; i++)
447 s.string[i] = tolower ((unsigned char) s.string[i]);
451 function MBLEN.BYTE (string s, idx)
453 if (idx < 0 || idx >= s.length || (int) idx != idx)
459 string function UPCASE (string s)
463 for (i = 0; i < s.length; i++)
464 s.string[i] = toupper ((unsigned char) s.string[i]);
468 absorb_miss string function LPAD (string s, n)
471 if (n < 0 || n > MAX_STRING || (int) n != n)
473 else if (s.length >= n)
477 struct substring t = alloc_string (e, n);
478 memset (t.string, ' ', n - s.length);
479 memcpy (&t.string[(int) n - s.length], s.string, s.length);
484 absorb_miss string function LPAD (string s, n, string c)
487 if (n < 0 || n > MAX_STRING || (int) n != n || c.length != 1)
489 else if (s.length >= n)
493 struct substring t = alloc_string (e, n);
494 memset (t.string, c.string[0], n - s.length);
495 memcpy (&t.string[(int) n - s.length], s.string, s.length);
500 absorb_miss string function RPAD (string s, n)
503 if (n < 0 || n > MAX_STRING || (int) n != n)
505 else if (s.length >= n)
509 struct substring t = alloc_string (e, n);
510 memcpy (t.string, s.string, s.length);
511 memset (&t.string[s.length], ' ', n - s.length);
516 absorb_miss string function RPAD (string s, n, string c)
519 if (n < 0 || n > MAX_STRING || (int) n != n || c.length != 1)
521 else if (s.length >= n)
525 struct substring t = alloc_string (e, n);
526 memcpy (t.string, s.string, s.length);
527 memset (&t.string[s.length], c.string[0], n - s.length);
532 string function LTRIM (string s)
534 while (s.length > 0 && s.string[0] == ' ')
542 string function LTRIM (string s, string c)
546 while (s.length > 0 && s.string[0] == c.string[0])
557 string function RTRIM (string s)
559 while (s.length > 0 && s.string[s.length - 1] == ' ')
564 string function RTRIM (string s, string c)
568 while (s.length > 0 && s.string[s.length - 1] == c.string[0])
576 function NUMBER (string s, ni_format f)
579 data_in (ss_head (s, f->w), f->type, f->d, 0, &out, 0);
583 absorb_miss string function STRING (x, no_format f)
587 struct substring dst;
590 dst = alloc_string (e, f->w);
591 assert (!fmt_is_string (f->type));
592 data_out (&v, f, dst.string);
596 absorb_miss string function SUBSTR (string s, ofs)
599 if (ofs >= 1 && ofs <= s.length && (int) ofs == ofs)
600 return copy_string (e, &s.string[(int) ofs - 1], s.length - ofs + 1);
605 absorb_miss string function SUBSTR (string s, ofs, cnt)
608 if (ofs >= 1 && ofs <= s.length && (int) ofs == ofs
609 && cnt >= 1 && cnt <= INT_MAX && (int) cnt == cnt)
611 int cnt_max = s.length - (int) ofs + 1;
612 return copy_string (e, &s.string[(int) ofs - 1],
613 cnt <= cnt_max ? cnt : cnt_max);
619 absorb_miss no_opt no_abbrev string function VALUELABEL (var v)
623 const char *label = var_lookup_value_label (v, case_data (c, v));
625 return copy_string (e, label, strlen (label));
631 operator SQUARE (x) = x * x;
632 boolean operator NUM_TO_BOOLEAN (x)
634 if (x == 0. || x == 1. || x == SYSMIS)
638 msg (SE, _("A number being treated as a Boolean in an "
639 "expression was found to have a value other than "
640 "0 (false), 1 (true), or the system-missing value. "
641 "The result was forced to 0."));
646 operator BOOLEAN_TO_NUM (boolean x) = x;
648 // Beta distribution.
649 function PDF.BETA (x >= 0 && x <= 1, a > 0, b > 0)
650 = gsl_ran_beta_pdf (x, a, b);
651 function CDF.BETA (x >= 0 && x <= 1, a > 0, b > 0) = gsl_cdf_beta_P (x, a, b);
652 function IDF.BETA (P >= 0 && P <= 1, a > 0, b > 0)
653 = gslextras_cdf_beta_Pinv (P, a, b);
654 no_opt function RV.BETA (a > 0, b > 0) = gsl_ran_beta (get_rng (), a, b);
655 function NCDF.BETA (x >= 0, a > 0, b > 0, lambda > 0)
656 = ncdf_beta (x, a, b, lambda);
657 function NPDF.BETA (x >= 0, a > 0, b > 0, lambda > 0)
658 = npdf_beta (x, a, b, lambda);
660 // Bivariate normal distribution.
661 function CDF.BVNOR (x0, x1, r >= -1 && r <= 1) = cdf_bvnor (x0, x1, r);
662 function PDF.BVNOR (x0, x1, r >= -1 && r <= 1)
663 = gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
665 // Cauchy distribution.
666 function CDF.CAUCHY (x, a, b > 0) = gsl_cdf_cauchy_P ((x - a) / b, 1);
667 function IDF.CAUCHY (P > 0 && P < 1, a, b > 0)
668 = a + b * gsl_cdf_cauchy_Pinv (P, 1);
669 function PDF.CAUCHY (x, a, b > 0) = gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
670 no_opt function RV.CAUCHY (a, b > 0) = a + b * gsl_ran_cauchy (get_rng (), 1);
672 // Chi-square distribution.
673 function CDF.CHISQ (x >= 0, df > 0) = gsl_cdf_chisq_P (x, df);
674 function IDF.CHISQ (P >= 0 && P < 1, df > 0) = gsl_cdf_chisq_Pinv (P, df);
675 function PDF.CHISQ (x >= 0, df > 0) = gsl_ran_chisq_pdf (x, df);
676 no_opt function RV.CHISQ (df > 0) = gsl_ran_chisq (get_rng (), df);
677 function NCDF.CHISQ (x >= 0, df > 0, c) = unimplemented;
678 function NPDF.CHISQ (x >= 0, df > 0, c) = unimplemented;
679 function SIG.CHISQ (x >= 0, df > 0) = gsl_cdf_chisq_Q (x, df);
681 // Exponential distribution.
682 function CDF.EXP (x >= 0, a > 0) = gsl_cdf_exponential_P (x, 1. / a);
683 function IDF.EXP (P >= 0 && P < 1, a > 0)
684 = gsl_cdf_exponential_Pinv (P, 1. / a);
685 function PDF.EXP (x >= 0, a > 0) = gsl_ran_exponential_pdf (x, 1. / a);
686 no_opt function RV.EXP (a > 0) = gsl_ran_exponential (get_rng (), 1. / a);
688 // Exponential power distribution.
689 extension function PDF.XPOWER (x, a > 0, b >= 0)
690 = gsl_ran_exppow_pdf (x, a, b);
691 no_opt extension function RV.XPOWER (a > 0, b >= 0)
692 = gsl_ran_exppow (get_rng (), a, b);
695 function CDF.F (x >= 0, df1 > 0, df2 > 0) = gsl_cdf_fdist_P (x, df1, df2);
696 function IDF.F (P >= 0 && P < 1, df1 > 0, df2 > 0) = idf_fdist (P, df1, df2);
697 function PDF.F (x >= 0, df1 > 0, df2 > 0) = gsl_ran_fdist_pdf (x, df1, df2);
698 no_opt function RV.F (df1 > 0, df2 > 0) = gsl_ran_fdist (get_rng (), df1, df2);
699 function NCDF.F (x >= 0, df1 > 0, df2 > 0, lambda >= 0) = unimplemented;
700 function NPDF.F (x >= 0, df1 > 0, df2 > 0, lmabda >= 0) = unimplemented;
701 function SIG.F (x >= 0, df1 > 0, df2 > 0) = gsl_cdf_fdist_Q (x, df1, df2);
703 // Gamma distribution.
704 function CDF.GAMMA (x >= 0, a > 0, b > 0) = gsl_cdf_gamma_P (x, a, 1. / b);
705 function IDF.GAMMA (P >= 0 && P <= 1, a > 0, b > 0)
706 = gsl_cdf_gamma_Pinv (P, a, 1. / b);
707 function PDF.GAMMA (x >= 0, a > 0, b > 0) = gsl_ran_gamma_pdf (x, a, 1. / b);
708 no_opt function RV.GAMMA (a > 0, b > 0)
709 = gsl_ran_gamma (get_rng (), a, 1. / b);
711 // Half-normal distribution.
712 function CDF.HALFNRM (x, a, b > 0) = unimplemented;
713 function IDF.HALFNRM (P > 0 && P < 1, a, b > 0) = unimplemented;
714 function PDF.HALFNRM (x, a, b > 0) = unimplemented;
715 no_opt function RV.HALFNRM (a, b > 0) = unimplemented;
717 // Inverse Gaussian distribution.
718 function CDF.IGAUSS (x > 0, a > 0, b > 0) = unimplemented;
719 function IDF.IGAUSS (P >= 0 && P < 1, a > 0, b > 0) = unimplemented;
720 function PDF.IGAUSS (x > 0, a > 0, b > 0) = unimplemented;
721 no_opt function RV.IGAUSS (a > 0, b > 0) = unimplemented;
723 // Landau distribution.
724 extension function PDF.LANDAU (x) = gsl_ran_landau_pdf (x);
725 no_opt extension function RV.LANDAU () = gsl_ran_landau (get_rng ());
727 // Laplace distribution.
728 function CDF.LAPLACE (x, a, b > 0) = gsl_cdf_laplace_P ((x - a) / b, 1);
729 function IDF.LAPLACE (P > 0 && P < 1, a, b > 0)
730 = a + b * gsl_cdf_laplace_Pinv (P, 1);
731 function PDF.LAPLACE (x, a, b > 0) = gsl_ran_laplace_pdf ((x - a) / b, 1) / b;
732 no_opt function RV.LAPLACE (a, b > 0)
733 = a + b * gsl_ran_laplace (get_rng (), 1);
735 // Levy alpha-stable distribution.
736 no_opt extension function RV.LEVY (c, alpha > 0 && alpha <= 2)
737 = gsl_ran_levy (get_rng (), c, alpha);
739 // Levy skew alpha-stable distribution.
740 no_opt extension function RV.LVSKEW (c, alpha > 0 && alpha <= 2,
741 beta >= -1 && beta <= 1)
742 = gsl_ran_levy_skew (get_rng (), c, alpha, beta);
744 // Logistic distribution.
745 function CDF.LOGISTIC (x, a, b > 0) = gsl_cdf_logistic_P ((x - a) / b, 1);
746 function IDF.LOGISTIC (P > 0 && P < 1, a, b > 0)
747 = a + b * gsl_cdf_logistic_Pinv (P, 1);
748 function PDF.LOGISTIC (x, a, b > 0)
749 = gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
750 no_opt function RV.LOGISTIC (a, b > 0)
751 = a + b * gsl_ran_logistic (get_rng (), 1);
753 // Lognormal distribution.
754 function CDF.LNORMAL (x >= 0, m > 0, s > 0)
755 = gsl_cdf_lognormal_P (x, log (m), s);
756 function IDF.LNORMAL (P >= 0 && P < 1, m > 0, s > 0)
757 = gsl_cdf_lognormal_Pinv (P, log (m), s);
758 function PDF.LNORMAL (x >= 0, m > 0, s > 0)
759 = gsl_ran_lognormal_pdf (x, log (m), s);
760 no_opt function RV.LNORMAL (m > 0, s > 0)
761 = gsl_ran_lognormal (get_rng (), log (m), s);
763 // Normal distribution.
764 function CDF.NORMAL (x, u, s > 0) = gsl_cdf_gaussian_P (x - u, s);
765 function IDF.NORMAL (P > 0 && P < 1, u, s > 0)
766 = u + gsl_cdf_gaussian_Pinv (P, s);
767 function PDF.NORMAL (x, u, s > 0) = gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
768 no_opt function RV.NORMAL (u, s > 0) = u + gsl_ran_gaussian (get_rng (), s);
769 function CDFNORM (x) = gsl_cdf_ugaussian_P (x);
770 function PROBIT (P > 0 && P < 1) = gsl_cdf_ugaussian_Pinv (P);
771 no_opt function NORMAL (s > 0) = gsl_ran_gaussian (get_rng (), s);
773 // Normal tail distribution.
774 function PDF.NTAIL (x, a > 0, sigma > 0)
775 = gsl_ran_gaussian_tail_pdf (x, a, sigma);
776 no_opt function RV.NTAIL (a > 0, sigma > 0)
777 = gsl_ran_gaussian_tail (get_rng (), a, sigma);
779 // Pareto distribution.
780 function CDF.PARETO (x >= a, a > 0, b > 0) = gsl_cdf_pareto_P (x, b, a);
781 function IDF.PARETO (P >= 0 && P < 1, a > 0, b > 0)
782 = gsl_cdf_pareto_Pinv (P, b, a);
783 function PDF.PARETO (x >= a, a > 0, b > 0) = gsl_ran_pareto_pdf (x, b, a);
784 no_opt function RV.PARETO (a > 0, b > 0) = gsl_ran_pareto (get_rng (), b, a);
786 // Rayleigh distribution.
787 extension function CDF.RAYLEIGH (x, sigma > 0) = gsl_cdf_rayleigh_P (x, sigma);
788 extension function IDF.RAYLEIGH (P >= 0 && P <= 1, sigma > 0)
789 = gsl_cdf_rayleigh_Pinv (P, sigma);
790 extension function PDF.RAYLEIGH (x, sigma > 0)
791 = gsl_ran_rayleigh_pdf (x, sigma);
792 no_opt extension function RV.RAYLEIGH (sigma > 0)
793 = gsl_ran_rayleigh (get_rng (), sigma);
795 // Rayleigh tail distribution.
796 extension function PDF.RTAIL (x, a, sigma)
797 = gsl_ran_rayleigh_tail_pdf (x, a, sigma);
798 no_opt extension function RV.RTAIL (a, sigma)
799 = gsl_ran_rayleigh_tail (get_rng (), a, sigma);
801 // Studentized maximum modulus distribution.
802 function CDF.SMOD (x > 0, a >= 1, b >= 1) = unimplemented;
803 function IDF.SMOD (P >= 0 && P < 1, a >= 1, b >= 1) = unimplemented;
805 // Studentized range distribution.
806 function CDF.SRANGE (x > 0, a >= 1, b >= 1) = unimplemented;
807 function IDF.SRANGE (P >= 0 && P < 1, a >= 1, b >= 1) = unimplemented;
809 // Student t distribution.
810 function CDF.T (x, df > 0) = gsl_cdf_tdist_P (x, df);
811 function IDF.T (P > 0 && P < 1, df > 0) = gsl_cdf_tdist_Pinv (P, df);
812 function PDF.T (x, df > 0) = gsl_ran_tdist_pdf (x, df);
813 no_opt function RV.T (df > 0) = gsl_ran_tdist (get_rng (), df);
814 function NCDF.T (x, df > 0, nc) = unimplemented;
815 function NPDF.T (x, df > 0, nc) = unimplemented;
817 // Type-1 Gumbel distribution.
818 extension function CDF.T1G (x, a, b) = gsl_cdf_gumbel1_P (x, a, b);
819 extension function IDF.T1G (P >= 0 && P <= 1, a, b)
820 = gsl_cdf_gumbel1_P (P, a, b);
821 extension function PDF.T1G (x, a, b) = gsl_ran_gumbel1_pdf (x, a, b);
822 no_opt extension function RV.T1G (a, b) = gsl_ran_gumbel1 (get_rng (), a, b);
824 // Type-2 Gumbel distribution.
825 extension function CDF.T2G (x, a, b) = gsl_cdf_gumbel2_P (x, a, b);
826 extension function IDF.T2G (P >= 0 && P <= 1, a, b)
827 = gsl_cdf_gumbel2_P (P, a, b);
828 extension function PDF.T2G (x, a, b) = gsl_ran_gumbel2_pdf (x, a, b);
829 no_opt extension function RV.T2G (a, b) = gsl_ran_gumbel2 (get_rng (), a, b);
831 // Uniform distribution.
832 function CDF.UNIFORM (x <= b, a <= x, b) = gsl_cdf_flat_P (x, a, b);
833 function IDF.UNIFORM (P >= 0 && P <= 1, a <= b, b)
834 = gsl_cdf_flat_Pinv (P, a, b);
835 function PDF.UNIFORM (x <= b, a <= x, b) = gsl_ran_flat_pdf (x, a, b);
836 no_opt function RV.UNIFORM (a <= b, b) = gsl_ran_flat (get_rng (), a, b);
837 no_opt function UNIFORM (b >= 0) = gsl_ran_flat (get_rng (), 0, b);
839 // Weibull distribution.
840 function CDF.WEIBULL (x >= 0, a > 0, b > 0) = gsl_cdf_weibull_P (x, a, b);
841 function IDF.WEIBULL (P >= 0 && P < 1, a > 0, b > 0)
842 = gsl_cdf_weibull_Pinv (P, a, b);
843 function PDF.WEIBULL (x >= 0, a > 0, b > 0) = gsl_ran_weibull_pdf (x, a, b);
844 no_opt function RV.WEIBULL (a > 0, b > 0) = gsl_ran_weibull (get_rng (), a, b);
846 // Bernoulli distribution.
847 function CDF.BERNOULLI (k == 0 || k == 1, p >= 0 && p <= 1)
849 function PDF.BERNOULLI (k == 0 || k == 1, p >= 0 && p <= 1)
850 = gsl_ran_bernoulli_pdf (k, p);
851 no_opt function RV.BERNOULLI (p >= 0 && p <= 1)
852 = gsl_ran_bernoulli (get_rng (), p);
854 // Binomial distribution.
855 function CDF.BINOM (k, n > 0 && n == floor (n), p >= 0 && p <= 1)
856 = gslextras_cdf_binomial_P (k, p, n);
857 function PDF.BINOM (k >= 0 && k == floor (k) && k <= n,
858 n > 0 && n == floor (n),
860 = gsl_ran_binomial_pdf (k, p, n);
861 no_opt function RV.BINOM (p > 0 && p == floor (p), n >= 0 && n <= 1)
862 = gsl_ran_binomial (get_rng (), p, n);
864 // Geometric distribution.
865 function CDF.GEOM (k >= 1 && k == floor (k), p >= 0 && p <= 1)
866 = gslextras_cdf_geometric_P (k, p);
867 function PDF.GEOM (k >= 1 && k == floor (k),
869 = gsl_ran_geometric_pdf (k, p);
870 no_opt function RV.GEOM (p >= 0 && p <= 1) = gsl_ran_geometric (get_rng (), p);
872 // Hypergeometric distribution.
873 function CDF.HYPER (k >= 0 && k == floor (k) && k <= c,
874 a > 0 && a == floor (a),
875 b > 0 && b == floor (b) && b <= a,
876 c > 0 && c == floor (c) && c <= a)
877 = gslextras_cdf_hypergeometric_P (k, c, a - c, b);
878 function PDF.HYPER (k >= 0 && k == floor (k) && k <= c,
879 a > 0 && a == floor (a),
880 b > 0 && b == floor (b) && b <= a,
881 c > 0 && c == floor (c) && c <= a)
882 = gsl_ran_hypergeometric_pdf (k, c, a - c, b);
883 no_opt function RV.HYPER (a > 0 && a == floor (a),
884 b > 0 && b == floor (b) && b <= a,
885 c > 0 && c == floor (c) && c <= a)
886 = gsl_ran_hypergeometric (get_rng (), c, a - c, b);
888 // Logarithmic distribution.
889 extension function PDF.LOG (k >= 1, p > 0 && p <= 1)
890 = gsl_ran_logarithmic_pdf (k, p);
891 no_opt extension function RV.LOG (p > 0 && p <= 1)
892 = gsl_ran_logarithmic (get_rng (), p);
894 // Negative binomial distribution.
895 function CDF.NEGBIN (k >= 1, n == floor (n), p > 0 && p <= 1)
896 = gslextras_cdf_negative_binomial_P (k, p, n);
897 function PDF.NEGBIN (k >= 1, n == floor (n), p > 0 && p <= 1)
898 = gsl_ran_negative_binomial_pdf (k, p, n);
899 no_opt function RV.NEGBIN (n == floor (n), p > 0 && p <= 1)
900 = gsl_ran_negative_binomial (get_rng (), p, n);
902 // Poisson distribution.
903 function CDF.POISSON (k >= 0 && k == floor (k), mu > 0)
904 = gslextras_cdf_poisson_P (k, mu);
905 function PDF.POISSON (k >= 0 && k == floor (k), mu > 0)
906 = gsl_ran_poisson_pdf (k, mu);
907 no_opt function RV.POISSON (mu > 0) = gsl_ran_poisson (get_rng (), mu);
910 absorb_miss boolean function MISSING (x) = x == SYSMIS || !finite (x);
911 absorb_miss boolean function SYSMIS (x) = x == SYSMIS || !finite (x);
912 no_opt boolean function SYSMIS (num_var v)
915 return case_num (c, v) == SYSMIS;
917 no_opt boolean function VALUE (num_var v)
920 return case_num (c, v);
923 no_opt operator VEC_ELEM_NUM (idx)
927 if (idx >= 1 && idx <= vector_get_var_cnt (v))
929 const struct variable *var = vector_get_var (v, (size_t) idx - 1);
930 double value = case_num (c, var);
931 return !var_is_num_user_missing (var, value) ? value : SYSMIS;
936 msg (SE, _("SYSMIS is not a valid index value for vector "
937 "%s. The result will be set to SYSMIS."),
938 vector_get_name (v));
940 msg (SE, _("%g is not a valid index value for vector %s. "
941 "The result will be set to SYSMIS."),
942 idx, vector_get_name (v));
947 absorb_miss no_opt string operator VEC_ELEM_STR (idx)
952 if (idx >= 1 && idx <= vector_get_var_cnt (v))
954 struct variable *var = vector_get_var (v, (size_t) idx - 1);
955 return copy_string (e, case_str (c, var), var_get_width (var));
960 msg (SE, _("SYSMIS is not a valid index value for vector "
961 "%s. The result will be set to the empty string."),
962 vector_get_name (v));
964 msg (SE, _("%g is not a valid index value for vector %s. "
965 "The result will be set to the empty string."),
966 idx, vector_get_name (v));
973 no_opt operator NUM_VAR ()
977 double d = case_num (c, v);
978 return !var_is_num_user_missing (v, d) ? d : SYSMIS;
981 no_opt string operator STR_VAR ()
986 struct substring s = alloc_string (e, var_get_width (v));
987 memcpy (s.string, case_str (c, v), var_get_width (v));
991 no_opt perm_only function LAG (num_var v, pos_int n_before)
994 struct ccase *c = lagged_case (ds, n_before);
997 double x = case_num (c, v);
998 return !var_is_num_user_missing (v, x) ? x : SYSMIS;
1004 no_opt perm_only function LAG (num_var v)
1007 struct ccase *c = lagged_case (ds, 1);
1010 double x = case_num (c, v);
1011 return !var_is_num_user_missing (v, x) ? x : SYSMIS;
1017 no_opt perm_only string function LAG (str_var v, pos_int n_before)
1021 struct ccase *c = lagged_case (ds, n_before);
1023 return copy_string (e, case_str (c, v), var_get_width (v));
1025 return empty_string;
1028 no_opt perm_only string function LAG (str_var v)
1032 struct ccase *c = lagged_case (ds, 1);
1034 return copy_string (e, case_str (c, v), var_get_width (v));
1036 return empty_string;
1039 no_opt operator NUM_SYS ()
1043 return case_num (c, v) == SYSMIS;
1046 no_opt operator NUM_VAL ()
1050 return case_num (c, v);
1053 no_opt operator CASENUM ()