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);
115 draw_tick (ch, TICK_ABSCISSA,
116 x_pos + width / 2.0, "%g", (upper + lower) / 2.0);
122 /* Plots a histogram of the data in HIST with the given LABEL.
123 Labels the histogram with each of N, MEAN, and STDDEV that is
124 not SYSMIS. If all three are not SYSMIS and SHOW_NORMAL is
125 true, also draws a normal curve on the histogram. */
127 histogram_plot (const struct histogram *hist,
129 double n, double mean, double stddev,
135 struct chart *ch = chart_create ();
137 chart_write_title (ch, _("HISTOGRAM"));
139 chart_write_ylabel (ch, _("Frequency"));
140 chart_write_xlabel (ch, label);
142 if ( ! hist ) /* If this happens, probably all values are SYSMIS */
149 bins = gsl_histogram_bins (hist->gsl_hist);
152 chart_write_yscale (ch, 0, gsl_histogram_max_val (hist->gsl_hist), 5);
154 for ( i = 0 ; i < bins ; ++i )
155 hist_draw_bar (ch, hist, i);
157 histogram_write_legend (ch, n, mean, stddev);
159 if (show_normal && n != SYSMIS && mean != SYSMIS && stddev != SYSMIS)
161 /* Draw the normal curve */
164 double x_min, x_max, not_used ;
165 double abscissa_scale ;
166 double ordinate_scale ;
169 gsl_histogram_get_range (hist->gsl_hist, 0, &x_min, ¬_used);
170 range = not_used - x_min;
171 gsl_histogram_get_range (hist->gsl_hist, bins - 1, ¬_used, &x_max);
173 abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
174 ordinate_scale = (ch->data_top - ch->data_bottom) /
175 gsl_histogram_max_val (hist->gsl_hist) ;
177 pl_move_r (ch->lp, ch->data_left, ch->data_bottom);
178 for ( d = ch->data_left;
179 d <= ch->data_right ;
180 d += (ch->data_right - ch->data_left) / 100.0)
182 const double x = (d - ch->data_left) / abscissa_scale + x_min ;
183 const double y = n * range *
184 gsl_ran_gaussian_pdf (x - mean, stddev);
186 pl_fcont_r (ch->lp, d, ch->data_bottom + y * ordinate_scale);
189 pl_endpath_r (ch->lp);