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