1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/*
2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more
3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements.  See the NOTICE file distributed with
4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership.
5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0
6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with
7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License.  You may obtain a copy of the License at
8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *
9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *      http://www.apache.org/licenses/LICENSE-2.0
10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *
11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software
12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS,
13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and
15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License.
16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */
17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.stat.descriptive.moment;
18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.Serializable;
20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.stat.descriptive.AbstractStorelessUnivariateStatistic;
22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath;
23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/**
25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Computes the skewness of the available values.
26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>
27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * We use the following (unbiased) formula to define skewness:</p>
28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>
29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * skewness = [n / (n -1) (n - 2)] sum[(x_i - mean)^3] / std^3 </p>
30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>
31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * where n is the number of values, mean is the {@link Mean} and std is the
32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link StandardDeviation} </p>
33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>
34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <strong>Note that this implementation is not synchronized.</strong> If
35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * multiple threads access an instance of this class concurrently, and at least
36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * one of the threads invokes the <code>increment()</code> or
37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <code>clear()</code> method, it must be synchronized externally. </p>
38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond *
39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1006299 $ $Date: 2010-10-10 16:47:17 +0200 (dim. 10 oct. 2010) $
40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */
41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class Skewness extends AbstractStorelessUnivariateStatistic implements Serializable {
42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /** Serializable version identifier */
44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    private static final long serialVersionUID = 7101857578996691352L;
45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /** Third moment on which this statistic is based */
47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    protected ThirdMoment moment = null;
48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     /**
50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Determines whether or not this statistic can be incremented or cleared.
51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * <p>
52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Statistics based on (constructed from) external moments cannot
53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * be incremented or cleared.</p>
54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    */
55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    protected boolean incMoment;
56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Constructs a Skewness
59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public Skewness() {
61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        incMoment = true;
62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        moment = new ThirdMoment();
63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Constructs a Skewness with an external moment
67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @param m3 external moment
68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public Skewness(final ThirdMoment m3) {
70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        incMoment = false;
71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        this.moment = m3;
72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Copy constructor, creates a new {@code Skewness} identical
76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * to the {@code original}
77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     *
78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @param original the {@code Skewness} instance to copy
79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public Skewness(Skewness original) {
81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        copy(original, this);
82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * {@inheritDoc}
86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    @Override
88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public void increment(final double d) {
89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        if (incMoment) {
90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            moment.increment(d);
91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        }
92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Returns the value of the statistic based on the values that have been added.
96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * <p>
97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * See {@link Skewness} for the definition used in the computation.</p>
98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     *
99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @return the skewness of the available values.
100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    @Override
102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public double getResult() {
103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        if (moment.n < 3) {
105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            return Double.NaN;
106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        }
107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        double variance = moment.m2 / (moment.n - 1);
108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        if (variance < 10E-20) {
109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            return 0.0d;
110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        } else {
111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            double n0 = moment.getN();
112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            return  (n0 * moment.m3) /
113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            ((n0 - 1) * (n0 -2) * FastMath.sqrt(variance) * variance);
114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        }
115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * {@inheritDoc}
119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public long getN() {
121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        return moment.getN();
122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * {@inheritDoc}
126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    @Override
128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public void clear() {
129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        if (incMoment) {
130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            moment.clear();
131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        }
132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Returns the Skewness of the entries in the specifed portion of the
136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * input array.
137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * <p>
138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * See {@link Skewness} for the definition used in the computation.</p>
139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * <p>
140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Throws <code>IllegalArgumentException</code> if the array is null.</p>
141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     *
142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @param values the input array
143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @param begin the index of the first array element to include
144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @param length the number of elements to include
145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @return the skewness of the values or Double.NaN if length is less than
146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * 3
147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @throws IllegalArgumentException if the array is null or the array index
148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     *  parameters are not valid
149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    @Override
151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public double evaluate(final double[] values,final int begin,
152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            final int length) {
153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        // Initialize the skewness
155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        double skew = Double.NaN;
156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        if (test(values, begin, length) && length > 2 ){
158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            Mean mean = new Mean();
159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            // Get the mean and the standard deviation
160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            double m = mean.evaluate(values, begin, length);
161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            // Calc the std, this is implemented here instead
163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            // of using the standardDeviation method eliminate
164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            // a duplicate pass to get the mean
165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            double accum = 0.0;
166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            double accum2 = 0.0;
167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            for (int i = begin; i < begin + length; i++) {
168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond                final double d = values[i] - m;
169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond                accum  += d * d;
170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond                accum2 += d;
171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            }
172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            final double variance = (accum - (accum2 * accum2 / length)) / (length - 1);
173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            double accum3 = 0.0;
175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            for (int i = begin; i < begin + length; i++) {
176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond                final double d = values[i] - m;
177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond                accum3 += d * d * d;
178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            }
179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            accum3 /= variance * FastMath.sqrt(variance);
180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            // Get N
182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            double n0 = length;
183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            // Calculate skewness
185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond            skew = (n0 / ((n0 - 1) * (n0 - 2))) * accum3;
186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        }
187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        return skew;
188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * {@inheritDoc}
192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    @Override
194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public Skewness copy() {
195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        Skewness result = new Skewness();
196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        copy(this, result);
197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        return result;
198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond
200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    /**
201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * Copies source to dest.
202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * <p>Neither source nor dest can be null.</p>
203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     *
204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @param source Skewness to copy
205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @param dest Skewness to copy to
206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     * @throws NullPointerException if either source or dest is null
207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond     */
208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    public static void copy(Skewness source, Skewness dest) {
209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        dest.setData(source.getDataRef());
210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        dest.moment = new ThirdMoment(source.moment.copy());
211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond        dest.incMoment = source.incMoment;
212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond    }
213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond}
214