src/output/charts/plot-hist.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2004, 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
  18 #include <config.h>
  19
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <gsl/gsl_histogram.h>
  23 #include <gsl/gsl_randist.h>
  24 #include <assert.h>
  25
  26 #include <output/charts/plot-hist.h>
  27 #include <output/charts/plot-chart.h>
  28 #include <output/chart-provider.h>
  29
  30 #include <data/variable.h>
  31 #include <libpspp/cast.h>
  32 #include <libpspp/hash.h>
  33 #include <output/chart.h>
  34 #include <math/histogram.h>
  35 #include <math/moments.h>
  36
  37 #include "gettext.h"
  38 #define _(msgid) gettext (msgid)
  39
  40 static const struct chart_class histogram_chart_class;
  41
  42 /* Write the legend of the chart */
  43 static void
  44 histogram_write_legend (cairo_t *cr, const struct chart_geometry *geom,
  45                         double n, double mean, double stddev)
  46 {
  47   double y = geom->data_bottom;
  48   cairo_save (cr);
  49
  50   if (n != SYSMIS)
  51     {
  52       char *buf = xasprintf ("N = %.2f", n);
  53       cairo_move_to (cr, geom->legend_left, y);
  54       chart_label (cr, 'l', 'b', geom->font_size, buf);
  55       y += geom->font_size * 1.5;
  56       free (buf);
  57     }
  58
  59   if (mean != SYSMIS)
  60     {
  61       char *buf = xasprintf ("Mean = %.1f", mean);
  62       cairo_move_to (cr,geom->legend_left, y);
  63       chart_label (cr, 'l', 'b', geom->font_size, buf);
  64       y += geom->font_size * 1.5;
  65       free (buf);
  66     }
  67
  68   if (stddev != SYSMIS)
  69     {
  70       char *buf = xasprintf ("Std. Dev = %.2f", stddev);
  71       cairo_move_to (cr, geom->legend_left, y);
  72       chart_label (cr, 'l', 'b', geom->font_size, buf);
  73       free (buf);
  74     }
  75
  76   cairo_restore (cr);
  77 }
  78
  79 static void
  80 hist_draw_bar (cairo_t *cr, const struct chart_geometry *geom,
  81                const gsl_histogram *h, int bar)
  82 {
  83   double upper;
  84   double lower;
  85   double height;
  86
  87   const size_t bins = gsl_histogram_bins (h);
  88   const double x_pos = (geom->data_right - geom->data_left) * bar / (double) bins ;
  89   const double width = (geom->data_right - geom->data_left) / (double) bins ;
  90
  91   assert ( 0 == gsl_histogram_get_range (h, bar, &lower, &upper));
  92
  93   assert ( upper >= lower);
  94
  95   height = gsl_histogram_get (h, bar) *
  96     (geom->data_top - geom->data_bottom) / gsl_histogram_max_val (h);
  97
  98   cairo_rectangle (cr, geom->data_left + x_pos, geom->data_bottom,
  99                    width, height);
 100   cairo_save (cr);
 101   cairo_set_source_rgb (cr,
 102                         geom->fill_colour.red / 255.0,
 103                         geom->fill_colour.green / 255.0,
 104                         geom->fill_colour.blue / 255.0);
 105   cairo_fill_preserve (cr);
 106   cairo_restore (cr);
 107   cairo_stroke (cr);
 108
 109   draw_tick (cr, geom, TICK_ABSCISSA,
 110              x_pos + width / 2.0, "%g", (upper + lower) / 2.0);
 111 }
 112
 113 struct histogram_chart
 114   {
 115     struct chart chart;
 116     gsl_histogram *gsl_hist;
 117     char *label;
 118     double n;
 119     double mean;
 120     double stddev;
 121     bool show_normal;
 122   };
 123
 124 /* Plots a histogram of the data in HIST with the given LABEL.
 125    Labels the histogram with each of N, MEAN, and STDDEV that is
 126    not SYSMIS.  If all three are not SYSMIS and SHOW_NORMAL is
 127    true, also draws a normal curve on the histogram. */
 128 struct chart *
 129 histogram_chart_create (const struct histogram *hist, const char *label,
 130                         double n, double mean, double stddev,
 131                         bool show_normal)
 132 {
 133   struct histogram_chart *h;
 134
 135   h = xmalloc (sizeof *h);
 136   chart_init (&h->chart, &histogram_chart_class);
 137   h->gsl_hist = hist->gsl_hist ? gsl_histogram_clone (hist->gsl_hist) : NULL;
 138   h->label = xstrdup (label);
 139   h->n = n;
 140   h->mean = mean;
 141   h->stddev = stddev;
 142   h->show_normal = show_normal;
 143   return &h->chart;
 144 }
 145
 146 static void
 147 histogram_chart_draw (const struct chart *chart, cairo_t *cr,
 148                       struct chart_geometry *geom)
 149 {
 150   struct histogram_chart *h = UP_CAST (chart, struct histogram_chart, chart);
 151   int i;
 152   int bins;
 153
 154   chart_write_title (cr, geom, _("HISTOGRAM"));
 155
 156   chart_write_ylabel (cr, geom, _("Frequency"));
 157   chart_write_xlabel (cr, geom, h->label);
 158
 159   if (h->gsl_hist == NULL)
 160     {
 161       /* Probably all values are SYSMIS. */
 162       return;
 163     }
 164
 165   bins = gsl_histogram_bins (h->gsl_hist);
 166
 167   chart_write_yscale (cr, geom, 0, gsl_histogram_max_val (h->gsl_hist), 5);
 168
 169   for (i = 0; i < bins; i++)
 170     hist_draw_bar (cr, geom, h->gsl_hist, i);
 171
 172   histogram_write_legend (cr, geom, h->n, h->mean, h->stddev);
 173
 174   if (h->show_normal
 175       && h->n != SYSMIS && h->mean != SYSMIS && h->stddev != SYSMIS)
 176     {
 177       /* Draw the normal curve */
 178       double d;
 179       double x_min, x_max, not_used;
 180       double abscissa_scale;
 181       double ordinate_scale;
 182       double range;
 183
 184       gsl_histogram_get_range (h->gsl_hist, 0, &x_min, &not_used);
 185       range = not_used - x_min;
 186       gsl_histogram_get_range (h->gsl_hist, bins - 1, &not_used, &x_max);
 187
 188       abscissa_scale = (geom->data_right - geom->data_left) / (x_max - x_min);
 189       ordinate_scale = (geom->data_top - geom->data_bottom) /
 190         gsl_histogram_max_val (h->gsl_hist);
 191
 192       cairo_move_to (cr, geom->data_left, geom->data_bottom);
 193       for (d = geom->data_left;
 194            d <= geom->data_right;
 195            d += (geom->data_right - geom->data_left) / 100.0)
 196         {
 197           const double x = (d - geom->data_left) / abscissa_scale + x_min;
 198           const double y = h->n * range *
 199             gsl_ran_gaussian_pdf (x - h->mean, h->stddev);
 200
 201           cairo_line_to (cr, d, geom->data_bottom  + y * ordinate_scale);
 202
 203         }
 204       cairo_stroke (cr);
 205     }
 206 }
 207
 208
 209 static void
 210 histogram_chart_destroy (struct chart *chart)
 211 {
 212   struct histogram_chart *h = UP_CAST (chart, struct histogram_chart, chart);
 213   if (h->gsl_hist != NULL)
 214     gsl_histogram_free (h->gsl_hist);
 215   free (h->label);
 216   free (h);
 217 }
 218
 219 static const struct chart_class histogram_chart_class =
 220   {
 221     histogram_chart_draw,
 222     histogram_chart_destroy
 223   };