1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2004, 2009 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)
50 pl_savestate_r (ch->lp);
54 sprintf (buf, "N = %.2f", n);
55 pl_move_r (ch->lp, ch->legend_left, y);
56 pl_alabel_r (ch->lp, 0, 'b', buf);
57 y += ch->font_size * 1.5;
62 sprintf (buf, "Mean = %.1f", mean);
63 pl_fmove_r (ch->lp,ch->legend_left, y);
64 pl_alabel_r (ch->lp, 0, 'b', buf);
65 y += ch->font_size * 1.5;
70 sprintf (buf, "Std. Dev = %.2f", stddev);
71 pl_fmove_r (ch->lp, ch->legend_left, y);
72 pl_alabel_r (ch->lp, 0, 'b', buf);
75 pl_restorestate_r (ch->lp);
78 static void hist_draw_bar (struct chart *ch, const struct histogram *hist, int bar);
82 hist_draw_bar (struct chart *ch, const struct histogram *hist, int bar)
92 const size_t bins = gsl_histogram_bins (hist->gsl_hist);
93 const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
94 const double width = (ch->data_right - ch->data_left) / (double) bins ;
96 assert ( 0 == gsl_histogram_get_range (hist->gsl_hist, bar, &lower, &upper));
98 assert ( upper >= lower);
100 height = gsl_histogram_get (hist->gsl_hist, bar) *
101 (ch->data_top - ch->data_bottom) / gsl_histogram_max_val (hist->gsl_hist);
103 pl_savestate_r (ch->lp);
104 pl_move_r (ch->lp,ch->data_left, ch->data_bottom);
105 pl_fillcolorname_r (ch->lp, ch->fill_colour);
106 pl_filltype_r (ch->lp,1);
109 pl_fboxrel_r (ch->lp,
111 x_pos + width, height);
113 pl_restorestate_r (ch->lp);
117 snprintf (buf,5,"%g", (upper + lower) / 2.0);
118 draw_tick (ch, TICK_ABSCISSA,
119 x_pos + width / 2.0, buf);
126 /* Plots a histogram of the data in HIST with the given LABEL.
127 Labels the histogram with each of N, MEAN, and STDDEV that is
128 not SYSMIS. If all three are not SYSMIS and SHOW_NORMAL is
129 true, also draws a normal curve on the histogram. */
131 histogram_plot (const struct histogram *hist,
133 double n, double mean, double stddev,
139 struct chart *ch = chart_create ();
141 chart_write_title (ch, _("HISTOGRAM"));
143 chart_write_ylabel (ch, _("Frequency"));
144 chart_write_xlabel (ch, label);
146 if ( ! hist ) /* If this happens, probably all values are SYSMIS */
153 bins = gsl_histogram_bins (hist->gsl_hist);
156 chart_write_yscale (ch, 0, gsl_histogram_max_val (hist->gsl_hist), 5);
158 for ( i = 0 ; i < bins ; ++i )
159 hist_draw_bar (ch, hist, i);
161 histogram_write_legend (ch, n, mean, stddev);
163 if (show_normal && n != SYSMIS && mean != SYSMIS && stddev != SYSMIS)
165 /* Draw the normal curve */
168 double x_min, x_max, not_used ;
169 double abscissa_scale ;
170 double ordinate_scale ;
173 gsl_histogram_get_range (hist->gsl_hist, 0, &x_min, ¬_used);
174 range = not_used - x_min;
175 gsl_histogram_get_range (hist->gsl_hist, bins - 1, ¬_used, &x_max);
177 abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
178 ordinate_scale = (ch->data_top - ch->data_bottom) /
179 gsl_histogram_max_val (hist->gsl_hist) ;
181 pl_move_r (ch->lp, ch->data_left, ch->data_bottom);
182 for ( d = ch->data_left;
183 d <= ch->data_right ;
184 d += (ch->data_right - ch->data_left) / 100.0)
186 const double x = (d - ch->data_left) / abscissa_scale + x_min ;
187 const double y = n * range *
188 gsl_ran_gaussian_pdf (x - mean, stddev);
190 pl_fcont_r (ch->lp, d, ch->data_bottom + y * ordinate_scale);
193 pl_endpath_r (ch->lp);