src/language/xforms/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
  22 #include <gsl/gsl_rng.h>
  23 #include <limits.h>
  24 #include <stdio.h>
  25 #include <math.h>
  26
  27 #include <data/procedure.h>
  28 #include <data/variable.h>
  29 #include <language/command.h>
  30 #include <language/lexer/lexer.h>
  31 #include <libpspp/alloc.h>
  32 #include <libpspp/compiler.h>
  33 #include <libpspp/message.h>
  34 #include <libpspp/str.h>
  35 #include <math/random.h>
  36
  37 #include "gettext.h"
  38 #define _(msgid) gettext (msgid)
  39
  40 /* The two different types of samples. */
  41 enum
  42   {
  43     TYPE_A_FROM_B,              /* 5 FROM 10 */
  44     TYPE_FRACTION               /* 0.5 */
  45   };
  46
  47 /* SAMPLE transformation. */
  48 struct sample_trns
  49   {
  50     int type;                   /* One of TYPE_*. */
  51     int n, N;                   /* TYPE_A_FROM_B: n from N. */
  52     int m, t;                   /* TYPE_A_FROM_B: # picked so far; # so far. */
  53     unsigned frac;              /* TYPE_FRACTION: a fraction of UINT_MAX. */
  54   };
  55
  56 static trns_proc_func sample_trns_proc;
  57 static trns_free_func sample_trns_free;
  58
  59 int
  60 cmd_sample (struct lexer *lexer, struct dataset *ds)
  61 {
  62   struct sample_trns *trns;
  63
  64   int type;
  65   int a, b;
  66   unsigned frac;
  67
  68   if (!lex_force_num (lexer))
  69     return CMD_FAILURE;
  70   if (!lex_is_integer (lexer))
  71     {
  72       unsigned long min = gsl_rng_min (get_rng ());
  73       unsigned long max = gsl_rng_max (get_rng ());
  74
  75       type = TYPE_FRACTION;
  76       if (lex_tokval (lexer) <= 0 || lex_tokval (lexer) >= 1)
  77         {
  78           msg (SE, _("The sampling factor must be between 0 and 1 "
  79                      "exclusive."));
  80           return CMD_FAILURE;
  81         }
  82
  83       frac = lex_tokval (lexer) * (max - min) + min;
  84       a = b = 0;
  85     }
  86   else
  87     {
  88       type = TYPE_A_FROM_B;
  89       a = lex_integer (lexer);
  90       lex_get (lexer);
  91       if (!lex_force_match_id (lexer, "FROM"))
  92         return CMD_FAILURE;
  93       if (!lex_force_int (lexer))
  94         return CMD_FAILURE;
  95       b = lex_integer (lexer);
  96       if (a >= b)
  97         {
  98           msg (SE, _("Cannot sample %d observations from a population of "
  99                      "%d."),
 100                a, b);
 101           return CMD_FAILURE;
 102         }
 103
 104       frac = 0;
 105     }
 106   lex_get (lexer);
 107
 108   trns = xmalloc (sizeof *trns);
 109   trns->type = type;
 110   trns->n = a;
 111   trns->N = b;
 112   trns->m = trns->t = 0;
 113   trns->frac = frac;
 114   add_transformation (ds, sample_trns_proc, sample_trns_free, trns);
 115
 116   return lex_end_of_command (lexer);
 117 }
 118
 119 /* Executes a SAMPLE transformation. */
 120 static int
 121 sample_trns_proc (void *t_, struct ccase *c UNUSED,
 122                   casenumber case_num UNUSED)
 123 {
 124   struct sample_trns *t = t_;
 125   double U;
 126
 127   if (t->type == TYPE_FRACTION)
 128     {
 129       if (gsl_rng_get (get_rng ()) <= t->frac)
 130         return TRNS_CONTINUE;
 131       else
 132         return TRNS_DROP_CASE;
 133     }
 134
 135   if (t->m >= t->n)
 136     return TRNS_DROP_CASE;
 137
 138   U = gsl_rng_uniform (get_rng ());
 139   if ((t->N - t->t) * U >= t->n - t->m)
 140     {
 141       t->t++;
 142       return TRNS_DROP_CASE;
 143     }
 144   else
 145     {
 146       t->m++;
 147       t->t++;
 148       return TRNS_CONTINUE;
 149     }
 150 }
 151
 152 static bool
 153 sample_trns_free (void *t_)
 154 {
 155   struct sample_trns *t = t_;
 156   free (t);
 157   return true;
 158 }