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., 51 Franklin Street, Fifth Floor, Boston, MA
  18    02110-1301, USA. */
  19
  20 #include <config.h>
  21 #include <gsl/gsl_rng.h>
  22 #include <limits.h>
  23 #include <stdio.h>
  24 #include <math.h>
  25 #include "alloc.h"
  26 #include "command.h"
  27 #include "error.h"
  28 #include "lexer.h"
  29 #include "settings.h"
  30 #include "str.h"
  31 #include "var.h"
  32
  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: # picked so far; # so far. */
  49     unsigned frac;              /* TYPE_FRACTION: a fraction of UINT_MAX. */
  50   };
  51
  52 static trns_proc_func sample_trns_proc;
  53
  54 int
  55 cmd_sample (void)
  56 {
  57   struct sample_trns *trns;
  58
  59   int type;
  60   int a, b;
  61   unsigned frac;
  62
  63   if (!lex_force_num ())
  64     return CMD_FAILURE;
  65   if (!lex_is_integer ())
  66     {
  67       unsigned long min = gsl_rng_min (get_rng ());
  68       unsigned long max = gsl_rng_max (get_rng ());
  69
  70       type = TYPE_FRACTION;
  71       if (tokval <= 0 || tokval >= 1)
  72         {
  73           msg (SE, _("The sampling factor must be between 0 and 1 "
  74                      "exclusive."));
  75           return CMD_FAILURE;
  76         }
  77
  78       frac = tokval * (max - min) + min;
  79       a = b = 0;
  80     }
  81   else
  82     {
  83       type = TYPE_A_FROM_B;
  84       a = lex_integer ();
  85       lex_get ();
  86       if (!lex_force_match_id ("FROM"))
  87         return CMD_FAILURE;
  88       if (!lex_force_int ())
  89         return CMD_FAILURE;
  90       b = lex_integer ();
  91       if (a >= b)
  92         {
  93           msg (SE, _("Cannot sample %d observations from a population of "
  94                      "%d."),
  95                a, b);
  96           return CMD_FAILURE;
  97         }
  98
  99       frac = 0;
 100     }
 101   lex_get ();
 102
 103   trns = xmalloc (sizeof *trns);
 104   trns->h.proc = sample_trns_proc;
 105   trns->h.free = NULL;
 106   trns->type = type;
 107   trns->n = a;
 108   trns->N = b;
 109   trns->m = trns->t = 0;
 110   trns->frac = frac;
 111   add_transformation ((struct trns_header *) trns);
 112
 113   return lex_end_of_command ();
 114 }
 115
 116 /* Executes a SAMPLE transformation. */
 117 static int
 118 sample_trns_proc (struct trns_header * trns, struct ccase *c UNUSED,
 119                   int case_num UNUSED)
 120 {
 121   struct sample_trns *t = (struct sample_trns *) trns;
 122   double U;
 123
 124   if (t->type == TYPE_FRACTION)
 125     {
 126       if (gsl_rng_get (get_rng ()) <= t->frac)
 127         return -1;
 128       else
 129         return -2;
 130     }
 131
 132   if (t->m >= t->n)
 133     return -2;
 134
 135   U = gsl_rng_uniform (get_rng ());
 136   if ((t->N - t->t) * U >= t->n - t->m)
 137     {
 138       t->t++;
 139       return -2;
 140     }
 141   else
 142     {
 143       t->m++;
 144       t->t++;
 145       return -1;
 146     }
 147 }