1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2004 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
23 #include <gsl/gsl_histogram.h>
24 #include <gsl/gsl_randist.h>
27 #include <output/charts/plot-hist.h>
28 #include <output/charts/plot-chart.h>
30 #include <data/variable.h>
31 #include <libpspp/hash.h>
32 #include <output/chart.h>
33 #include <math/histogram.h>
34 #include <math/moments.h>
37 #define _(msgid) gettext (msgid)
39 /* Write the legend of the chart */
41 histogram_write_legend (struct chart *ch, double n, double mean, double stddev)
48 pl_savestate_r (ch->lp);
50 sprintf (buf, "N = %.2f", n);
51 pl_move_r (ch->lp, ch->legend_left, ch->data_bottom);
52 pl_alabel_r (ch->lp, 0, 'b', buf);
54 sprintf (buf, "Mean = %.1f", mean);
55 pl_fmove_r (ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5);
56 pl_alabel_r (ch->lp, 0, 'b', buf);
58 sprintf (buf, "Std. Dev = %.2f", stddev);
59 pl_fmove_r (ch->lp, ch->legend_left, ch->data_bottom + ch->font_size * 1.5 * 2);
60 pl_alabel_r (ch->lp, 0, 'b', buf);
62 pl_restorestate_r (ch->lp);
65 static void hist_draw_bar (struct chart *ch, const struct histogram *hist, int bar);
69 hist_draw_bar (struct chart *ch, const struct histogram *hist, int bar)
79 const size_t bins = gsl_histogram_bins (hist->gsl_hist);
80 const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
81 const double width = (ch->data_right - ch->data_left) / (double) bins ;
83 assert ( 0 == gsl_histogram_get_range (hist->gsl_hist, bar, &lower, &upper));
85 assert ( upper >= lower);
87 height = gsl_histogram_get (hist->gsl_hist, bar) *
88 (ch->data_top - ch->data_bottom) / gsl_histogram_max_val (hist->gsl_hist);
90 pl_savestate_r (ch->lp);
91 pl_move_r (ch->lp,ch->data_left, ch->data_bottom);
92 pl_fillcolorname_r (ch->lp, ch->fill_colour);
93 pl_filltype_r (ch->lp,1);
98 x_pos + width, height);
100 pl_restorestate_r (ch->lp);
104 snprintf (buf,5,"%g", (upper + lower) / 2.0);
105 draw_tick (ch, TICK_ABSCISSA,
106 x_pos + width / 2.0, buf);
114 histogram_plot (const struct histogram *hist,
116 const struct moments1 *m)
120 moments1_calculate (m, &n, &mean, &var, NULL, NULL);
122 histogram_plot_n (hist, label, n, mean, sqrt(var), m);
126 /* This function is deprecated. Don't use it in new code */
128 histogram_plot_n (const struct histogram *hist,
130 double n, double mean, double stddev,
136 struct chart *ch = chart_create ();
138 chart_write_title (ch, _("HISTOGRAM"));
140 chart_write_ylabel (ch, _("Frequency"));
141 chart_write_xlabel (ch, label);
143 if ( ! hist ) /* If this happens, probably all values are SYSMIS */
150 bins = gsl_histogram_bins (hist->gsl_hist);
153 chart_write_yscale (ch, 0, gsl_histogram_max_val (hist->gsl_hist), 5);
155 for ( i = 0 ; i < bins ; ++i )
156 hist_draw_bar (ch, hist, i);
158 histogram_write_legend (ch, n, mean, stddev);
162 /* Draw the normal curve */
165 double x_min, x_max, not_used ;
166 double abscissa_scale ;
167 double ordinate_scale ;
170 gsl_histogram_get_range (hist->gsl_hist, 0, &x_min, ¬_used);
171 range = not_used - x_min;
172 gsl_histogram_get_range (hist->gsl_hist, bins - 1, ¬_used, &x_max);
174 abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
175 ordinate_scale = (ch->data_top - ch->data_bottom) /
176 gsl_histogram_max_val (hist->gsl_hist) ;
178 pl_move_r (ch->lp, ch->data_left, ch->data_bottom);
179 for ( d = ch->data_left;
180 d <= ch->data_right ;
181 d += (ch->data_right - ch->data_left) / 100.0)
183 const double x = (d - ch->data_left) / abscissa_scale + x_min ;
184 const double y = n * range *
185 gsl_ran_gaussian_pdf (x - mean, stddev);
187 pl_fcont_r (ch->lp, d, ch->data_bottom + y * ordinate_scale);
190 pl_endpath_r (ch->lp);