1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include <gsl/gsl_rng.h>
24 #include <data/procedure.h>
25 #include <data/variable.h>
26 #include <language/command.h>
27 #include <language/lexer/lexer.h>
28 #include <libpspp/alloc.h>
29 #include <libpspp/compiler.h>
30 #include <libpspp/message.h>
31 #include <libpspp/str.h>
32 #include <math/random.h>
35 #define _(msgid) gettext (msgid)
37 /* The two different types of samples. */
40 TYPE_A_FROM_B, /* 5 FROM 10 */
41 TYPE_FRACTION /* 0.5 */
44 /* SAMPLE transformation. */
47 int type; /* One of TYPE_*. */
48 int n, N; /* TYPE_A_FROM_B: n from N. */
49 int m, t; /* TYPE_A_FROM_B: # picked so far; # so far. */
50 unsigned frac; /* TYPE_FRACTION: a fraction of UINT_MAX. */
53 static trns_proc_func sample_trns_proc;
54 static trns_free_func sample_trns_free;
57 cmd_sample (struct lexer *lexer, struct dataset *ds)
59 struct sample_trns *trns;
65 if (!lex_force_num (lexer))
67 if (!lex_is_integer (lexer))
69 unsigned long min = gsl_rng_min (get_rng ());
70 unsigned long max = gsl_rng_max (get_rng ());
73 if (lex_tokval (lexer) <= 0 || lex_tokval (lexer) >= 1)
75 msg (SE, _("The sampling factor must be between 0 and 1 "
80 frac = lex_tokval (lexer) * (max - min) + min;
86 a = lex_integer (lexer);
88 if (!lex_force_match_id (lexer, "FROM"))
90 if (!lex_force_int (lexer))
92 b = lex_integer (lexer);
95 msg (SE, _("Cannot sample %d observations from a population of "
105 trns = xmalloc (sizeof *trns);
109 trns->m = trns->t = 0;
111 add_transformation (ds, sample_trns_proc, sample_trns_free, trns);
113 return lex_end_of_command (lexer);
116 /* Executes a SAMPLE transformation. */
118 sample_trns_proc (void *t_, struct ccase *c UNUSED,
119 casenumber case_num UNUSED)
121 struct sample_trns *t = t_;
124 if (t->type == TYPE_FRACTION)
126 if (gsl_rng_get (get_rng ()) <= t->frac)
127 return TRNS_CONTINUE;
129 return TRNS_DROP_CASE;
133 return TRNS_DROP_CASE;
135 U = gsl_rng_uniform (get_rng ());
136 if ((t->N - t->t) * U >= t->n - t->m)
139 return TRNS_DROP_CASE;
145 return TRNS_CONTINUE;
150 sample_trns_free (void *t_)
152 struct sample_trns *t = t_;