src/plot-hist.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 2004 Free Software Foundation, Inc.
   3    Written by John Darrington <john@darrington.wattle.id.au>
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18    02110-1301, USA. */
  19
  20 /* If you add/modify any public symbols in this file, don't forget to
  21    change the stubs in dummy-chart.c */
  22
  23 #include <config.h>
  24
  25 #include <stdio.h>
  26 #include <plot.h>
  27 #include <math.h>
  28 #include <gsl/gsl_histogram.h>
  29 #include <gsl/gsl_randist.h>
  30 #include <assert.h>
  31 #include "hash.h"
  32 #include "var.h"
  33 #include "chart.h"
  34
  35 #include "gettext.h"
  36 #define _(msgid) gettext (msgid)
  37
  38 /* Write the legend of the chart */
  39 void
  40 histogram_write_legend(struct chart *ch, const struct normal_curve *norm)
  41 {
  42   char buf[100];
  43   if ( !ch )
  44           return ;
  45
  46   pl_savestate_r(ch->lp);
  47
  48   sprintf(buf,"N = %.2f",norm->N);
  49   pl_move_r(ch->lp, ch->legend_left, ch->data_bottom);
  50   pl_alabel_r(ch->lp,0,'b',buf);
  51
  52   sprintf(buf,"Mean = %.1f",norm->mean);
  53   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5);
  54   pl_alabel_r(ch->lp,0,'b',buf);
  55
  56   sprintf(buf,"Std. Dev = %.2f",norm->stddev);
  57   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5 * 2);
  58   pl_alabel_r(ch->lp,0,'b',buf);
  59
  60   pl_restorestate_r(ch->lp);
  61 }
  62
  63 static void hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar);
  64
  65
  66 static void
  67 hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar)
  68 {
  69   if ( !ch )
  70     return ;
  71
  72
  73   {
  74     double upper;
  75     double lower;
  76     double height;
  77
  78     const size_t bins = gsl_histogram_bins(hist);
  79     const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
  80     const double width = (ch->data_right - ch->data_left) / (double) bins ;
  81
  82
  83     assert ( 0 == gsl_histogram_get_range(hist, bar, &lower, &upper));
  84
  85     assert( upper >= lower);
  86
  87     height = gsl_histogram_get(hist, bar) *
  88       (ch->data_top - ch->data_bottom) / gsl_histogram_max_val(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
 114 void
 115 histogram_plot(const gsl_histogram *hist,
 116                const char *factorname,
 117                const struct normal_curve *norm, short show_normal)
 118 {
 119   int i;
 120   int bins;
 121
 122   struct chart *ch;
 123
 124   bins = gsl_histogram_bins(hist);
 125
 126   ch = chart_create();
 127   chart_write_title(ch, _("HISTOGRAM"));
 128
 129   chart_write_ylabel(ch, _("Frequency"));
 130   chart_write_xlabel(ch, factorname);
 131
 132   chart_write_yscale(ch, 0, gsl_histogram_max_val(hist), 5);
 133
 134   for ( i = 0 ; i < bins ; ++i )
 135       hist_draw_bar(ch, hist, i);
 136
 137   histogram_write_legend(ch, norm);
 138
 139   if ( show_normal  )
 140   {
 141     /* Draw the normal curve */
 142
 143     double d ;
 144     double x_min, x_max, not_used ;
 145     double abscissa_scale ;
 146     double ordinate_scale ;
 147     double range ;
 148
 149     gsl_histogram_get_range(hist, 0, &x_min, &not_used);
 150     range = not_used - x_min;
 151     gsl_histogram_get_range(hist, bins - 1, &not_used, &x_max);
 152     assert(range == x_max - not_used);
 153
 154     abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
 155     ordinate_scale = (ch->data_top - ch->data_bottom) /
 156       gsl_histogram_max_val(hist) ;
 157
 158     pl_move_r(ch->lp, ch->data_left, ch->data_bottom);
 159     for( d = ch->data_left;
 160          d <= ch->data_right ;
 161          d += (ch->data_right - ch->data_left) / 100.0)
 162       {
 163         const double x = (d - ch->data_left) / abscissa_scale + x_min ;
 164         const double y = norm->N * range *
 165           gsl_ran_gaussian_pdf(x - norm->mean, norm->stddev);
 166
 167         pl_fcont_r(ch->lp,  d,  ch->data_bottom  + y * ordinate_scale);
 168
 169       }
 170     pl_endpath_r(ch->lp);
 171
 172   }
 173   chart_submit(ch);
 174 }
 175