Fixed a bug which crept into the npplot function
[pspp] / src / examine.q
index 697176f3269f23cb6af8c2e9b67efc023896a823..4444f9652a8ca874b06e08a76809659c2f1af7bb 100644 (file)
@@ -41,6 +41,8 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 #include "hash.h"
 #include "casefile.h"
 #include "factor_stats.h"
+#include "moments.h"
+
 /* (headers) */
 #include "chart.h"
 
@@ -80,7 +82,7 @@ struct factor
   /* Hash table of factor stats indexed by 2 values */
   struct hsh_table *fstats;
 
-  /* The hash table after it's been crunched */
+  /* The hash table after it has been crunched */
   struct factor_statistics **fs;
 
   struct factor *next;
@@ -120,8 +122,6 @@ print_factors(void)
          */
 
                 
-         printf("Sum is %g; ",(*fs)->m[0].sum);
-         printf("N is %g; ",(*fs)->m[0].n);
          printf("Mean is %g\n",(*fs)->m[0].mean);
 
          fs++ ;
@@ -193,6 +193,9 @@ cmd_examine(void)
 
   multipass_procedure_with_splits (run_examine, &cmd);
 
+  if ( totals ) 
+    free(totals);
+
   return CMD_SUCCESS;
 };
 
@@ -205,7 +208,7 @@ output_examine(void)
   struct factor *fctr;
 
   /* Show totals if appropriate */
-  if ( ! cmd.sbc_nototal )
+  if ( ! cmd.sbc_nototal || factors == 0 )
     {
       show_summary(dependent_vars, n_dependent_vars, 0);
 
@@ -638,12 +641,11 @@ run_examine(const struct casefile *cf, void *cmd_ )
       fctr = fctr->next;
     }
 
-  /* 
-  print_factors();
-  */
-
   output_examine();
 
+  for ( v = 0 ; v < n_dependent_vars ; ++v ) 
+    hsh_destroy(totals[v].ordered_data);
+
 }
 
 
@@ -1021,7 +1023,6 @@ populate_extremes(struct tab_table *t,
   int extremity;
   int idx=0;
 
-  const int n_data = hsh_count(m->ordered_data);
 
   tab_text(t, col, row,
           TAB_RIGHT | TAT_TITLE ,
@@ -1053,10 +1054,10 @@ populate_extremes(struct tab_table *t,
 
 
   /* Lowest */
-  for (idx = 0, extremity = 0; extremity < n && idx < n_data ; ++idx ) 
+  for (idx = 0, extremity = 0; extremity < n && idx < m->n_data ; ++idx ) 
     {
       int j;
-      const struct weighted_value *wv = &m->wv[idx];
+      const struct weighted_value *wv = m->wvp[idx];
       struct case_node *cn = wv->case_nos;
 
       
@@ -1083,10 +1084,10 @@ populate_extremes(struct tab_table *t,
 
 
   /* Highest */
-  for (idx = n_data - 1, extremity = 0; extremity < n && idx >= 0; --idx ) 
+  for (idx = m->n_data - 1, extremity = 0; extremity < n && idx >= 0; --idx ) 
     {
       int j;
-      const struct weighted_value *wv = &m->wv[idx];
+      const struct weighted_value *wv = m->wvp[idx];
       struct case_node *cn = wv->case_nos;
 
       for (j = 0 ; j < wv->w ; ++j  )
@@ -1288,7 +1289,7 @@ populate_descriptives(struct tab_table *tbl, int col, int row,
   tab_float (tbl, col + 3,
             row,
             TAB_CENTER,
-            m->stderr,
+            m->se_mean,
             8,3);
   
 
@@ -1306,7 +1307,7 @@ populate_descriptives(struct tab_table *tbl, int col, int row,
   tab_float (tbl, col + 2,
             row + 1,
             TAB_CENTER,
-            m->mean - t * m->stderr
+            m->mean - t * m->se_mean
             8,3);
 
   tab_text (tbl, col + 1,  
@@ -1318,7 +1319,7 @@ populate_descriptives(struct tab_table *tbl, int col, int row,
   tab_float (tbl, col + 2,
             row + 2,
             TAB_CENTER,
-            m->mean + t * m->stderr
+            m->mean + t * m->se_mean
             8,3);
 
   tab_text (tbl, col, 
@@ -1407,10 +1408,41 @@ populate_descriptives(struct tab_table *tbl, int col, int row,
            TAB_LEFT | TAT_TITLE,
            _("Skewness"));
 
+
+  tab_float (tbl, col + 2,
+            row + 11,
+            TAB_CENTER,
+            m->skewness,
+            8,3);
+
+  /* stderr of skewness */
+  tab_float (tbl, col + 3,
+            row + 11,
+            TAB_CENTER,
+            calc_seskew(m->n),
+            8,3);
+
+
   tab_text (tbl, col, 
            row + 12,
            TAB_LEFT | TAT_TITLE,
            _("Kurtosis"));
+
+
+  tab_float (tbl, col + 2,
+            row + 12,
+            TAB_CENTER,
+            m->kurtosis,
+            8,3);
+
+  /* stderr of kurtosis */
+  tab_float (tbl, col + 3,
+            row + 12,
+            TAB_CENTER,
+            calc_sekurt(m->n),
+            8,3);
+
+
 }
 
 
@@ -1431,15 +1463,12 @@ np_plot(const struct metrics *m, const char *factorname)
   /* 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;
 
   /* Cowardly refuse to plot an empty data set */
-  if ( n_data == 0 ) 
+  if ( m->n_data == 0 ) 
     return ; 
 
   chart_initialise(&np_chart);
@@ -1453,8 +1482,8 @@ np_plot(const struct metrics *m, const char *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));
+  yfirst = gsl_cdf_ugaussian_Pinv (m->wvp[0]->rank / ( m->n + 1));
+  ylast =  gsl_cdf_ugaussian_Pinv (m->wvp[m->n_data-1]->rank / ( m->n + 1));
 
   {
     /* Need to make sure that both the scatter plot and the ideal fit into the
@@ -1479,16 +1508,16 @@ np_plot(const struct metrics *m, const char *factorname)
   /* 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));
+  d_data = xmalloc (m->n_data * sizeof(double));
   double d_max = -DBL_MAX;
   double d_min = DBL_MAX;
-  for ( i = 0 ; i < n_data; ++i ) 
+  for ( i = 0 ; i < m->n_data; ++i ) 
     {
-      const double ns = gsl_cdf_ugaussian_Pinv (wv[i].rank / ( m->n + 1));
+      const double ns = gsl_cdf_ugaussian_Pinv (m->wvp[i]->rank / ( m->n + 1));
 
-      chart_datum(&np_chart, 0, wv[i].v.f, ns);
+      chart_datum(&np_chart, 0, m->wvp[i]->v.f, ns);
 
-      d_data[i] = (wv[i].v.f - m->mean) / m->stddev  - ns;
+      d_data[i] = (m->wvp[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];
@@ -1497,8 +1526,8 @@ np_plot(const struct metrics *m, const char *factorname)
   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]);
+  for ( i = 0 ; i < m->n_data; ++i ) 
+      chart_datum(&dnp_chart, 0, m->wvp[i]->v.f, d_data[i]);
 
   free(d_data);
   }