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 <plot.h>
  22 #include <math.h>
  23 #include <gsl/gsl_histogram.h>
  24 #include <gsl/gsl_randist.h>
  25 #include <assert.h>
  26
  27 #include <output/charts/plot-hist.h>
  28 #include <output/charts/plot-chart.h>
  29
  30 #include <data/variable.h>
  31 #include <libpspp/hash.h>
  32 #include <output/chart.h>
  33 #include <math/histogram.h>
  34 #include <math/moments.h>
  35
  36 #include "gettext.h"
  37 #define _(msgid) gettext (msgid)
  38
  39 /* Write the legend of the chart */
  40 static void
  41 histogram_write_legend (struct chart *ch, double n, double mean, double stddev)
  42 {
  43   double y;
  44   char buf[100];
  45
  46   if (!ch)
  47     return ;
  48
  49   y = ch->data_bottom;
  50   pl_savestate_r (ch->lp);
  51
  52   if (n != SYSMIS)
  53     {
  54       sprintf (buf, "N = %.2f", n);
  55       pl_move_r (ch->lp, ch->legend_left, y);
  56       pl_alabel_r (ch->lp, 0, 'b', buf);
  57       y += ch->font_size * 1.5;
  58     }
  59
  60   if (mean != SYSMIS)
  61     {
  62       sprintf (buf, "Mean = %.1f", mean);
  63       pl_fmove_r (ch->lp,ch->legend_left, y);
  64       pl_alabel_r (ch->lp, 0, 'b', buf);
  65       y += ch->font_size * 1.5;
  66     }
  67
  68   if (stddev != SYSMIS)
  69     {
  70       sprintf (buf, "Std. Dev = %.2f", stddev);
  71       pl_fmove_r (ch->lp, ch->legend_left, y);
  72       pl_alabel_r (ch->lp, 0, 'b', buf);
  73     }
  74
  75   pl_restorestate_r (ch->lp);
  76 }
  77
  78 static void hist_draw_bar (struct chart *ch, const struct histogram *hist, int bar);
  79
  80
  81 static void
  82 hist_draw_bar (struct chart *ch, const struct histogram *hist, int bar)
  83 {
  84   if (!ch)
  85     return ;
  86
  87   {
  88     double upper;
  89     double lower;
  90     double height;
  91
  92     const size_t bins = gsl_histogram_bins (hist->gsl_hist);
  93     const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
  94     const double width = (ch->data_right - ch->data_left) / (double) bins ;
  95
  96     assert ( 0 == gsl_histogram_get_range (hist->gsl_hist, bar, &lower, &upper));
  97
  98     assert ( upper >= lower);
  99
 100     height = gsl_histogram_get (hist->gsl_hist, bar) *
 101      (ch->data_top - ch->data_bottom) / gsl_histogram_max_val (hist->gsl_hist);
 102
 103     pl_savestate_r (ch->lp);
 104     pl_move_r (ch->lp,ch->data_left, ch->data_bottom);
 105     pl_fillcolorname_r (ch->lp, ch->fill_colour);
 106     pl_filltype_r (ch->lp,1);
 107
 108
 109     pl_fboxrel_r (ch->lp,
 110                  x_pos, 0,
 111                  x_pos + width, height);
 112
 113     pl_restorestate_r (ch->lp);
 114
 115     draw_tick (ch, TICK_ABSCISSA,
 116                x_pos + width / 2.0, "%g", (upper + lower) / 2.0);
 117   }
 118 }
 119
 120
 121
 122 /* Plots a histogram of the data in HIST with the given LABEL.
 123    Labels the histogram with each of N, MEAN, and STDDEV that is
 124    not SYSMIS.  If all three are not SYSMIS and SHOW_NORMAL is
 125    true, also draws a normal curve on the histogram. */
 126 void
 127 histogram_plot (const struct histogram *hist,
 128                 const char *label,
 129                 double n, double mean, double stddev,
 130                 bool show_normal)
 131 {
 132   int i;
 133   int bins;
 134
 135   struct chart *ch = chart_create ();
 136
 137   chart_write_title (ch, _("HISTOGRAM"));
 138
 139   chart_write_ylabel (ch, _("Frequency"));
 140   chart_write_xlabel (ch, label);
 141
 142   if ( ! hist ) /* If this happens, probably all values are SYSMIS */
 143     {
 144       chart_submit (ch);
 145       return;
 146     }
 147   else
 148     {
 149       bins = gsl_histogram_bins (hist->gsl_hist);
 150     }
 151
 152   chart_write_yscale (ch, 0, gsl_histogram_max_val (hist->gsl_hist), 5);
 153
 154   for ( i = 0 ; i < bins ; ++i )
 155     hist_draw_bar (ch, hist, i);
 156
 157   histogram_write_legend (ch, n, mean, stddev);
 158
 159   if (show_normal && n != SYSMIS && mean != SYSMIS && stddev != SYSMIS)
 160     {
 161       /* Draw the normal curve */
 162
 163       double d ;
 164       double x_min, x_max, not_used ;
 165       double abscissa_scale ;
 166       double ordinate_scale ;
 167       double range ;
 168
 169       gsl_histogram_get_range (hist->gsl_hist, 0, &x_min, &not_used);
 170       range = not_used - x_min;
 171       gsl_histogram_get_range (hist->gsl_hist, bins - 1, &not_used, &x_max);
 172
 173       abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
 174       ordinate_scale = (ch->data_top - ch->data_bottom) /
 175         gsl_histogram_max_val (hist->gsl_hist) ;
 176
 177       pl_move_r (ch->lp, ch->data_left, ch->data_bottom);
 178       for ( d = ch->data_left;
 179             d <= ch->data_right ;
 180             d += (ch->data_right - ch->data_left) / 100.0)
 181         {
 182           const double x = (d - ch->data_left) / abscissa_scale + x_min ;
 183           const double y = n * range *
 184             gsl_ran_gaussian_pdf (x - mean, stddev);
 185
 186           pl_fcont_r (ch->lp,  d,  ch->data_bottom  + y * ordinate_scale);
 187
 188         }
 189       pl_endpath_r (ch->lp);
 190     }
 191
 192   chart_submit (ch);
 193 }
 194
 195