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 */
017package org.apache.commons.math4.legacy.stat.descriptive.moment;
018
019import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
020import org.apache.commons.math4.legacy.exception.NullArgumentException;
021import org.apache.commons.math4.legacy.stat.descriptive.AbstractStorelessUnivariateStatistic;
022import org.apache.commons.math4.core.jdkmath.JdkMath;
023import org.apache.commons.math4.legacy.core.MathArrays;
024
025
026/**
027 * Computes the Kurtosis of the available values.
028 * <p>
029 * We use the following (unbiased) formula to define kurtosis:</p>
030 * <p>
031 * kurtosis = { [n(n+1) / (n -1)(n - 2)(n-3)] sum[(x_i - mean)^4] / std^4 } - [3(n-1)^2 / (n-2)(n-3)]
032 * </p><p>
033 * where n is the number of values, mean is the {@link Mean} and std is the
034 * {@link StandardDeviation}</p>
035 * <p>
036 * Note that this statistic is undefined for {@code n < 4}.  <code>Double.Nan</code>
037 * is returned when there is not sufficient data to compute the statistic.
038 * Note that Double.NaN may also be returned if the input includes NaN
039 * and / or infinite values.</p>
040 * <p>
041 * <strong>Note that this implementation is not synchronized.</strong> If
042 * multiple threads access an instance of this class concurrently, and at least
043 * one of the threads invokes the <code>increment()</code> or
044 * <code>clear()</code> method, it must be synchronized externally.</p>
045 */
046public class Kurtosis extends AbstractStorelessUnivariateStatistic {
047    /**Fourth Moment on which this statistic is based. */
048    protected FourthMoment moment;
049
050    /**
051     * Determines whether or not this statistic can be incremented or cleared.
052     * <p>
053     * Statistics based on (constructed from) external moments cannot
054     * be incremented or cleared.</p>
055     */
056    protected boolean incMoment;
057
058    /**
059     * Construct a Kurtosis.
060     */
061    public Kurtosis() {
062        incMoment = true;
063        moment = new FourthMoment();
064    }
065
066    /**
067     * Construct a Kurtosis from an external moment.
068     *
069     * @param m4 external Moment
070     */
071    public Kurtosis(final FourthMoment m4) {
072        incMoment = false;
073        this.moment = m4;
074    }
075
076    /**
077     * Copy constructor, creates a new {@code Kurtosis} identical
078     * to the {@code original}.
079     *
080     * @param original the {@code Kurtosis} instance to copy
081     * @throws NullArgumentException if original is null
082     */
083    public Kurtosis(Kurtosis original) throws NullArgumentException {
084        copy(original, this);
085    }
086
087    /**
088     * {@inheritDoc}
089     * <p>Note that when {@link #Kurtosis(FourthMoment)} is used to
090     * create a Variance, this method does nothing. In that case, the
091     * FourthMoment should be incremented directly.</p>
092     */
093    @Override
094    public void increment(final double d) {
095        if (incMoment) {
096            moment.increment(d);
097        }
098    }
099
100    /**
101     * {@inheritDoc}
102     */
103    @Override
104    public double getResult() {
105        double kurtosis = Double.NaN;
106        if (moment.getN() > 3) {
107            double variance = moment.m2 / (moment.n - 1);
108                if (moment.n <= 3 || variance < 10E-20) {
109                    kurtosis = 0.0;
110                } else {
111                    double n = moment.n;
112                    kurtosis =
113                        (n * (n + 1) * moment.getResult() -
114                                3 * moment.m2 * moment.m2 * (n - 1)) /
115                                ((n - 1) * (n -2) * (n -3) * variance * variance);
116                }
117        }
118        return kurtosis;
119    }
120
121    /**
122     * {@inheritDoc}
123     */
124    @Override
125    public void clear() {
126        if (incMoment) {
127            moment.clear();
128        }
129    }
130
131    /**
132     * {@inheritDoc}
133     */
134    @Override
135    public long getN() {
136        return moment.getN();
137    }
138
139    /* UnivariateStatistic Approach  */
140
141    /**
142     * Returns the kurtosis of the entries in the specified portion of the
143     * input array.
144     * <p>
145     * See {@link Kurtosis} for details on the computing algorithm.</p>
146     * <p>
147     * Throws <code>IllegalArgumentException</code> if the array is null.</p>
148     *
149     * @param values the input array
150     * @param begin index of the first array element to include
151     * @param length the number of elements to include
152     * @return the kurtosis of the values or Double.NaN if length is less than 4
153     * @throws MathIllegalArgumentException if the input array is null or the array
154     * index parameters are not valid
155     */
156    @Override
157    public double evaluate(final double[] values, final int begin, final int length)
158        throws MathIllegalArgumentException {
159
160        // Initialize the kurtosis
161        double kurt = Double.NaN;
162
163        if (MathArrays.verifyValues(values, begin, length) && length > 3) {
164            // Compute the mean and standard deviation
165            Variance variance = new Variance();
166            variance.incrementAll(values, begin, length);
167            double mean = variance.moment.m1;
168            double stdDev = JdkMath.sqrt(variance.getResult());
169
170            // Sum the ^4 of the distance from the mean divided by the
171            // standard deviation
172            double accum3 = 0.0;
173            for (int i = begin; i < begin + length; i++) {
174                accum3 += JdkMath.pow(values[i] - mean, 4.0);
175            }
176            accum3 /= JdkMath.pow(stdDev, 4.0d);
177
178            // Get N
179            double n0 = length;
180
181            double coefficientOne =
182                (n0 * (n0 + 1)) / ((n0 - 1) * (n0 - 2) * (n0 - 3));
183            double termTwo =
184                (3 * JdkMath.pow(n0 - 1, 2.0)) / ((n0 - 2) * (n0 - 3));
185
186            // Calculate kurtosis
187            kurt = (coefficientOne * accum3) - termTwo;
188        }
189        return kurt;
190    }
191
192    /**
193     * {@inheritDoc}
194     */
195    @Override
196    public Kurtosis copy() {
197        Kurtosis result = new Kurtosis();
198        // No try-catch because args are guaranteed non-null
199        copy(this, result);
200        return result;
201    }
202
203    /**
204     * Copies source to dest.
205     * <p>Neither source nor dest can be null.</p>
206     *
207     * @param source Kurtosis to copy
208     * @param dest Kurtosis to copy to
209     * @throws NullArgumentException if either source or dest is null
210     */
211    public static void copy(Kurtosis source, Kurtosis dest)
212        throws NullArgumentException {
213        NullArgumentException.check(source);
214        NullArgumentException.check(dest);
215        dest.moment = source.moment.copy();
216        dest.incMoment = source.incMoment;
217    }
218}