001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math.stat.descriptive.moment;
018    
019    import java.io.Serializable;
020    
021    import org.apache.commons.math.exception.NullArgumentException;
022    import org.apache.commons.math.exception.util.LocalizedFormats;
023    import org.apache.commons.math.stat.descriptive.WeightedEvaluation;
024    import org.apache.commons.math.stat.descriptive.AbstractStorelessUnivariateStatistic;
025    
026    /**
027     * Computes the variance of the available values.  By default, the unbiased
028     * "sample variance" definitional formula is used:
029     * <p>
030     * variance = sum((x_i - mean)^2) / (n - 1) </p>
031     * <p>
032     * where mean is the {@link Mean} and <code>n</code> is the number
033     * of sample observations.</p>
034     * <p>
035     * The definitional formula does not have good numerical properties, so
036     * this implementation does not compute the statistic using the definitional
037     * formula. <ul>
038     * <li> The <code>getResult</code> method computes the variance using
039     * updating formulas based on West's algorithm, as described in
040     * <a href="http://doi.acm.org/10.1145/359146.359152"> Chan, T. F. and
041     * J. G. Lewis 1979, <i>Communications of the ACM</i>,
042     * vol. 22 no. 9, pp. 526-531.</a></li>
043     * <li> The <code>evaluate</code> methods leverage the fact that they have the
044     * full array of values in memory to execute a two-pass algorithm.
045     * Specifically, these methods use the "corrected two-pass algorithm" from
046     * Chan, Golub, Levesque, <i>Algorithms for Computing the Sample Variance</i>,
047     * American Statistician, vol. 37, no. 3 (1983) pp. 242-247.</li></ul>
048     * Note that adding values using <code>increment</code> or
049     * <code>incrementAll</code> and then executing <code>getResult</code> will
050     * sometimes give a different, less accurate, result than executing
051     * <code>evaluate</code> with the full array of values. The former approach
052     * should only be used when the full array of values is not available.</p>
053     * <p>
054     * The "population variance"  ( sum((x_i - mean)^2) / n ) can also
055     * be computed using this statistic.  The <code>isBiasCorrected</code>
056     * property determines whether the "population" or "sample" value is
057     * returned by the <code>evaluate</code> and <code>getResult</code> methods.
058     * To compute population variances, set this property to <code>false.</code>
059     * </p>
060     * <p>
061     * <strong>Note that this implementation is not synchronized.</strong> If
062     * multiple threads access an instance of this class concurrently, and at least
063     * one of the threads invokes the <code>increment()</code> or
064     * <code>clear()</code> method, it must be synchronized externally.</p>
065     *
066     * @version $Revision: 1006299 $ $Date: 2010-10-10 16:47:17 +0200 (dim. 10 oct. 2010) $
067     */
068    public class Variance extends AbstractStorelessUnivariateStatistic implements Serializable, WeightedEvaluation {
069    
070        /** Serializable version identifier */
071        private static final long serialVersionUID = -9111962718267217978L;
072    
073        /** SecondMoment is used in incremental calculation of Variance*/
074        protected SecondMoment moment = null;
075    
076        /**
077         * Boolean test to determine if this Variance should also increment
078         * the second moment, this evaluates to false when this Variance is
079         * constructed with an external SecondMoment as a parameter.
080         */
081        protected boolean incMoment = true;
082    
083        /**
084         * Determines whether or not bias correction is applied when computing the
085         * value of the statisic.  True means that bias is corrected.  See
086         * {@link Variance} for details on the formula.
087         */
088        private boolean isBiasCorrected = true;
089    
090        /**
091         * Constructs a Variance with default (true) <code>isBiasCorrected</code>
092         * property.
093         */
094        public Variance() {
095            moment = new SecondMoment();
096        }
097    
098        /**
099         * Constructs a Variance based on an external second moment.
100         *
101         * @param m2 the SecondMoment (Third or Fourth moments work
102         * here as well.)
103         */
104        public Variance(final SecondMoment m2) {
105            incMoment = false;
106            this.moment = m2;
107        }
108    
109        /**
110         * Constructs a Variance with the specified <code>isBiasCorrected</code>
111         * property
112         *
113         * @param isBiasCorrected  setting for bias correction - true means
114         * bias will be corrected and is equivalent to using the argumentless
115         * constructor
116         */
117        public Variance(boolean isBiasCorrected) {
118            moment = new SecondMoment();
119            this.isBiasCorrected = isBiasCorrected;
120        }
121    
122        /**
123         * Constructs a Variance with the specified <code>isBiasCorrected</code>
124         * property and the supplied external second moment.
125         *
126         * @param isBiasCorrected  setting for bias correction - true means
127         * bias will be corrected
128         * @param m2 the SecondMoment (Third or Fourth moments work
129         * here as well.)
130         */
131        public Variance(boolean isBiasCorrected, SecondMoment m2) {
132            incMoment = false;
133            this.moment = m2;
134            this.isBiasCorrected = isBiasCorrected;
135        }
136    
137        /**
138         * Copy constructor, creates a new {@code Variance} identical
139         * to the {@code original}
140         *
141         * @param original the {@code Variance} instance to copy
142         */
143        public Variance(Variance original) {
144            copy(original, this);
145        }
146    
147        /**
148         * {@inheritDoc}
149         * <p>If all values are available, it is more accurate to use
150         * {@link #evaluate(double[])} rather than adding values one at a time
151         * using this method and then executing {@link #getResult}, since
152         * <code>evaluate</code> leverages the fact that is has the full
153         * list of values together to execute a two-pass algorithm.
154         * See {@link Variance}.</p>
155         */
156        @Override
157        public void increment(final double d) {
158            if (incMoment) {
159                moment.increment(d);
160            }
161        }
162    
163        /**
164         * {@inheritDoc}
165         */
166        @Override
167        public double getResult() {
168                if (moment.n == 0) {
169                    return Double.NaN;
170                } else if (moment.n == 1) {
171                    return 0d;
172                } else {
173                    if (isBiasCorrected) {
174                        return moment.m2 / (moment.n - 1d);
175                    } else {
176                        return moment.m2 / (moment.n);
177                    }
178                }
179        }
180    
181        /**
182         * {@inheritDoc}
183         */
184        public long getN() {
185            return moment.getN();
186        }
187    
188        /**
189         * {@inheritDoc}
190         */
191        @Override
192        public void clear() {
193            if (incMoment) {
194                moment.clear();
195            }
196        }
197    
198        /**
199         * Returns the variance of the entries in the input array, or
200         * <code>Double.NaN</code> if the array is empty.
201         * <p>
202         * See {@link Variance} for details on the computing algorithm.</p>
203         * <p>
204         * Returns 0 for a single-value (i.e. length = 1) sample.</p>
205         * <p>
206         * Throws <code>IllegalArgumentException</code> if the array is null.</p>
207         * <p>
208         * Does not change the internal state of the statistic.</p>
209         *
210         * @param values the input array
211         * @return the variance of the values or Double.NaN if length = 0
212         * @throws IllegalArgumentException if the array is null
213         */
214        @Override
215        public double evaluate(final double[] values) {
216            if (values == null) {
217                throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
218            }
219            return evaluate(values, 0, values.length);
220        }
221    
222        /**
223         * Returns the variance of the entries in the specified portion of
224         * the input array, or <code>Double.NaN</code> if the designated subarray
225         * is empty.
226         * <p>
227         * See {@link Variance} for details on the computing algorithm.</p>
228         * <p>
229         * Returns 0 for a single-value (i.e. length = 1) sample.</p>
230         * <p>
231         * Does not change the internal state of the statistic.</p>
232         * <p>
233         * Throws <code>IllegalArgumentException</code> if the array is null.</p>
234         *
235         * @param values the input array
236         * @param begin index of the first array element to include
237         * @param length the number of elements to include
238         * @return the variance of the values or Double.NaN if length = 0
239         * @throws IllegalArgumentException if the array is null or the array index
240         *  parameters are not valid
241         */
242        @Override
243        public double evaluate(final double[] values, final int begin, final int length) {
244    
245            double var = Double.NaN;
246    
247            if (test(values, begin, length)) {
248                clear();
249                if (length == 1) {
250                    var = 0.0;
251                } else if (length > 1) {
252                    Mean mean = new Mean();
253                    double m = mean.evaluate(values, begin, length);
254                    var = evaluate(values, m, begin, length);
255                }
256            }
257            return var;
258        }
259    
260        /**
261         * <p>Returns the weighted variance of the entries in the specified portion of
262         * the input array, or <code>Double.NaN</code> if the designated subarray
263         * is empty.</p>
264         * <p>
265         * Uses the formula <pre>
266         *   &Sigma;(weights[i]*(values[i] - weightedMean)<sup>2</sup>)/(&Sigma;(weights[i]) - 1)
267         * </pre>
268         * where weightedMean is the weighted mean</p>
269         * <p>
270         * This formula will not return the same result as the unweighted variance when all
271         * weights are equal, unless all weights are equal to 1. The formula assumes that
272         * weights are to be treated as "expansion values," as will be the case if for example
273         * the weights represent frequency counts. To normalize weights so that the denominator
274         * in the variance computation equals the length of the input vector minus one, use <pre>
275         *   <code>evaluate(values, MathUtils.normalizeArray(weights, values.length)); </code>
276         * </pre>
277         * <p>
278         * Returns 0 for a single-value (i.e. length = 1) sample.</p>
279         * <p>
280         * Throws <code>IllegalArgumentException</code> if any of the following are true:
281         * <ul><li>the values array is null</li>
282         *     <li>the weights array is null</li>
283         *     <li>the weights array does not have the same length as the values array</li>
284         *     <li>the weights array contains one or more infinite values</li>
285         *     <li>the weights array contains one or more NaN values</li>
286         *     <li>the weights array contains negative values</li>
287         *     <li>the start and length arguments do not determine a valid array</li>
288         * </ul></p>
289         * <p>
290         * Does not change the internal state of the statistic.</p>
291         * <p>
292         * Throws <code>IllegalArgumentException</code> if either array is null.</p>
293         *
294         * @param values the input array
295         * @param weights the weights array
296         * @param begin index of the first array element to include
297         * @param length the number of elements to include
298         * @return the weighted variance of the values or Double.NaN if length = 0
299         * @throws IllegalArgumentException if the parameters are not valid
300         * @since 2.1
301         */
302        public double evaluate(final double[] values, final double[] weights,
303                               final int begin, final int length) {
304    
305            double var = Double.NaN;
306    
307            if (test(values, weights,begin, length)) {
308                clear();
309                if (length == 1) {
310                    var = 0.0;
311                } else if (length > 1) {
312                    Mean mean = new Mean();
313                    double m = mean.evaluate(values, weights, begin, length);
314                    var = evaluate(values, weights, m, begin, length);
315                }
316            }
317            return var;
318        }
319    
320        /**
321         * <p>
322         * Returns the weighted variance of the entries in the the input array.</p>
323         * <p>
324         * Uses the formula <pre>
325         *   &Sigma;(weights[i]*(values[i] - weightedMean)<sup>2</sup>)/(&Sigma;(weights[i]) - 1)
326         * </pre>
327         * where weightedMean is the weighted mean</p>
328         * <p>
329         * This formula will not return the same result as the unweighted variance when all
330         * weights are equal, unless all weights are equal to 1. The formula assumes that
331         * weights are to be treated as "expansion values," as will be the case if for example
332         * the weights represent frequency counts. To normalize weights so that the denominator
333         * in the variance computation equals the length of the input vector minus one, use <pre>
334         *   <code>evaluate(values, MathUtils.normalizeArray(weights, values.length)); </code>
335         * </pre>
336         * <p>
337         * Returns 0 for a single-value (i.e. length = 1) sample.</p>
338         * <p>
339         * Throws <code>IllegalArgumentException</code> if any of the following are true:
340         * <ul><li>the values array is null</li>
341         *     <li>the weights array is null</li>
342         *     <li>the weights array does not have the same length as the values array</li>
343         *     <li>the weights array contains one or more infinite values</li>
344         *     <li>the weights array contains one or more NaN values</li>
345         *     <li>the weights array contains negative values</li>
346         * </ul></p>
347         * <p>
348         * Does not change the internal state of the statistic.</p>
349         * <p>
350         * Throws <code>IllegalArgumentException</code> if either array is null.</p>
351         *
352         * @param values the input array
353         * @param weights the weights array
354         * @return the weighted variance of the values
355         * @throws IllegalArgumentException if the parameters are not valid
356         * @since 2.1
357         */
358        public double evaluate(final double[] values, final double[] weights) {
359            return evaluate(values, weights, 0, values.length);
360        }
361    
362        /**
363         * Returns the variance of the entries in the specified portion of
364         * the input array, using the precomputed mean value.  Returns
365         * <code>Double.NaN</code> if the designated subarray is empty.
366         * <p>
367         * See {@link Variance} for details on the computing algorithm.</p>
368         * <p>
369         * The formula used assumes that the supplied mean value is the arithmetic
370         * mean of the sample data, not a known population parameter.  This method
371         * is supplied only to save computation when the mean has already been
372         * computed.</p>
373         * <p>
374         * Returns 0 for a single-value (i.e. length = 1) sample.</p>
375         * <p>
376         * Throws <code>IllegalArgumentException</code> if the array is null.</p>
377         * <p>
378         * Does not change the internal state of the statistic.</p>
379         *
380         * @param values the input array
381         * @param mean the precomputed mean value
382         * @param begin index of the first array element to include
383         * @param length the number of elements to include
384         * @return the variance of the values or Double.NaN if length = 0
385         * @throws IllegalArgumentException if the array is null or the array index
386         *  parameters are not valid
387         */
388        public double evaluate(final double[] values, final double mean,
389                final int begin, final int length) {
390    
391            double var = Double.NaN;
392    
393            if (test(values, begin, length)) {
394                if (length == 1) {
395                    var = 0.0;
396                } else if (length > 1) {
397                    double accum = 0.0;
398                    double dev = 0.0;
399                    double accum2 = 0.0;
400                    for (int i = begin; i < begin + length; i++) {
401                        dev = values[i] - mean;
402                        accum += dev * dev;
403                        accum2 += dev;
404                    }
405                    double len = length;
406                    if (isBiasCorrected) {
407                        var = (accum - (accum2 * accum2 / len)) / (len - 1.0);
408                    } else {
409                        var = (accum - (accum2 * accum2 / len)) / len;
410                    }
411                }
412            }
413            return var;
414        }
415    
416        /**
417         * Returns the variance of the entries in the input array, using the
418         * precomputed mean value.  Returns <code>Double.NaN</code> if the array
419         * is empty.
420         * <p>
421         * See {@link Variance} for details on the computing algorithm.</p>
422         * <p>
423         * If <code>isBiasCorrected</code> is <code>true</code> the formula used
424         * assumes that the supplied mean value is the arithmetic mean of the
425         * sample data, not a known population parameter.  If the mean is a known
426         * population parameter, or if the "population" version of the variance is
427         * desired, set <code>isBiasCorrected</code> to <code>false</code> before
428         * invoking this method.</p>
429         * <p>
430         * Returns 0 for a single-value (i.e. length = 1) sample.</p>
431         * <p>
432         * Throws <code>IllegalArgumentException</code> if the array is null.</p>
433         * <p>
434         * Does not change the internal state of the statistic.</p>
435         *
436         * @param values the input array
437         * @param mean the precomputed mean value
438         * @return the variance of the values or Double.NaN if the array is empty
439         * @throws IllegalArgumentException if the array is null
440         */
441        public double evaluate(final double[] values, final double mean) {
442            return evaluate(values, mean, 0, values.length);
443        }
444    
445        /**
446         * Returns the weighted variance of the entries in the specified portion of
447         * the input array, using the precomputed weighted mean value.  Returns
448         * <code>Double.NaN</code> if the designated subarray is empty.
449         * <p>
450         * Uses the formula <pre>
451         *   &Sigma;(weights[i]*(values[i] - mean)<sup>2</sup>)/(&Sigma;(weights[i]) - 1)
452         * </pre></p>
453         * <p>
454         * The formula used assumes that the supplied mean value is the weighted arithmetic
455         * mean of the sample data, not a known population parameter. This method
456         * is supplied only to save computation when the mean has already been
457         * computed.</p>
458         * <p>
459         * This formula will not return the same result as the unweighted variance when all
460         * weights are equal, unless all weights are equal to 1. The formula assumes that
461         * weights are to be treated as "expansion values," as will be the case if for example
462         * the weights represent frequency counts. To normalize weights so that the denominator
463         * in the variance computation equals the length of the input vector minus one, use <pre>
464         *   <code>evaluate(values, MathUtils.normalizeArray(weights, values.length), mean); </code>
465         * </pre>
466         * <p>
467         * Returns 0 for a single-value (i.e. length = 1) sample.</p>
468         * <p>
469         * Throws <code>IllegalArgumentException</code> if any of the following are true:
470         * <ul><li>the values array is null</li>
471         *     <li>the weights array is null</li>
472         *     <li>the weights array does not have the same length as the values array</li>
473         *     <li>the weights array contains one or more infinite values</li>
474         *     <li>the weights array contains one or more NaN values</li>
475         *     <li>the weights array contains negative values</li>
476         *     <li>the start and length arguments do not determine a valid array</li>
477         * </ul></p>
478         * <p>
479         * Does not change the internal state of the statistic.</p>
480         *
481         * @param values the input array
482         * @param weights the weights array
483         * @param mean the precomputed weighted mean value
484         * @param begin index of the first array element to include
485         * @param length the number of elements to include
486         * @return the variance of the values or Double.NaN if length = 0
487         * @throws IllegalArgumentException if the parameters are not valid
488         * @since 2.1
489         */
490        public double evaluate(final double[] values, final double[] weights,
491                               final double mean, final int begin, final int length) {
492    
493            double var = Double.NaN;
494    
495            if (test(values, weights, begin, length)) {
496                if (length == 1) {
497                    var = 0.0;
498                } else if (length > 1) {
499                    double accum = 0.0;
500                    double dev = 0.0;
501                    double accum2 = 0.0;
502                    for (int i = begin; i < begin + length; i++) {
503                        dev = values[i] - mean;
504                        accum += weights[i] * (dev * dev);
505                        accum2 += weights[i] * dev;
506                    }
507    
508                    double sumWts = 0;
509                    for (int i = 0; i < weights.length; i++) {
510                        sumWts += weights[i];
511                    }
512    
513                    if (isBiasCorrected) {
514                        var = (accum - (accum2 * accum2 / sumWts)) / (sumWts - 1.0);
515                    } else {
516                        var = (accum - (accum2 * accum2 / sumWts)) / sumWts;
517                    }
518                }
519            }
520            return var;
521        }
522    
523        /**
524         * <p>Returns the weighted variance of the values in the input array, using
525         * the precomputed weighted mean value.</p>
526         * <p>
527         * Uses the formula <pre>
528         *   &Sigma;(weights[i]*(values[i] - mean)<sup>2</sup>)/(&Sigma;(weights[i]) - 1)
529         * </pre></p>
530         * <p>
531         * The formula used assumes that the supplied mean value is the weighted arithmetic
532         * mean of the sample data, not a known population parameter. This method
533         * is supplied only to save computation when the mean has already been
534         * computed.</p>
535         * <p>
536         * This formula will not return the same result as the unweighted variance when all
537         * weights are equal, unless all weights are equal to 1. The formula assumes that
538         * weights are to be treated as "expansion values," as will be the case if for example
539         * the weights represent frequency counts. To normalize weights so that the denominator
540         * in the variance computation equals the length of the input vector minus one, use <pre>
541         *   <code>evaluate(values, MathUtils.normalizeArray(weights, values.length), mean); </code>
542         * </pre>
543         * <p>
544         * Returns 0 for a single-value (i.e. length = 1) sample.</p>
545         * <p>
546         * Throws <code>IllegalArgumentException</code> if any of the following are true:
547         * <ul><li>the values array is null</li>
548         *     <li>the weights array is null</li>
549         *     <li>the weights array does not have the same length as the values array</li>
550         *     <li>the weights array contains one or more infinite values</li>
551         *     <li>the weights array contains one or more NaN values</li>
552         *     <li>the weights array contains negative values</li>
553         * </ul></p>
554         * <p>
555         * Does not change the internal state of the statistic.</p>
556         *
557         * @param values the input array
558         * @param weights the weights array
559         * @param mean the precomputed weighted mean value
560         * @return the variance of the values or Double.NaN if length = 0
561         * @throws IllegalArgumentException if the parameters are not valid
562         * @since 2.1
563         */
564        public double evaluate(final double[] values, final double[] weights, final double mean) {
565            return evaluate(values, weights, mean, 0, values.length);
566        }
567    
568        /**
569         * @return Returns the isBiasCorrected.
570         */
571        public boolean isBiasCorrected() {
572            return isBiasCorrected;
573        }
574    
575        /**
576         * @param biasCorrected The isBiasCorrected to set.
577         */
578        public void setBiasCorrected(boolean biasCorrected) {
579            this.isBiasCorrected = biasCorrected;
580        }
581    
582        /**
583         * {@inheritDoc}
584         */
585        @Override
586        public Variance copy() {
587            Variance result = new Variance();
588            copy(this, result);
589            return result;
590        }
591    
592        /**
593         * Copies source to dest.
594         * <p>Neither source nor dest can be null.</p>
595         *
596         * @param source Variance to copy
597         * @param dest Variance to copy to
598         * @throws NullPointerException if either source or dest is null
599         */
600        public static void copy(Variance source, Variance dest) {
601            if (source == null ||
602                dest == null) {
603                throw new NullArgumentException();
604            }
605            dest.setData(source.getDataRef());
606            dest.moment = source.moment.copy();
607            dest.isBiasCorrected = source.isBiasCorrected;
608            dest.incMoment = source.incMoment;
609        }
610    }