diff --git a/src/language/stats/descriptives.c b/src/language/stats/descriptives
index 8efb5f1..596b977 100644
--- a/src/language/stats/descriptives.c
+++ b/src/language/stats/descriptives.c
@@ -85,7 +85,7 @@ enum dsc_statistic
   {
     DSC_MEAN = 0, DSC_SEMEAN, DSC_STDDEV, DSC_VARIANCE, DSC_KURTOSIS,
     DSC_SEKURT, DSC_SKEWNESS, DSC_SESKEW, DSC_RANGE, DSC_MIN,
-    DSC_MAX, DSC_SUM, DSC_N_STATS,
+    DSC_MAX, DSC_SUM, DSC_HARMONIC, DSC_N_STATS,
 
     /* Only valid as sort criteria. */
     DSC_NAME = -2,              /* Sort by name. */
@@ -115,6 +115,7 @@ static const struct dsc_statistic_info dsc_info[DSC_N_STATS]
     {"MINIMUM", N_("Minimum"), MOMENT_NONE},
     {"MAXIMUM", N_("Maximum"), MOMENT_NONE},
     {"SUM", N_("Sum"), MOMENT_MEAN},
+    {"HARMONIC", N_("Harmonic Mean"), MOMENT_MEAN},
   };
 
 /* Statistics calculated by default if none are explicitly
@@ -884,9 +885,13 @@ calc_descriptives (struct dsc_proc *dsc, struct casereader 
       dv->valid = W = dsc->valid - dv->missing;
 
       if (dv->moments != NULL)
-        moments_calculate (dv->moments, NULL,
+        {
+          moments_calculate (dv->moments, NULL,
                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
+          moments_calc_additional (dv->moments, NULL, &dv->stats[DSC_HARMONIC])
+        }
+      
       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
diff --git a/src/language/stats/frequencies.q b/src/language/stats/frequencies.q
index b5ad8c4..d42031b 100644
--- a/src/language/stats/frequencies.q
+++ b/src/language/stats/frequencies.q
@@ -81,7 +81,7 @@
      +percentiles = double list;
      +statistics[st_]=mean,semean,median,mode,stddev,variance,
            kurtosis,skewness,range,minimum,maximum,sum,
-           default,seskewness,sekurtosis,all,none.
+           default,seskewness,sekurtosis,harmonic,all,none.
 */
 /* (declarations) */
 /* (functions) */
@@ -91,7 +91,7 @@ enum
   {
     FRQ_MEAN, FRQ_SEMEAN, FRQ_MEDIAN, FRQ_MODE, FRQ_STDDEV, FRQ_VARIANCE,
     FRQ_KURT, FRQ_SEKURT, FRQ_SKEW, FRQ_SESKEW, FRQ_RANGE, FRQ_MIN, FRQ_MAX,
-    FRQ_SUM, FRQ_N_STATS
+    FRQ_SUM, FRQ_HARMONIC, FRQ_N_STATS
   };
 
 /* Description of a statistic. */
@@ -118,6 +118,7 @@ static const struct frq_info st_name[FRQ_N_STATS + 1] =
   {FRQ_ST_MINIMUM, N_("Minimum")},
   {FRQ_ST_MAXIMUM, N_("Maximum")},
   {FRQ_ST_SUM, N_("Sum")},
+  {FRQ_ST_HARMONIC, N_("Harmonic Mean")},
   {-1, 0},
 };
 
@@ -139,7 +140,9 @@ struct percentile
         | BI(FRQ_MEAN) | BI(FRQ_SEMEAN) | BI(FRQ_STDDEV)       \
         | BI(FRQ_VARIANCE) | BI(FRQ_KURT) | BI(FRQ_SEKURT)     \
         | BI(FRQ_SKEW) | BI(FRQ_SESKEW) | BI(FRQ_RANGE)        \
-        | BI(FRQ_RANGE) | BI(FRQ_MODE) | BI(FRQ_MEDIAN))
+        | BI(FRQ_RANGE) | BI(FRQ_MODE) | BI(FRQ_MEDIAN)  \
+        | BI(FRQ_HARMONIC))
+
 
 struct frq_chart
   {
@@ -987,6 +990,7 @@ calc_stats (const struct var_freqs *vf, double d[FRQ_N_STATS
     moments_pass_two (m, f->value.f, f->count);
   moments_calculate (m, NULL, &d[FRQ_MEAN], &d[FRQ_VARIANCE],
                      &d[FRQ_SKEW], &d[FRQ_KURT]);
+  moments_calc_additional (m, NULL, &d[FRQ_HARMONIC]);
   moments_destroy (m);
 
   /* Formulas below are taken from _SPSS Statistical Algorithms_. */
diff --git a/src/math/moments.c b/src/math/moments.c
index 40180f7..f44d0b3 100644
--- a/src/math/moments.c
+++ b/src/math/moments.c
@@ -84,8 +84,11 @@ struct moments
     /* Pass one. */
     double w1;                  /* Total weight for pass 1, so far. */
     double sum;                 /* Sum of values so far. */
+    double sumreci;             /* Sum of reciprocals of values */
+    
     double mean;                /* Mean = sum / w1. */
-
+    double hmean;               /* Harmonic mean = w1 / sumreci, when all value
+    
     /* Pass two. */
     double w2;                  /* Total weight for pass 2, so far. */
     double d1;                  /* Sum of deviations from the mean. */
@@ -113,7 +116,7 @@ moments_clear (struct moments *m)
 {
   m->pass = 1;
   m->w1 = m->w2 = 0.;
-  m->sum = 0.;
+  m->sum = m->sumreci = 0.;
 }
 
 /* Creates and returns a data structure for calculating moment
@@ -146,6 +149,8 @@ moments_pass_one (struct moments *m, double value, double we
     {
       m->sum += value * weight;
       m->w1 += weight;
+      /* Sum of reciprocals only needed for variables with positive values */
+      m->sumreci = (value <= 0. || m->sumreci == SYSMIS) ? SYSMIS : m->sumreci 
     }
 }
 
@@ -160,6 +165,7 @@ moments_pass_two (struct moments *m, double value, double we
     {
       m->pass = 2;
       m->mean = m->w1 != 0. ? m->sum / m->w1 : 0.;
+      m->hmean = ( m->sumreci == SYSMIS ) ? SYSMIS : m->w1 / m->sumreci;
       m->d1 = m->d2 = m->d3 = m->d4 = 0.;
     }
 
@@ -239,6 +245,34 @@ moments_calculate (const struct moments *m,
     }
 }
 
+void
+moments_calc_additional (const struct moments *m, double *weight, double *hmean
+{
+  assert (m != NULL);
+  if (hmean != NULL)
+    *hmean = SYSMIS;
+  
+  if (weight != NULL)
+    *weight = m->w1;
+  
+  if (m->pass == 1)
+  {
+    if (hmean != NULL && m->w1 > 0.)
+      *hmean = ( m->sumreci == SYSMIS ) ? SYSMIS : m->w1 / m->sumreci;
+  }
+  else
+  {
+    assert (m->pass == 2);
+    assert (m->w1 == m->w2);
+    
+    if (m->w2 > 0.)
+    {
+      if (hmean != NULL)
+        *hmean = m->hmean;
+    }
+  }
+}
+
 /* Destroys a set of moments. */
 void
 moments_destroy (struct moments *m)
diff --git a/src/math/moments.h b/src/math/moments.h
index 8bf9165..c8d5863 100644
--- a/src/math/moments.h
+++ b/src/math/moments.h
@@ -42,6 +42,7 @@ void moments_calculate (const struct moments *,
                         double *weight,
                         double *mean, double *variance,
                         double *skewness, double *kurtosis);
+void moments_calc_additional (const struct moments *m, double *weight, double *
 void moments_destroy (struct moments *);
 
 /* Convenience functions for two-pass moments. */
diff --git a/src/ui/gui/psppire-dialog-action-descriptives.c b/src/ui/gui/psppir
index b83cea8..c35e5b9 100644
--- a/src/ui/gui/psppire-dialog-action-descriptives.c
+++ b/src/ui/gui/psppire-dialog-action-descriptives.c
@@ -45,7 +45,8 @@ G_DEFINE_TYPE (PsppireDialogActionDescriptives, psppire_dialog
   DS (SEMEAN, N_("Standard error"))             \
   DS (VARIANCE, N_("Variance"))                 \
   DS (KURTOSIS, N_("Kurtosis"))                 \
-  DS (SKEWNESS, N_("Skewness"))
+  DS (SKEWNESS, N_("Skewness"))                 \
+  DS (HARMONIC, N_("Harmonic mean"))
 
 enum
   {
diff --git a/src/ui/gui/psppire-dialog-action-frequencies.c b/src/ui/gui/psppire
index 6f80490..a786c22 100644
--- a/src/ui/gui/psppire-dialog-action-frequencies.c
+++ b/src/ui/gui/psppire-dialog-action-frequencies.c
@@ -50,9 +50,8 @@ This program is free software: you can redistribute it and/or 
   FS (KURTOSIS, N_("Kurtosis"))                                \
   FS (SEKURT, N_("Standard error of the kurtosis"))    \
   FS (MEDIAN, N_("Median"))                            \
-  FS (SUM, N_("Sum"))
-
-
+  FS (SUM, N_("Sum"))             \
+  FS (HARMONIC, N_("Harmonic mean"))
 
 enum
 {
