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