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 <stdio.h>
  21 #include <plot.h>
  22 #include <stdarg.h>
  23 #include <math.h>
  24 #include "hash.h"
  25
  26 #include "var.h"
  27 #include "chart.h"
  28
  29 /* Number of bins in which to divide data */
  30 #define BINS 15
  31
  32 /* The approximate no of ticks on the y axis */
  33 #define YTICKS 10
  34
  35 #define M_PI ( 22.0 / 7.0 )
  36
  37
  38 static double gaussian(double x, double mu, double sigma ) ;
  39
  40
  41 static double
  42 gaussian(double x, double mu, double sigma )
  43 {
  44   return (exp( - (( x - mu )* (x - mu) / (2.0 * sigma * sigma)  ))
  45     / ( sigma * sqrt( M_PI * 2.0) ))  ;
  46 }
  47
  48
  49 /* Adjust tick to be a sensible value */
  50 void adjust_tick(double *tick);
  51
  52
  53 /* Write the legend of the chart */
  54 static void
  55 write_legend(struct chart *ch, struct normal_curve *norm)
  56 {
  57   char buf[100];
  58   pl_savestate_r(ch->lp);
  59
  60   sprintf(buf,"N = %.2f",norm->N);
  61   pl_move_r(ch->lp, ch->legend_left, ch->data_bottom);
  62   pl_alabel_r(ch->lp,0,'b',buf);
  63
  64   sprintf(buf,"Mean = %.1f",norm->mean);
  65   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5);
  66   pl_alabel_r(ch->lp,0,'b',buf);
  67
  68   sprintf(buf,"Std. Dev = %.2f",norm->stddev);
  69   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5 * 2);
  70   pl_alabel_r(ch->lp,0,'b',buf);
  71
  72   pl_restorestate_r(ch->lp);
  73 }
  74
  75
  76
  77
  78 /* Draw a histogram.
  79    If show_normal is non zero then superimpose a normal curve
  80 */
  81 void
  82 draw_histogram(struct chart *ch,
  83                const struct variable *var,
  84                const char *title,
  85                struct normal_curve *norm,
  86                int show_normal)
  87 {
  88
  89   double d;
  90   int count;
  91
  92   double x_min = DBL_MAX;
  93   double x_max = -DBL_MAX;
  94   double y_min = DBL_MAX;
  95   double y_max = -DBL_MAX;
  96
  97   double y_tick ;
  98
  99
 100   double ordinate_values[BINS];
 101
 102   const struct freq_tab *frq_tab = &var->p.frq.tab ;
 103
 104   struct hsh_iterator hi;
 105   struct hsh_table *fh = frq_tab->data;
 106   struct freq *frq;
 107
 108   double interval_size = fabs(ch->data_right - ch->data_left) / ( BINS );
 109
 110   /* Find out the extremes of the x value
 111      FIXME: These don't need to be calculated here, since the
 112      calling routine should know them */
 113
 114   for ( frq = hsh_first(fh,&hi); frq != 0; frq = hsh_next(fh,&hi) )
 115     {
 116       if ( frq->v.f < x_min ) x_min = frq->v.f ;
 117       if ( frq->v.f > x_max ) x_max = frq->v.f ;
 118     }
 119
 120   double x_interval = fabs(x_max - x_min) / ( BINS - 1);
 121
 122
 123   double abscissa_scale =
 124     fabs( (ch->data_right -  ch->data_left) / (x_max - x_min) ) ;
 125
 126
 127   /* Find out the bin values */
 128   for ( count = 0, d = x_min; d <= x_max ; d += x_interval, ++count )
 129     {
 130
 131       double y = 0;
 132
 133       for ( frq = hsh_first(fh,&hi); frq != 0; frq = hsh_next(fh,&hi) )
 134         {
 135           if ( frq->v.f >= d && frq->v.f < d + x_interval )
 136             y += frq->c;
 137         }
 138
 139       ordinate_values[count] = y ;
 140
 141       if ( y > y_max ) y_max = y ;
 142       if ( y < y_min ) y_min = y;
 143     }
 144
 145   y_tick = ( y_max - y_min ) / (double) (YTICKS - 1) ;
 146
 147   adjust_tick(&y_tick);
 148
 149   y_min = floor( y_min / y_tick ) * y_tick ;
 150   y_max = ceil( y_max / y_tick ) * y_tick ;
 151
 152   double ordinate_scale =
 153     fabs(ch->data_top -  ch->data_bottom) / fabs(y_max - y_min) ;
 154
 155
 156   /* Move to data bottom-left */
 157   pl_move_r(ch->lp, ch->data_left, ch->data_bottom);
 158
 159   pl_savestate_r(ch->lp);
 160   pl_fillcolorname_r(ch->lp, ch->fill_colour);
 161   pl_filltype_r(ch->lp,1);
 162
 163   /* Draw the histogram */
 164   for ( count = 0, d = x_min; d <= x_max ; d += x_interval, ++count )
 165     {
 166       const double x = count * interval_size ;
 167       pl_savestate_r(ch->lp);
 168
 169       draw_tick (ch, TICK_ABSCISSA, x + (interval_size / 2.0 ), "%.1f",d);
 170
 171       pl_fboxrel_r(ch->lp, x, 0, x + interval_size,
 172                    ordinate_values[count] * ordinate_scale );
 173
 174       pl_restorestate_r(ch->lp);
 175
 176     }
 177   pl_restorestate_r(ch->lp);
 178
 179   /* Put the y axis on */
 180   for ( d = y_min; d <= y_max ; d += y_tick )
 181     {
 182       draw_tick (ch, TICK_ORDINATE, (d - y_min ) * ordinate_scale, "%g", d);
 183     }
 184
 185   /* Write the abscissa label */
 186   pl_move_r(ch->lp,ch->data_left, ch->abscissa_top);
 187   pl_alabel_r(ch->lp,0,'t',var->label ? var->label : var->name);
 188
 189
 190   /* Write the ordinate label */
 191   pl_savestate_r(ch->lp);
 192   pl_move_r(ch->lp,ch->data_bottom, ch->ordinate_right);
 193   pl_textangle_r(ch->lp,90);
 194   pl_alabel_r(ch->lp,0,0,"Frequency");
 195
 196   pl_restorestate_r(ch->lp);
 197
 198
 199   chart_write_title(ch, title);
 200
 201
 202   /* Write the legend */
 203   write_legend(ch,norm);
 204
 205
 206   if ( show_normal  )
 207   {
 208     /* Draw the normal curve */
 209     double d;
 210
 211     pl_move_r(ch->lp, ch->data_left, ch->data_bottom);
 212     for( d = ch->data_left;
 213          d <= ch->data_right ;
 214          d += (ch->data_right - ch->data_left) / 100.0)
 215       {
 216         const double x = (d  - ch->data_left  - interval_size / 2.0  ) /
 217           abscissa_scale  + x_min ;
 218
 219         pl_fcont_r(ch->lp,  d,
 220                    ch->data_bottom +
 221                    norm->N * gaussian(x, norm->mean, norm->stddev)
 222                    * ordinate_scale);
 223
 224       }
 225     pl_endpath_r(ch->lp);
 226   }
 227
 228 }
 229
 230
 231
 232 double
 233 log10(double x)
 234 {
 235   return log(x) / log(10.0) ;
 236 }
 237
 238
 239 /* Adjust tick to be a sensible value */
 240 void
 241 adjust_tick(double *tick)
 242 {
 243     int i;
 244     const double standard_ticks[] = {1, 2, 5};
 245
 246     const double factor = pow(10,ceil(log10(standard_ticks[0] / *tick))) ;
 247
 248     for (i = 2  ; i >=0 ; --i)
 249       {
 250         if ( *tick > standard_ticks[i] / factor )
 251           {
 252             *tick = standard_ticks[i] / factor ;
 253             break;
 254           }
 255       }
 256
 257   }
 258