src/histogram.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., 59 Temple Place - Suite 330, Boston, MA
  18    02111-1307, USA. */
  19
  20 #include <config.h>
  21
  22 #include <stdio.h>
  23 #include <plot.h>
  24 #include <math.h>
  25 #include <gsl/gsl_histogram.h>
  26 #include <assert.h>
  27
  28 #include "hash.h"
  29
  30 #include "var.h"
  31 #include "chart.h"
  32
  33 /* Number of bins in which to divide data */
  34 #define BINS 7
  35
  36 /* The approximate no of ticks on the y axis */
  37 #define YTICKS 10
  38
  39 #ifndef M_PI
  40 #define M_PI ( 22.0 / 7.0 )
  41 #endif
  42
  43
  44 static double gaussian(double x, double mu, double sigma ) ;
  45
  46
  47 static double
  48 gaussian(double x, double mu, double sigma )
  49 {
  50   return (exp( - (( x - mu )* (x - mu) / (2.0 * sigma * sigma)  ))
  51     / ( sigma * sqrt( M_PI * 2.0) ))  ;
  52 }
  53
  54
  55
  56 /* Write the legend of the chart */
  57 void
  58 histogram_write_legend(struct chart *ch, const struct normal_curve *norm)
  59 {
  60   char buf[100];
  61   pl_savestate_r(ch->lp);
  62
  63   sprintf(buf,"N = %.2f",norm->N);
  64   pl_move_r(ch->lp, ch->legend_left, ch->data_bottom);
  65   pl_alabel_r(ch->lp,0,'b',buf);
  66
  67   sprintf(buf,"Mean = %.1f",norm->mean);
  68   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5);
  69   pl_alabel_r(ch->lp,0,'b',buf);
  70
  71   sprintf(buf,"Std. Dev = %.2f",norm->stddev);
  72   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5 * 2);
  73   pl_alabel_r(ch->lp,0,'b',buf);
  74
  75   pl_restorestate_r(ch->lp);
  76 }
  77
  78 static void hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar);
  79
  80
  81 static void
  82 hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar)
  83 {
  84   double upper;
  85   double lower;
  86   double height;
  87
  88   const size_t bins = gsl_histogram_bins(hist);
  89   const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
  90   const double width = (ch->data_right - ch->data_left) / (double) bins ;
  91
  92
  93   assert ( 0 == gsl_histogram_get_range(hist, bar, &lower, &upper));
  94
  95   assert( upper >= lower);
  96
  97   height = gsl_histogram_get(hist, bar) *
  98     (ch->data_top - ch->data_bottom) / gsl_histogram_max_val(hist);
  99
 100   pl_savestate_r(ch->lp);
 101   pl_move_r(ch->lp,ch->data_left, ch->data_bottom);
 102   pl_fillcolorname_r(ch->lp, ch->fill_colour);
 103   pl_filltype_r(ch->lp,1);
 104
 105
 106   pl_fboxrel_r(ch->lp,
 107                x_pos, 0,
 108                x_pos + width, height);
 109
 110   pl_restorestate_r(ch->lp);
 111
 112   {
 113   char buf[5];
 114   snprintf(buf,5,"%g",(upper + lower) / 2.0);
 115   draw_tick(ch, TICK_ABSCISSA,
 116             x_pos + width / 2.0, buf);
 117   }
 118
 119
 120 }
 121
 122
 123
 124
 125 void
 126 histogram_plot(const gsl_histogram *hist,
 127                const char *factorname,
 128                const struct normal_curve *norm, short show_normal)
 129 {
 130   int i;
 131   int bins;
 132
 133   struct chart ch;
 134
 135   bins = gsl_histogram_bins(hist);
 136
 137   chart_initialise(&ch);
 138   chart_write_title(&ch, _("HISTOGRAM"));
 139
 140   chart_write_ylabel(&ch, _("Frequency"));
 141   chart_write_xlabel(&ch, factorname);
 142
 143   chart_write_yscale(&ch, 0, gsl_histogram_max_val(hist), 5);
 144
 145   for ( i = 0 ; i < bins ; ++i )
 146       hist_draw_bar(&ch, hist, i);
 147
 148   histogram_write_legend(&ch, norm);
 149
 150   if ( show_normal  )
 151   {
 152     /* Draw the normal curve */
 153
 154     double d ;
 155     double x_min, x_max, not_used ;
 156     double abscissa_scale ;
 157     double ordinate_scale ;
 158
 159
 160     gsl_histogram_get_range(hist, 0, &x_min, &not_used);
 161     gsl_histogram_get_range(hist, bins - 1, &x_max, &not_used);
 162
 163     abscissa_scale = (ch.data_right - ch.data_left) / (x_max - x_min);
 164     ordinate_scale = (ch.data_top - ch.data_bottom) /
 165       gsl_histogram_max_val(hist) ;
 166
 167     pl_move_r(ch.lp, ch.data_left, ch.data_bottom);
 168     for( d = ch.data_left;
 169          d <= ch.data_right ;
 170          d += (ch.data_right - ch.data_left) / 100.0)
 171       {
 172         const double x = (d - ch.data_left) / abscissa_scale + x_min ;
 173
 174         pl_fcont_r(ch.lp,  d,
 175                    ch.data_bottom  + norm->N * ordinate_scale *
 176                    gaussian(x, norm->mean, norm->stddev)
 177                    );
 178       }
 179     pl_endpath_r(ch.lp);
 180
 181   }
 182   chart_finalise(&ch);
 183 }
 184