b721f2db9c1228464ed2b2c30486a63cfbf8cc08
[pspp] / src / sample.c
1 /* PSPP - computes sample statistics.
2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3    Written by Ben Pfaff <blp@gnu.org>.
4
5    This program is free software; you can redistribute it and/or
6    modify it under the terms of the GNU General Public License as
7    published by the Free Software Foundation; either version 2 of the
8    License, or (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful, but
11    WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13    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, write to the Free Software
17    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
18    02111-1307, USA. */
19
20 #include <config.h>
21 #include <limits.h>
22 #include <stdio.h>
23 #include <math.h>
24 #include "alloc.h"
25 #include "command.h"
26 #include "error.h"
27 #include "lexer.h"
28 #include "random.h"
29 #include "str.h"
30 #include "var.h"
31
32 #include "debug-print.h"
33
34 /* The two different types of samples. */
35 enum
36   {
37     TYPE_A_FROM_B,              /* 5 FROM 10 */
38     TYPE_FRACTION               /* 0.5 */
39   };
40
41 /* SAMPLE transformation. */
42 struct sample_trns
43   {
44     struct trns_header h;
45     int type;                   /* One of TYPE_*. */
46     int n, N;                   /* TYPE_A_FROM_B: n from N. */
47     int m, t;                   /* TYPE_A_FROM_B: # picked so far; # so far. */
48     unsigned frac;              /* TYPE_FRACTION: a fraction of UINT_MAX. */
49   };
50
51 int sample_trns_proc (struct trns_header *, struct ccase *);
52
53 int
54 cmd_sample (void)
55 {
56   struct sample_trns *trns;
57
58   int type;
59   int a, b;
60   unsigned frac;
61
62   lex_match_id ("SAMPLE");
63
64   if (!lex_force_num ())
65     return CMD_FAILURE;
66   if (!lex_integer_p ())
67     {
68       type = TYPE_FRACTION;
69       if (tokval <= 0 || tokval >= 1)
70         {
71           msg (SE, _("The sampling factor must be between 0 and 1 "
72                      "exclusive."));
73           return CMD_FAILURE;
74         }
75           
76       frac = tokval * UINT_MAX;
77       a = b = 0;
78     }
79   else
80     {
81       type = TYPE_A_FROM_B;
82       a = lex_integer ();
83       lex_get ();
84       if (!lex_force_match_id ("FROM"))
85         return CMD_FAILURE;
86       if (!lex_force_int ())
87         return CMD_FAILURE;
88       b = lex_integer ();
89       if (a >= b)
90         {
91           msg (SE, _("Cannot sample %d observations from a population of "
92                      "%d."),
93                a, b);
94           return CMD_FAILURE;
95         }
96       
97       frac = 0;
98     }
99   lex_get ();
100
101 #if DEBUGGING
102   if (type == TYPE_FRACTION)
103     printf ("SAMPLE %g.\n", frac / (double) UINT_MAX);
104   else
105     printf ("SAMPLE %d FROM %d.\n", a, b);
106 #endif
107
108   trns = xmalloc (sizeof *trns);
109   trns->h.proc = sample_trns_proc;
110   trns->h.free = NULL;
111   trns->type = type;
112   trns->n = a;
113   trns->N = b;
114   trns->m = trns->t = 0;
115   trns->frac = frac;
116   add_transformation ((struct trns_header *) trns);
117
118   return lex_end_of_command ();
119 }
120
121 int
122 sample_trns_proc (struct trns_header * trns, struct ccase *c unused)
123 {
124   struct sample_trns *t = (struct sample_trns *) trns;
125   double U;
126
127   if (t->type == TYPE_FRACTION) 
128     {
129       if (rng_get_unsigned (pspp_rng ()) <= t->frac)
130         return -1;
131       else
132         return -2;
133     }
134
135   if (t->m >= t->n)
136     return -2;
137
138   U = rng_get_double (pspp_rng ());
139   if ((t->N - t->t) * U >= t->n - t->m)
140     {
141       t->t++;
142       return -2;
143     }
144   else
145     {
146       t->m++;
147       t->t++;
148       return -1;
149     }
150 }