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