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/>. */
22 #include <gsl/gsl_histogram.h>
23 #include <gsl/gsl_randist.h>
26 #include <output/charts/plot-hist.h>
27 #include <output/charts/plot-chart.h>
28 #include <output/chart-provider.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 static const struct chart_class histogram_chart_class;
41 /* Write the legend of the chart */
43 histogram_write_legend (plPlotter *lp, const struct chart_geometry *geom,
44 double n, double mean, double stddev)
46 double y = geom->data_bottom;
51 char *buf = xasprintf ("N = %.2f", n);
52 pl_move_r (lp, geom->legend_left, y);
53 pl_alabel_r (lp, 0, 'b', buf);
54 y += geom->font_size * 1.5;
60 char *buf = xasprintf ("Mean = %.1f", mean);
61 pl_fmove_r (lp,geom->legend_left, y);
62 pl_alabel_r (lp, 0, 'b', buf);
63 y += geom->font_size * 1.5;
69 char *buf = xasprintf ("Std. Dev = %.2f", stddev);
70 pl_fmove_r (lp, geom->legend_left, y);
71 pl_alabel_r (lp, 0, 'b', buf);
75 pl_restorestate_r (lp);
79 hist_draw_bar (plPlotter *lp, const struct chart_geometry *geom,
80 const gsl_histogram *h, int bar)
86 const size_t bins = gsl_histogram_bins (h);
87 const double x_pos = (geom->data_right - geom->data_left) * bar / (double) bins ;
88 const double width = (geom->data_right - geom->data_left) / (double) bins ;
90 assert ( 0 == gsl_histogram_get_range (h, bar, &lower, &upper));
92 assert ( upper >= lower);
94 height = gsl_histogram_get (h, bar) *
95 (geom->data_top - geom->data_bottom) / gsl_histogram_max_val (h);
98 pl_move_r (lp,geom->data_left, geom->data_bottom);
100 geom->fill_colour.red * 257,
101 geom->fill_colour.green * 257,
102 geom->fill_colour.blue * 257);
103 pl_filltype_r (lp,1);
108 x_pos + width, height);
110 pl_restorestate_r (lp);
112 draw_tick (lp, geom, TICK_ABSCISSA,
113 x_pos + width / 2.0, "%g", (upper + lower) / 2.0);
116 struct histogram_chart
119 gsl_histogram *gsl_hist;
127 /* Plots a histogram of the data in HIST with the given LABEL.
128 Labels the histogram with each of N, MEAN, and STDDEV that is
129 not SYSMIS. If all three are not SYSMIS and SHOW_NORMAL is
130 true, also draws a normal curve on the histogram. */
132 histogram_chart_create (const struct histogram *hist, const char *label,
133 double n, double mean, double stddev,
136 struct histogram_chart *h;
138 h = xmalloc (sizeof *h);
139 chart_init (&h->chart, &histogram_chart_class);
140 h->gsl_hist = hist->gsl_hist ? gsl_histogram_clone (hist->gsl_hist) : NULL;
141 h->label = xstrdup (label);
145 h->show_normal = show_normal;
150 histogram_chart_draw (const struct chart *chart, plPlotter *lp)
152 struct histogram_chart *h = (struct histogram_chart *) chart;
153 struct chart_geometry geom;
157 chart_geometry_init (lp, &geom);
159 chart_write_title (lp, &geom, _("HISTOGRAM"));
161 chart_write_ylabel (lp, &geom, _("Frequency"));
162 chart_write_xlabel (lp, &geom, h->label);
164 if (h->gsl_hist == NULL)
166 /* Probably all values are SYSMIS. */
170 bins = gsl_histogram_bins (h->gsl_hist);
172 chart_write_yscale (lp, &geom, 0, gsl_histogram_max_val (h->gsl_hist), 5);
174 for (i = 0; i < bins; i++)
175 hist_draw_bar (lp, &geom, h->gsl_hist, i);
177 histogram_write_legend (lp, &geom, h->n, h->mean, h->stddev);
180 && h->n != SYSMIS && h->mean != SYSMIS && h->stddev != SYSMIS)
182 /* Draw the normal curve */
184 double x_min, x_max, not_used;
185 double abscissa_scale;
186 double ordinate_scale;
189 gsl_histogram_get_range (h->gsl_hist, 0, &x_min, ¬_used);
190 range = not_used - x_min;
191 gsl_histogram_get_range (h->gsl_hist, bins - 1, ¬_used, &x_max);
193 abscissa_scale = (geom.data_right - geom.data_left) / (x_max - x_min);
194 ordinate_scale = (geom.data_top - geom.data_bottom) /
195 gsl_histogram_max_val (h->gsl_hist);
197 pl_move_r (lp, geom.data_left, geom.data_bottom);
198 for (d = geom.data_left;
199 d <= geom.data_right;
200 d += (geom.data_right - geom.data_left) / 100.0)
202 const double x = (d - geom.data_left) / abscissa_scale + x_min;
203 const double y = h->n * range *
204 gsl_ran_gaussian_pdf (x - h->mean, h->stddev);
206 pl_fcont_r (lp, d, geom.data_bottom + y * ordinate_scale);
212 chart_geometry_free (lp);
217 histogram_chart_destroy (struct chart *chart)
219 struct histogram_chart *h = (struct histogram_chart *) chart;
220 if (h->gsl_hist != NULL)
221 gsl_histogram_free (h->gsl_hist);
226 static const struct chart_class histogram_chart_class =
228 histogram_chart_draw,
229 histogram_chart_destroy