Added code to calculate the histogram limits
[pspp-builds.git] / 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