+
+
+
+/* Plot the normal and detrended normal plots for m
+ Label the plots with factorname */
+void
+np_plot(const struct metrics *m, const char *factorname)
+{
+ int i;
+ double yfirst=0, ylast=0;
+
+ /* Normal Plot */
+ struct chart np_chart;
+
+ /* Detrended Normal Plot */
+ struct chart dnp_chart;
+
+ const struct weighted_value *wv = m->wv;
+ const int n_data = hsh_count(m->ordered_data) ;
+
+ /* The slope and intercept of the ideal normal probability line */
+ const double slope = 1.0 / m->stddev;
+ const double intercept = - m->mean / m->stddev;
+
+ chart_initialise(&np_chart);
+ chart_write_title(&np_chart, _("Normal Q-Q Plot of %s"), factorname);
+ chart_write_xlabel(&np_chart, _("Observed Value"));
+ chart_write_ylabel(&np_chart, _("Expected Normal"));
+
+ chart_initialise(&dnp_chart);
+ chart_write_title(&dnp_chart, _("Detrended Normal Q-Q Plot of %s"),
+ factorname);
+ chart_write_xlabel(&dnp_chart, _("Observed Value"));
+ chart_write_ylabel(&dnp_chart, _("Dev from Normal"));
+
+ yfirst = gsl_cdf_ugaussian_Pinv (wv[0].rank / ( m->n + 1));
+ ylast = gsl_cdf_ugaussian_Pinv (wv[n_data-1].rank / ( m->n + 1));
+
+ {
+ /* Need to make sure that both the scatter plot and the ideal fit into the
+ plot */
+ double x_lower = min(m->min, (yfirst - intercept) / slope) ;
+ double x_upper = max(m->max, (ylast - intercept) / slope) ;
+ double slack = (x_upper - x_lower) * 0.05 ;
+
+ chart_write_xscale(&np_chart, x_lower - slack, x_upper + slack,
+ chart_rounded_tick((m->max - m->min) / 5.0));
+
+
+ chart_write_xscale(&dnp_chart, m->min, m->max,
+ chart_rounded_tick((m->max - m->min) / 5.0));
+
+ }
+
+ chart_write_yscale(&np_chart, yfirst, ylast,
+ chart_rounded_tick((ylast - yfirst)/5.0) );
+
+ {
+ /* We have to cache the detrended data, beacause we need to
+ find its limits before we can plot it */
+ double *d_data;
+ d_data = xmalloc (n_data * sizeof(double));
+ double d_max = -DBL_MAX;
+ double d_min = DBL_MAX;
+ for ( i = 0 ; i < n_data; ++i )
+ {
+ const double ns = gsl_cdf_ugaussian_Pinv (wv[i].rank / ( m->n + 1));
+
+ chart_datum(&np_chart, 0, wv[i].v.f, ns);
+
+ d_data[i] = (wv[i].v.f - m->mean) / m->stddev - ns;
+
+ if ( d_data[i] < d_min ) d_min = d_data[i];
+ if ( d_data[i] > d_max ) d_max = d_data[i];
+ }
+
+ chart_write_yscale(&dnp_chart, d_min, d_max,
+ chart_rounded_tick((d_max - d_min) / 5.0));
+
+ for ( i = 0 ; i < n_data; ++i )
+ chart_datum(&dnp_chart, 0, wv[i].v.f, d_data[i]);
+
+ free(d_data);
+ }
+
+ chart_line(&np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
+ chart_line(&dnp_chart, 0, 0, m->min, m->max , CHART_DIM_X);
+
+ chart_finalise(&np_chart);
+ chart_finalise(&dnp_chart);
+
+}