src/math/linreg/predict.c

   1 /*
   2    lib/linreg/predict.c
   3
   4    Copyright (C) 2005 Free Software Foundation, Inc. Written by Jason H. Stover.
   5
   6    This program is free software; you can redistribute it and/or modify it under
   7    the terms of the GNU General Public License as published by the Free
   8    Software Foundation; either version 2 of the License, or (at your option)
   9    any later version.
  10
  11    This program is distributed in the hope that it will be useful, but WITHOUT
  12    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
  14    more details.
  15
  16    You should have received a copy of the GNU General Public License along with
  17    this program; if not, write to the Free Software Foundation, Inc., 51
  18    Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
  19  */
  20
  21 #include <math/linreg/linreg.h>
  22 #include <math/coefficient.h>
  23 #include <gl/xalloc.h>
  24
  25 /*
  26   Predict the value of the dependent variable with the
  27   new set of predictors. PREDICTORS must point to a list
  28   of variables, each of whose values are stored in VALS,
  29   in the same order.
  30  */
  31 double
  32 pspp_linreg_predict (const struct variable **predictors,
  33                      const union value **vals, const void *c_, int n_vals)
  34 {
  35   const pspp_linreg_cache *c = c_;
  36   int i;
  37   int j;
  38   const struct pspp_coeff **found;
  39   const struct pspp_coeff *coe;
  40   double result;
  41   double tmp;
  42
  43   if (predictors == NULL || vals == NULL || c == NULL)
  44     {
  45       return GSL_NAN;
  46     }
  47   if (c->coeff == NULL)
  48     {
  49       /* The stupid model: just guess the mean. */
  50       return c->depvar_mean;
  51     }
  52   found = xnmalloc (c->n_coeffs, sizeof (*found));
  53   *found = c->coeff[0];
  54   result = c->coeff[0]->estimate;       /* Intercept. */
  55
  56   /*
  57      The loops guard against the possibility that the caller passed us
  58      inadequate information, such as too few or too many values, or
  59      a redundant list of variable names.
  60    */
  61   for (j = 0; j < n_vals; j++)
  62     {
  63       coe = pspp_linreg_get_coeff (c, predictors[j], vals[j]);
  64       i = 1;
  65       while (found[i] == coe && i < c->n_coeffs)
  66         {
  67           i++;
  68         }
  69       if (i < c->n_coeffs)
  70         {
  71           found[i] = coe;
  72           tmp = pspp_coeff_get_est (coe);
  73           if (predictors[j]->type == NUMERIC)
  74             {
  75               tmp *= vals[j]->f;
  76             }
  77           result += tmp;
  78         }
  79     }
  80   free (found);
  81
  82   return result;
  83 }
  84
  85 double
  86 pspp_linreg_residual (const struct variable **predictors,
  87                       const union value **vals,
  88                       const union value *obs, const void *c, int n_vals)
  89 {
  90   double pred;
  91   double result;
  92
  93   if (predictors == NULL || vals == NULL || c == NULL || obs == NULL)
  94     {
  95       return GSL_NAN;
  96     }
  97   pred = pspp_linreg_predict (predictors, vals, c, n_vals);
  98
  99   result = gsl_isnan (pred) ? GSL_NAN : (obs->f - pred);
 100   return result;
 101 }