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 static trns_proc_func sample_trns_proc;
  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 /* Executes a SAMPLE transformation. */
 122 static int
 123 sample_trns_proc (struct trns_header * trns, struct ccase *c UNUSED,
 124                   int case_num UNUSED)
 125 {
 126   struct sample_trns *t = (struct sample_trns *) trns;
 127   double U;
 128
 129   if (t->type == TYPE_FRACTION)
 130     {
 131       if (rng_get_unsigned (pspp_rng ()) <= t->frac)
 132         return -1;
 133       else
 134         return -2;
 135     }
 136
 137   if (t->m >= t->n)
 138     return -2;
 139
 140   U = rng_get_double (pspp_rng ());
 141   if ((t->N - t->t) * U >= t->n - t->m)
 142     {
 143       t->t++;
 144       return -2;
 145     }
 146   else
 147     {
 148       t->m++;
 149       t->t++;
 150       return -1;
 151     }
 152 }