1 /* PSPP - computes sample statistics.
2 Copyright (C) 2004 Free Software Foundation, Inc.
3 Written by John Darrington <john@darrington.wattle.id.au>
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.
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.
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
29 /* Number of bins in which to divide data */
32 /* The approximate no of ticks on the y axis */
35 #define M_PI ( 22.0 / 7.0 )
38 static double gaussian(double x, double mu, double sigma ) ;
42 gaussian(double x, double mu, double sigma )
44 return (exp( - (( x - mu )* (x - mu) / (2.0 * sigma * sigma) ))
45 / ( sigma * sqrt( M_PI * 2.0) )) ;
49 /* Adjust tick to be a sensible value */
50 void adjust_tick(double *tick);
53 /* Write the legend of the chart */
55 write_legend(struct chart *ch, struct normal_curve *norm)
58 pl_savestate_r(ch->lp);
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);
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);
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);
72 pl_restorestate_r(ch->lp);
79 If show_normal is non zero then superimpose a normal curve
82 draw_histogram(struct chart *ch,
83 const struct variable *var,
85 struct normal_curve *norm,
92 double x_min = DBL_MAX;
93 double x_max = -DBL_MAX;
94 double y_min = DBL_MAX;
95 double y_max = -DBL_MAX;
100 double ordinate_values[BINS];
102 const struct freq_tab *frq_tab = &var->p.frq.tab ;
104 struct hsh_iterator hi;
105 struct hsh_table *fh = frq_tab->data;
108 double interval_size = fabs(ch->data_right - ch->data_left) / ( BINS );
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 */
114 for ( frq = hsh_first(fh,&hi); frq != 0; frq = hsh_next(fh,&hi) )
116 if ( frq->v.f < x_min ) x_min = frq->v.f ;
117 if ( frq->v.f > x_max ) x_max = frq->v.f ;
120 double x_interval = fabs(x_max - x_min) / ( BINS - 1);
123 double abscissa_scale =
124 fabs( (ch->data_right - ch->data_left) / (x_max - x_min) ) ;
127 /* Find out the bin values */
128 for ( count = 0, d = x_min; d <= x_max ; d += x_interval, ++count )
133 for ( frq = hsh_first(fh,&hi); frq != 0; frq = hsh_next(fh,&hi) )
135 if ( frq->v.f >= d && frq->v.f < d + x_interval )
139 ordinate_values[count] = y ;
141 if ( y > y_max ) y_max = y ;
142 if ( y < y_min ) y_min = y;
145 y_tick = ( y_max - y_min ) / (double) (YTICKS - 1) ;
147 adjust_tick(&y_tick);
149 y_min = floor( y_min / y_tick ) * y_tick ;
150 y_max = ceil( y_max / y_tick ) * y_tick ;
152 double ordinate_scale =
153 fabs(ch->data_top - ch->data_bottom) / fabs(y_max - y_min) ;
156 /* Move to data bottom-left */
157 pl_move_r(ch->lp, ch->data_left, ch->data_bottom);
159 pl_savestate_r(ch->lp);
160 pl_fillcolorname_r(ch->lp, ch->fill_colour);
161 pl_filltype_r(ch->lp,1);
163 /* Draw the histogram */
164 for ( count = 0, d = x_min; d <= x_max ; d += x_interval, ++count )
166 const double x = count * interval_size ;
167 pl_savestate_r(ch->lp);
169 draw_tick (ch, TICK_ABSCISSA, x + (interval_size / 2.0 ), "%.1f",d);
171 pl_fboxrel_r(ch->lp, x, 0, x + interval_size,
172 ordinate_values[count] * ordinate_scale );
174 pl_restorestate_r(ch->lp);
177 pl_restorestate_r(ch->lp);
179 /* Put the y axis on */
180 for ( d = y_min; d <= y_max ; d += y_tick )
182 draw_tick (ch, TICK_ORDINATE, (d - y_min ) * ordinate_scale, "%g", d);
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);
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");
196 pl_restorestate_r(ch->lp);
199 chart_write_title(ch, title);
202 /* Write the legend */
203 write_legend(ch,norm);
208 /* Draw the normal curve */
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)
216 const double x = (d - ch->data_left - interval_size / 2.0 ) /
217 abscissa_scale + x_min ;
219 pl_fcont_r(ch->lp, d,
221 norm->N * gaussian(x, norm->mean, norm->stddev)
225 pl_endpath_r(ch->lp);
235 return log(x) / log(10.0) ;
239 /* Adjust tick to be a sensible value */
241 adjust_tick(double *tick)
244 const double standard_ticks[] = {1, 2, 5};
246 const double factor = pow(10,ceil(log10(standard_ticks[0] / *tick))) ;
248 for (i = 2 ; i >=0 ; --i)
250 if ( *tick > standard_ticks[i] / factor )
252 *tick = standard_ticks[i] / factor ;