src/output/charts/np-plot.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2004, 2008, 2009 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 <output/charts/np-plot.h>
  20
  21 #include <gsl/gsl_cdf.h>
  22
  23 #include <data/casereader.h>
  24 #include <data/casewriter.h>
  25 #include <libpspp/message.h>
  26 #include <math/np.h>
  27 #include <output/chart-provider.h>
  28 #include <output/charts/cartesian.h>
  29 #include <output/charts/plot-chart.h>
  30
  31 #include "gl/minmax.h"
  32
  33 #include "gettext.h"
  34 #define _(msgid) gettext (msgid)
  35
  36 /* An NP or DNP plot. */
  37 struct np_plot_chart
  38   {
  39     struct chart chart;
  40     char *label;
  41     struct casereader *data;
  42
  43     /* Copied directly from struct np. */
  44     double y_min, y_max;
  45     double dns_min, dns_max;
  46
  47     /* Calculated. */
  48     double slope, intercept;
  49     double y_first, y_last;
  50     double x_lower, x_upper;
  51     double slack;
  52   };
  53
  54 static const struct chart_class np_plot_chart_class;
  55 static const struct chart_class dnp_plot_chart_class;
  56
  57 static struct chart *
  58 make_np_plot (const struct chart_class *class,
  59               const struct np *np, const struct casereader *reader,
  60               const char *label)
  61 {
  62   struct np_plot_chart *npp;
  63
  64   if (np->n < 1.0)
  65     return NULL;
  66
  67   npp = xmalloc (sizeof *npp);
  68   chart_init (&npp->chart, class);
  69   npp->label = xstrdup (label);
  70   npp->data = casereader_clone (reader);
  71   npp->y_min = np->y_min;
  72   npp->y_max = np->y_max;
  73   npp->dns_min = np->dns_min;
  74   npp->dns_max = np->dns_max;
  75
  76   /* Slope and intercept of the ideal normal probability line. */
  77   npp->slope = 1.0 / np->stddev;
  78   npp->intercept = -np->mean / np->stddev;
  79
  80   npp->y_first = gsl_cdf_ugaussian_Pinv (1 / (np->n + 1));
  81   npp->y_last = gsl_cdf_ugaussian_Pinv (np->n / (np->n + 1));
  82
  83   /* Need to make sure that both the scatter plot and the ideal fit into the
  84      plot. */
  85   npp->x_lower = MIN (np->y_min, (npp->y_first - npp->intercept) / npp->slope);
  86   npp->x_upper = MAX (np->y_max, (npp->y_last  - npp->intercept) / npp->slope);
  87   npp->slack = (npp->x_upper - npp->x_lower) * 0.05;
  88
  89   return &npp->chart;
  90 }
  91
  92 /* Creates and returns a normal probability plot corresponding to
  93    the calculations in NP and the data in READER, and label the
  94    plot with LABEL.  The data in READER must have Y-values in
  95    value index NP_IDX_Y and NS-values in value index NP_IDX_NS.
  96
  97    Returns a null pointer if the data set is empty.
  98
  99    The caller retains ownership of NP and READER. */
 100 struct chart *
 101 np_plot_create (const struct np *np, const struct casereader *reader,
 102                 const char *label)
 103 {
 104   return make_np_plot (&np_plot_chart_class, np, reader, label);
 105 }
 106
 107 /* Creates and returns a detrended normal probability plot
 108    corresponding to the calculations in NP and the data in
 109    READER, and label the plot with LABEL.  The data in READER
 110    must have Y-values in value index NP_IDX_Y and DNS-values in
 111    value index NP_IDX_DNS.
 112
 113    Returns a null pointer if the data set is empty.
 114
 115    The caller retains ownership of NP and READER. */
 116 struct chart *
 117 dnp_plot_create (const struct np *np, const struct casereader *reader,
 118                  const char *label)
 119 {
 120   return make_np_plot (&dnp_plot_chart_class, np, reader, label);
 121 }
 122
 123 static void
 124 np_plot_chart_draw (const struct chart *chart, plPlotter *lp)
 125 {
 126   const struct np_plot_chart *npp = (struct np_plot_chart *) chart;
 127   struct chart_geometry geom;
 128   struct casereader *data;
 129   struct ccase *c;
 130
 131   chart_geometry_init (lp, &geom);
 132   chart_write_title (lp, &geom, _("Normal Q-Q Plot of %s"), npp->label);
 133   chart_write_xlabel (lp, &geom, _("Observed Value"));
 134   chart_write_ylabel (lp, &geom, _("Expected Normal"));
 135   chart_write_xscale (lp, &geom,
 136                       npp->x_lower - npp->slack,
 137                       npp->x_upper + npp->slack, 5);
 138   chart_write_yscale (lp, &geom, npp->y_first, npp->y_last, 5);
 139
 140   data = casereader_clone (npp->data);
 141   for (; (c = casereader_read (data)) != NULL; case_unref (c))
 142     chart_datum (lp, &geom, 0,
 143                  case_data_idx (c, NP_IDX_Y)->f,
 144                  case_data_idx (c, NP_IDX_NS)->f);
 145   casereader_destroy (data);
 146
 147   chart_line (lp, &geom, npp->slope, npp->intercept,
 148               npp->y_first, npp->y_last, CHART_DIM_Y);
 149
 150   chart_geometry_free (lp);
 151 }
 152
 153 static void
 154 dnp_plot_chart_draw (const struct chart *chart, plPlotter *lp)
 155 {
 156   const struct np_plot_chart *dnpp = (struct np_plot_chart *) chart;
 157   struct chart_geometry geom;
 158   struct casereader *data;
 159   struct ccase *c;
 160
 161   chart_geometry_init (lp, &geom);
 162   chart_write_title (lp, &geom, _("Detrended Normal Q-Q Plot of %s"),
 163                      dnpp->label);
 164   chart_write_xlabel (lp, &geom, _("Observed Value"));
 165   chart_write_ylabel (lp, &geom, _("Dev from Normal"));
 166   chart_write_xscale (lp, &geom, dnpp->y_min, dnpp->y_max, 5);
 167   chart_write_yscale (lp, &geom, dnpp->dns_min, dnpp->dns_max, 5);
 168
 169   data = casereader_clone (dnpp->data);
 170   for (; (c = casereader_read (data)) != NULL; case_unref (c))
 171     chart_datum (lp, &geom, 0, case_data_idx (c, NP_IDX_Y)->f,
 172                  case_data_idx (c, NP_IDX_DNS)->f);
 173   casereader_destroy (data);
 174
 175   chart_line (lp, &geom, 0, 0, dnpp->y_min, dnpp->y_max, CHART_DIM_X);
 176
 177   chart_geometry_free (lp);
 178 }
 179
 180 static void
 181 np_plot_chart_destroy (struct chart *chart)
 182 {
 183   struct np_plot_chart *npp = (struct np_plot_chart *) chart;
 184
 185   casereader_destroy (npp->data);
 186   free (npp->label);
 187   free (npp);
 188 }
 189
 190 static const struct chart_class np_plot_chart_class =
 191   {
 192     np_plot_chart_draw,
 193     np_plot_chart_destroy
 194   };
 195
 196 static const struct chart_class dnp_plot_chart_class =
 197   {
 198     dnp_plot_chart_draw,
 199     np_plot_chart_destroy
 200   };