src/output/charts/plot-hist.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright  (C) 2004 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   char buf[100];
  44
  45   if (!ch)
  46     return ;
  47
  48   pl_savestate_r (ch->lp);
  49
  50   sprintf (buf, "N = %.2f", n);
  51   pl_move_r (ch->lp, ch->legend_left, ch->data_bottom);
  52   pl_alabel_r (ch->lp, 0, 'b', buf);
  53
  54   sprintf (buf, "Mean = %.1f", mean);
  55   pl_fmove_r (ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5);
  56   pl_alabel_r (ch->lp, 0, 'b', buf);
  57
  58   sprintf (buf, "Std. Dev = %.2f", stddev);
  59   pl_fmove_r (ch->lp, ch->legend_left, ch->data_bottom + ch->font_size * 1.5 * 2);
  60   pl_alabel_r (ch->lp, 0, 'b', buf);
  61
  62   pl_restorestate_r (ch->lp);
  63 }
  64
  65 static void hist_draw_bar (struct chart *ch, const struct histogram *hist, int bar);
  66
  67
  68 static void
  69 hist_draw_bar (struct chart *ch, const struct histogram *hist, int bar)
  70 {
  71   if (!ch)
  72     return ;
  73
  74   {
  75     double upper;
  76     double lower;
  77     double height;
  78
  79     const size_t bins = gsl_histogram_bins (hist->gsl_hist);
  80     const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
  81     const double width = (ch->data_right - ch->data_left) / (double) bins ;
  82
  83     assert ( 0 == gsl_histogram_get_range (hist->gsl_hist, bar, &lower, &upper));
  84
  85     assert ( upper >= lower);
  86
  87     height = gsl_histogram_get (hist->gsl_hist, bar) *
  88      (ch->data_top - ch->data_bottom) / gsl_histogram_max_val (hist->gsl_hist);
  89
  90     pl_savestate_r (ch->lp);
  91     pl_move_r (ch->lp,ch->data_left, ch->data_bottom);
  92     pl_fillcolorname_r (ch->lp, ch->fill_colour);
  93     pl_filltype_r (ch->lp,1);
  94
  95
  96     pl_fboxrel_r (ch->lp,
  97                  x_pos, 0,
  98                  x_pos + width, height);
  99
 100     pl_restorestate_r (ch->lp);
 101
 102     {
 103       char buf[5];
 104       snprintf (buf,5,"%g", (upper + lower) / 2.0);
 105       draw_tick (ch, TICK_ABSCISSA,
 106                 x_pos + width / 2.0, buf);
 107     }
 108   }
 109 }
 110
 111
 112
 113 void
 114 histogram_plot (const struct histogram *hist,
 115                 const char *label,
 116                 const struct moments1 *m)
 117 {
 118   double mean, var, n;
 119
 120   moments1_calculate (m, &n, &mean, &var, NULL,  NULL);
 121
 122   histogram_plot_n (hist, label, n, mean, sqrt(var), m);
 123 }
 124
 125
 126 /* This function is deprecated.  Don't use it in new code */
 127 void
 128 histogram_plot_n (const struct histogram *hist,
 129                   const char *label,
 130                   double n, double mean, double stddev,
 131                   bool show_normal)
 132 {
 133   int i;
 134   int bins;
 135
 136   struct chart *ch = chart_create ();
 137
 138   chart_write_title (ch, _("HISTOGRAM"));
 139
 140   chart_write_ylabel (ch, _("Frequency"));
 141   chart_write_xlabel (ch, label);
 142
 143   if ( ! hist ) /* If this happens, probably all values are SYSMIS */
 144     {
 145       chart_submit (ch);
 146       return;
 147     }
 148   else
 149     {
 150       bins = gsl_histogram_bins (hist->gsl_hist);
 151     }
 152
 153   chart_write_yscale (ch, 0, gsl_histogram_max_val (hist->gsl_hist), 5);
 154
 155   for ( i = 0 ; i < bins ; ++i )
 156     hist_draw_bar (ch, hist, i);
 157
 158   histogram_write_legend (ch, n, mean, stddev);
 159
 160   if (show_normal)
 161     {
 162       /* Draw the normal curve */
 163
 164       double d ;
 165       double x_min, x_max, not_used ;
 166       double abscissa_scale ;
 167       double ordinate_scale ;
 168       double range ;
 169
 170       gsl_histogram_get_range (hist->gsl_hist, 0, &x_min, &not_used);
 171       range = not_used - x_min;
 172       gsl_histogram_get_range (hist->gsl_hist, bins - 1, &not_used, &x_max);
 173
 174       abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
 175       ordinate_scale = (ch->data_top - ch->data_bottom) /
 176         gsl_histogram_max_val (hist->gsl_hist) ;
 177
 178       pl_move_r (ch->lp, ch->data_left, ch->data_bottom);
 179       for ( d = ch->data_left;
 180             d <= ch->data_right ;
 181             d += (ch->data_right - ch->data_left) / 100.0)
 182         {
 183           const double x = (d - ch->data_left) / abscissa_scale + x_min ;
 184           const double y = n * range *
 185             gsl_ran_gaussian_pdf (x - mean, stddev);
 186
 187           pl_fcont_r (ch->lp,  d,  ch->data_bottom  + y * ordinate_scale);
 188
 189         }
 190       pl_endpath_r (ch->lp);
 191     }
 192
 193   chart_submit (ch);
 194 }
 195
 196