84f2e8992f65340a5062be3afc7f7afe682dffcb
[pspp-builds.git] / 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 <stdio.h>
22 #include <math.h>
23 #include "alloc.h"
24 #include "command.h"
25 #include "error.h"
26 #include "lexer.h"
27 #include "random.h"
28 #include "str.h"
29 #include "var.h"
30
31 #undef DEBUGGING
32 /*#define DEBUGGING 1 */
33 #include "debug-print.h"
34
35 /* The two different types of samples. */
36 enum
37   {
38     TYPE_A_FROM_B,              /* 5 FROM 10 */
39     TYPE_FRACTION               /* 0.5 */
40   };
41
42 /* SAMPLE transformation. */
43 struct sample_trns
44   {
45     struct trns_header h;
46     int type;                   /* One of TYPE_*. */
47     int n, N;                   /* TYPE_A_FROM_B: n from N. */
48     int m, t;                   /* TYPE_A_FROM_B: # selected so far; # so far. */
49     int frac;                   /* TYPE_FRACTION: a fraction out of 65536. */
50   };
51
52 int sample_trns_proc (struct trns_header *, struct ccase *);
53
54 int
55 cmd_sample (void)
56 {
57   struct sample_trns *trns;
58
59   int type;
60   int a, b;
61   int frac;
62
63   lex_match_id ("SAMPLE");
64
65   if (!lex_force_num ())
66     return CMD_FAILURE;
67   if (!lex_integer_p ())
68     {
69       type = TYPE_FRACTION;
70       if (tokval <= 0 || tokval >= 1)
71         {
72           msg (SE, _("The sampling factor must be between 0 and 1 "
73                      "exclusive."));
74           return CMD_FAILURE;
75         }
76           
77       frac = tokval * 65536;
78       a = b = 0;
79     }
80   else
81     {
82       type = TYPE_A_FROM_B;
83       a = lex_integer ();
84       lex_get ();
85       if (!lex_force_match_id ("FROM"))
86         return CMD_FAILURE;
87       if (!lex_force_int ())
88         return CMD_FAILURE;
89       b = lex_integer ();
90       if (a >= b)
91         {
92           msg (SE, _("Cannot sample %d observations from a population of "
93                      "%d."),
94                a, b);
95           return CMD_FAILURE;
96         }
97       
98       frac = 0;
99     }
100   lex_get ();
101
102 #if DEBUGGING
103   if (type == TYPE_FRACTION)
104     printf ("SAMPLE %g.\n", frac / 65536.);
105   else
106     printf ("SAMPLE %d FROM %d.\n", a, b);
107 #endif
108
109   trns = xmalloc (sizeof *trns);
110   trns->h.proc = sample_trns_proc;
111   trns->h.free = NULL;
112   trns->type = type;
113   trns->n = a;
114   trns->N = b;
115   trns->m = trns->t = 0;
116   trns->frac = frac;
117   add_transformation ((struct trns_header *) trns);
118
119   return lex_end_of_command ();
120 }
121
122 int
123 sample_trns_proc (struct trns_header * trns, struct ccase *c unused)
124 {
125   struct sample_trns *t = (struct sample_trns *) trns;
126   double U;
127
128   if (t->type == TYPE_FRACTION)
129     return (rand_simple (0x10000) <= t->frac) - 2;
130
131   if (t->m >= t->n)
132     return -2;
133
134   U = rand_uniform (1);
135   if ((t->N - t->t) * U >= t->n - t->m)
136     {
137       t->t++;
138       return -2;
139     }
140   else
141     {
142       t->m++;
143       t->t++;
144       return -1;
145     }
146 }