src/language/xforms/sample.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <gsl/gsl_rng.h>
  20 #include <limits.h>
  21 #include <stdio.h>
  22 #include <math.h>
  23
  24 #include <data/procedure.h>
  25 #include <data/variable.h>
  26 #include <language/command.h>
  27 #include <language/lexer/lexer.h>
  28 #include <libpspp/alloc.h>
  29 #include <libpspp/compiler.h>
  30 #include <libpspp/message.h>
  31 #include <libpspp/str.h>
  32 #include <math/random.h>
  33
  34 #include "gettext.h"
  35 #define _(msgid) gettext (msgid)
  36
  37 /* The two different types of samples. */
  38 enum
  39   {
  40     TYPE_A_FROM_B,              /* 5 FROM 10 */
  41     TYPE_FRACTION               /* 0.5 */
  42   };
  43
  44 /* SAMPLE transformation. */
  45 struct sample_trns
  46   {
  47     int type;                   /* One of TYPE_*. */
  48     int n, N;                   /* TYPE_A_FROM_B: n from N. */
  49     int m, t;                   /* TYPE_A_FROM_B: # picked so far; # so far. */
  50     unsigned frac;              /* TYPE_FRACTION: a fraction of UINT_MAX. */
  51   };
  52
  53 static trns_proc_func sample_trns_proc;
  54 static trns_free_func sample_trns_free;
  55
  56 int
  57 cmd_sample (struct lexer *lexer, struct dataset *ds)
  58 {
  59   struct sample_trns *trns;
  60
  61   int type;
  62   int a, b;
  63   unsigned frac;
  64
  65   if (!lex_force_num (lexer))
  66     return CMD_FAILURE;
  67   if (!lex_is_integer (lexer))
  68     {
  69       unsigned long min = gsl_rng_min (get_rng ());
  70       unsigned long max = gsl_rng_max (get_rng ());
  71
  72       type = TYPE_FRACTION;
  73       if (lex_tokval (lexer) <= 0 || lex_tokval (lexer) >= 1)
  74         {
  75           msg (SE, _("The sampling factor must be between 0 and 1 "
  76                      "exclusive."));
  77           return CMD_FAILURE;
  78         }
  79
  80       frac = lex_tokval (lexer) * (max - min) + min;
  81       a = b = 0;
  82     }
  83   else
  84     {
  85       type = TYPE_A_FROM_B;
  86       a = lex_integer (lexer);
  87       lex_get (lexer);
  88       if (!lex_force_match_id (lexer, "FROM"))
  89         return CMD_FAILURE;
  90       if (!lex_force_int (lexer))
  91         return CMD_FAILURE;
  92       b = lex_integer (lexer);
  93       if (a >= b)
  94         {
  95           msg (SE, _("Cannot sample %d observations from a population of "
  96                      "%d."),
  97                a, b);
  98           return CMD_FAILURE;
  99         }
 100
 101       frac = 0;
 102     }
 103   lex_get (lexer);
 104
 105   trns = xmalloc (sizeof *trns);
 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 (ds, sample_trns_proc, sample_trns_free, trns);
 112
 113   return lex_end_of_command (lexer);
 114 }
 115
 116 /* Executes a SAMPLE transformation. */
 117 static int
 118 sample_trns_proc (void *t_, struct ccase *c UNUSED,
 119                   casenumber case_num UNUSED)
 120 {
 121   struct sample_trns *t = t_;
 122   double U;
 123
 124   if (t->type == TYPE_FRACTION)
 125     {
 126       if (gsl_rng_get (get_rng ()) <= t->frac)
 127         return TRNS_CONTINUE;
 128       else
 129         return TRNS_DROP_CASE;
 130     }
 131
 132   if (t->m >= t->n)
 133     return TRNS_DROP_CASE;
 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 TRNS_DROP_CASE;
 140     }
 141   else
 142     {
 143       t->m++;
 144       t->t++;
 145       return TRNS_CONTINUE;
 146     }
 147 }
 148
 149 static bool
 150 sample_trns_free (void *t_)
 151 {
 152   struct sample_trns *t = t_;
 153   free (t);
 154   return true;
 155 }