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);
99 pl_fillcolorname_r (lp, geom->fill_colour);
100 pl_filltype_r (lp,1);
105 x_pos + width, height);
107 pl_restorestate_r (lp);
109 draw_tick (lp, geom, TICK_ABSCISSA,
110 x_pos + width / 2.0, "%g", (upper + lower) / 2.0);
113 struct histogram_chart
116 gsl_histogram *gsl_hist;
124 /* Plots a histogram of the data in HIST with the given LABEL.
125 Labels the histogram with each of N, MEAN, and STDDEV that is
126 not SYSMIS. If all three are not SYSMIS and SHOW_NORMAL is
127 true, also draws a normal curve on the histogram. */
129 histogram_chart_create (const struct histogram *hist, const char *label,
130 double n, double mean, double stddev,
133 struct histogram_chart *h;
135 h = xmalloc (sizeof *h);
136 chart_init (&h->chart, &histogram_chart_class);
137 h->gsl_hist = hist->gsl_hist ? gsl_histogram_clone (hist->gsl_hist) : NULL;
138 h->label = xstrdup (label);
142 h->show_normal = show_normal;
147 histogram_chart_draw (const struct chart *chart, plPlotter *lp)
149 struct histogram_chart *h = (struct histogram_chart *) chart;
150 struct chart_geometry geom;
154 chart_geometry_init (lp, &geom);
156 chart_write_title (lp, &geom, _("HISTOGRAM"));
158 chart_write_ylabel (lp, &geom, _("Frequency"));
159 chart_write_xlabel (lp, &geom, h->label);
161 if (h->gsl_hist == NULL)
163 /* Probably all values are SYSMIS. */
167 bins = gsl_histogram_bins (h->gsl_hist);
169 chart_write_yscale (lp, &geom, 0, gsl_histogram_max_val (h->gsl_hist), 5);
171 for (i = 0; i < bins; i++)
172 hist_draw_bar (lp, &geom, h->gsl_hist, i);
174 histogram_write_legend (lp, &geom, h->n, h->mean, h->stddev);
177 && h->n != SYSMIS && h->mean != SYSMIS && h->stddev != SYSMIS)
179 /* Draw the normal curve */
181 double x_min, x_max, not_used;
182 double abscissa_scale;
183 double ordinate_scale;
186 gsl_histogram_get_range (h->gsl_hist, 0, &x_min, ¬_used);
187 range = not_used - x_min;
188 gsl_histogram_get_range (h->gsl_hist, bins - 1, ¬_used, &x_max);
190 abscissa_scale = (geom.data_right - geom.data_left) / (x_max - x_min);
191 ordinate_scale = (geom.data_top - geom.data_bottom) /
192 gsl_histogram_max_val (h->gsl_hist);
194 pl_move_r (lp, geom.data_left, geom.data_bottom);
195 for (d = geom.data_left;
196 d <= geom.data_right;
197 d += (geom.data_right - geom.data_left) / 100.0)
199 const double x = (d - geom.data_left) / abscissa_scale + x_min;
200 const double y = h->n * range *
201 gsl_ran_gaussian_pdf (x - h->mean, h->stddev);
203 pl_fcont_r (lp, d, geom.data_bottom + y * ordinate_scale);
209 chart_geometry_free (lp);
214 histogram_chart_destroy (struct chart *chart)
216 struct histogram_chart *h = (struct histogram_chart *) chart;
217 if (h->gsl_hist != NULL)
218 gsl_histogram_free (h->gsl_hist);
223 static const struct chart_class histogram_chart_class =
225 histogram_chart_draw,
226 histogram_chart_destroy