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     {
 116       char buf[5];
 117       snprintf (buf,5,"%g", (upper + lower) / 2.0);
 118       draw_tick (ch, TICK_ABSCISSA,
 119                 x_pos + width / 2.0, buf);
 120     }
 121   }
 122 }
 123
 124
 125
 126 /* Plots a histogram of the data in HIST with the given LABEL.
 127    Labels the histogram with each of N, MEAN, and STDDEV that is
 128    not SYSMIS.  If all three are not SYSMIS and SHOW_NORMAL is
 129    true, also draws a normal curve on the histogram. */
 130 void
 131 histogram_plot (const struct histogram *hist,
 132                 const char *label,
 133                 double n, double mean, double stddev,
 134                 bool show_normal)
 135 {
 136   int i;
 137   int bins;
 138
 139   struct chart *ch = chart_create ();
 140
 141   chart_write_title (ch, _("HISTOGRAM"));
 142
 143   chart_write_ylabel (ch, _("Frequency"));
 144   chart_write_xlabel (ch, label);
 145
 146   if ( ! hist ) /* If this happens, probably all values are SYSMIS */
 147     {
 148       chart_submit (ch);
 149       return;
 150     }
 151   else
 152     {
 153       bins = gsl_histogram_bins (hist->gsl_hist);
 154     }
 155
 156   chart_write_yscale (ch, 0, gsl_histogram_max_val (hist->gsl_hist), 5);
 157
 158   for ( i = 0 ; i < bins ; ++i )
 159     hist_draw_bar (ch, hist, i);
 160
 161   histogram_write_legend (ch, n, mean, stddev);
 162
 163   if (show_normal && n != SYSMIS && mean != SYSMIS && stddev != SYSMIS)
 164     {
 165       /* Draw the normal curve */
 166
 167       double d ;
 168       double x_min, x_max, not_used ;
 169       double abscissa_scale ;
 170       double ordinate_scale ;
 171       double range ;
 172
 173       gsl_histogram_get_range (hist->gsl_hist, 0, &x_min, &not_used);
 174       range = not_used - x_min;
 175       gsl_histogram_get_range (hist->gsl_hist, bins - 1, &not_used, &x_max);
 176
 177       abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
 178       ordinate_scale = (ch->data_top - ch->data_bottom) /
 179         gsl_histogram_max_val (hist->gsl_hist) ;
 180
 181       pl_move_r (ch->lp, ch->data_left, ch->data_bottom);
 182       for ( d = ch->data_left;
 183             d <= ch->data_right ;
 184             d += (ch->data_right - ch->data_left) / 100.0)
 185         {
 186           const double x = (d - ch->data_left) / abscissa_scale + x_min ;
 187           const double y = n * range *
 188             gsl_ran_gaussian_pdf (x - mean, stddev);
 189
 190           pl_fcont_r (ch->lp,  d,  ch->data_bottom  + y * ordinate_scale);
 191
 192         }
 193       pl_endpath_r (ch->lp);
 194     }
 195
 196   chart_submit (ch);
 197 }
 198
 199