+\f
+/* One-pass moments. */
+
+/* A set of one-pass moments. */
+struct moments1
+ {
+ enum moment max_moment; /* Highest-order moment we're computing. */
+ double w; /* Total weight so far. */
+ double d1; /* Sum of deviations from the mean. */
+ double d2; /* Sum of squared deviations from the mean. */
+ double d3; /* Sum of cubed deviations from the mean. */
+ double d4; /* Sum of (deviations from the mean)**4. */
+ };
+
+/* Initializes one-pass moments M for calculating moment
+ MAX_MOMENT and lower moments. */
+static void
+init_moments1 (struct moments1 *m, enum moment max_moment)
+{
+ assert (m != NULL);
+ assert (max_moment == MOMENT_MEAN || max_moment == MOMENT_VARIANCE
+ || max_moment == MOMENT_SKEWNESS || max_moment == MOMENT_KURTOSIS);
+ m->max_moment = max_moment;
+ moments1_clear (m);
+}
+
+/* Clears out a set of one-pass moments so that it can be reused
+ for a new set of values. The moments to be calculated are not
+ changed. */
+void
+moments1_clear (struct moments1 *m)
+{
+ m->w = 0.;
+ m->d1 = m->d2 = m->d3 = m->d4 = 0.;
+}
+
+/* Creates and returns a data structure for calculating moment
+ MAX_MOMENT and lower moments on a data series in a single
+ pass. The user should call moments1_add() for each value in
+ the series. The user may call moments1_calculate() to obtain
+ the current moments at any time. Call moments1_destroy() when
+ the moments are no longer needed.
+
+ One-pass moments should only be used when two passes over the
+ data are impractical. */
+struct moments1 *
+moments1_create (enum moment max_moment)
+{
+ struct moments1 *m = xmalloc (sizeof *m);
+ init_moments1 (m, max_moment);
+ return m;
+}
+
+/* Adds VALUE with the given WEIGHT to the calculation of
+ one-pass moments. */
+void
+moments1_add (struct moments1 *m, double value, double weight)
+{
+ assert (m != NULL);
+
+ if (value != SYSMIS && weight > 0.)
+ {
+ double prev_w, v1;
+
+ prev_w = m->w;
+ m->w += weight;
+ v1 = (weight / m->w) * (value - m->d1);
+ m->d1 += v1;
+
+ if (m->max_moment >= MOMENT_VARIANCE)
+ {
+ double v2 = v1 * v1;
+ double w_prev_w = m->w * prev_w;
+ double prev_m2 = m->d2;
+
+ m->d2 += w_prev_w / weight * v2;
+ if (m->max_moment >= MOMENT_SKEWNESS)
+ {
+ double w2 = weight * weight;
+ double v3 = v2 * v1;
+ double prev_m3 = m->d3;
+
+ m->d3 += (-3. * v1 * prev_m2
+ + w_prev_w / w2 * (m->w - 2. * weight) * v3);
+ if (m->max_moment >= MOMENT_KURTOSIS)
+ {
+ double w3 = w2 * weight;
+ double v4 = v2 * v2;
+
+ m->d4 += (-4. * v1 * prev_m3
+ + 6. * v2 * prev_m2
+ + ((pow2 (m->w) - 3. * weight * prev_w)
+ * v4 * w_prev_w / w3));
+ }
+ }
+ }
+ }
+}
+
+/* Calculates one-pass moments based on the input data. Stores
+ the total weight in *WEIGHT, the mean in *MEAN, the variance
+ in *VARIANCE, the skewness in *SKEWNESS, and the kurtosis in
+ *KURTOSIS. Any of these result parameters may be null
+ pointers, in which case the values are not calculated. If any
+ result cannot be calculated, either because they are undefined
+ based on the input data or because their moments are higher
+ than the maximum requested on moments_create(), then SYSMIS is
+ stored into that result. */
+void
+moments1_calculate (const struct moments1 *m,
+ double *weight,
+ double *mean, double *variance,
+ double *skewness, double *kurtosis)
+{
+ assert (m != NULL);