/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.apache.commons.math4.legacy.stat.inference;
  
  import java.util.ArrayList;
  import java.util.Collection;
  
  import org.apache.commons.statistics.distribution.FDistribution;
  import org.apache.commons.math4.legacy.exception.ConvergenceException;
  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  import org.apache.commons.math4.legacy.exception.NullArgumentException;
  import org.apache.commons.math4.legacy.exception.OutOfRangeException;
  import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  import org.apache.commons.math4.legacy.stat.descriptive.SummaryStatistics;
  
  /**
   * Implements one-way ANOVA (analysis of variance) statistics.
   *
   * <p> Tests for differences between two or more categories of univariate data
   * (for example, the body mass index of accountants, lawyers, doctors and
   * computer programmers).  When two categories are given, this is equivalent to
   * the {@link org.apache.commons.math4.legacy.stat.inference.TTest}.
   * </p><p>
   * Uses the {@link org.apache.commons.statistics.distribution.FDistribution
   * commons-math F Distribution implementation} to estimate exact p-values.</p>
   * <p>This implementation is based on a description at
   * http://faculty.vassar.edu/lowry/ch13pt1.html</p>;
   * <pre>
   * Abbreviations: bg = between groups,
   *                wg = within groups,
   *                ss = sum squared deviations
   * </pre>
   *
   * @since 1.2
   */
  public class OneWayAnova {
  
      /**
       * Default constructor.
       */
      public OneWayAnova() {
      }
  
      /**
       * Computes the ANOVA F-value for a collection of <code>double[]</code>
       * arrays.
       *
       * <p><strong>Preconditions</strong>: <ul>
       * <li>The categoryData <code>Collection</code> must contain
       * <code>double[]</code> arrays.</li>
       * <li> There must be at least two <code>double[]</code> arrays in the
       * <code>categoryData</code> collection and each of these arrays must
       * contain at least two values.</li></ul><p>
       * This implementation computes the F statistic using the definitional
       * formula<pre>
       *   F = msbg/mswg</pre>
       * where<pre>
       *  msbg = between group mean square
       *  mswg = within group mean square</pre>
       * are as defined <a href="http://faculty.vassar.edu/lowry/ch13pt1.html">
       * here</a>
       *
       * @param categoryData <code>Collection</code> of <code>double[]</code>
       * arrays each containing data for one category
       * @return Fvalue
       * @throws NullArgumentException if <code>categoryData</code> is <code>null</code>
       * @throws DimensionMismatchException if the length of the <code>categoryData</code>
       * array is less than 2 or a contained <code>double[]</code> array does not have
       * at least two values
       */
      public double anovaFValue(final Collection<double[]> categoryData)
          throws NullArgumentException, DimensionMismatchException {
  
          AnovaStats a = anovaStats(categoryData);
          return a.f;
      }
  
      /**
       * Computes the ANOVA P-value for a collection of <code>double[]</code>
       * arrays.
       *
       * <p><strong>Preconditions</strong>: <ul>
       * <li>The categoryData <code>Collection</code> must contain
       * <code>double[]</code> arrays.</li>
      * <li> There must be at least two <code>double[]</code> arrays in the
      * <code>categoryData</code> collection and each of these arrays must
      * contain at least two values.</li></ul><p>
      * This implementation uses the
      * {@link org.apache.commons.statistics.distribution.FDistribution
      * commons-math F Distribution implementation} to estimate the exact
      * p-value, using the formula<pre>
      *   p = survivalProbability(F)</pre>
      * where <code>F</code> is the F value and <code>survivalProbability = 1 - cumulativeProbability</code>
      * is the commons-statistics implementation of the F distribution.
      *
      * @param categoryData <code>Collection</code> of <code>double[]</code>
      * arrays each containing data for one category
      * @return Pvalue
      * @throws NullArgumentException if <code>categoryData</code> is <code>null</code>
      * @throws DimensionMismatchException if the length of the <code>categoryData</code>
      * array is less than 2 or a contained <code>double[]</code> array does not have
      * at least two values
      * @throws ConvergenceException if the p-value can not be computed due to a convergence error
      * @throws MaxCountExceededException if the maximum number of iterations is exceeded
      */
     public double anovaPValue(final Collection<double[]> categoryData)
         throws NullArgumentException, DimensionMismatchException,
         ConvergenceException, MaxCountExceededException {
 
         final AnovaStats a = anovaStats(categoryData);
         // No try-catch or advertised exception because args are valid
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
         final FDistribution fdist = FDistribution.of(a.dfbg, a.dfwg);
         return fdist.survivalProbability(a.f);
     }
 
     /**
      * Computes the ANOVA P-value for a collection of {@link SummaryStatistics}.
      *
      * <p><strong>Preconditions</strong>: <ul>
      * <li>The categoryData <code>Collection</code> must contain
      * {@link SummaryStatistics}.</li>
      * <li> There must be at least two {@link SummaryStatistics} in the
      * <code>categoryData</code> collection and each of these statistics must
      * contain at least two values.</li></ul><p>
      * This implementation uses the
      * {@link org.apache.commons.statistics.distribution.FDistribution
      * commons-math F Distribution implementation} to estimate the exact
      * p-value, using the formula<pre>
      *   p = survivalProbability(F)</pre>
      * where <code>F</code> is the F value and <code>survivalProbability = 1 - cumulativeProbability</code>
      * is the commons-statistics implementation of the F distribution.
      *
      * @param categoryData <code>Collection</code> of {@link SummaryStatistics}
      * each containing data for one category
      * @param allowOneElementData if true, allow computation for one catagory
      * only or for one data element per category
      * @return Pvalue
      * @throws NullArgumentException if <code>categoryData</code> is <code>null</code>
      * @throws DimensionMismatchException if the length of the <code>categoryData</code>
      * array is less than 2 or a contained {@link SummaryStatistics} does not have
      * at least two values
      * @throws ConvergenceException if the p-value can not be computed due to a convergence error
      * @throws MaxCountExceededException if the maximum number of iterations is exceeded
      * @since 3.2
      */
     public double anovaPValue(final Collection<SummaryStatistics> categoryData,
                               final boolean allowOneElementData)
         throws NullArgumentException, DimensionMismatchException,
                ConvergenceException, MaxCountExceededException {
 
         final AnovaStats a = anovaStats(categoryData, allowOneElementData);
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
         final FDistribution fdist = FDistribution.of(a.dfbg, a.dfwg);
         return fdist.survivalProbability(a.f);
     }
 
     /**
      * This method calls the method that actually does the calculations (except
      * P-value).
      *
      * @param categoryData
      *            <code>Collection</code> of <code>double[]</code> arrays each
      *            containing data for one category
      * @return computed AnovaStats
      * @throws NullArgumentException
      *             if <code>categoryData</code> is <code>null</code>
      * @throws DimensionMismatchException
      *             if the length of the <code>categoryData</code> array is less
      *             than 2 or a contained <code>double[]</code> array does not
      *             contain at least two values
      */
     private AnovaStats anovaStats(final Collection<double[]> categoryData)
         throws NullArgumentException, DimensionMismatchException {
 
         NullArgumentException.check(categoryData);
 
         final Collection<SummaryStatistics> categoryDataSummaryStatistics =
                 new ArrayList<>(categoryData.size());
 
         // convert arrays to SummaryStatistics
         for (final double[] data : categoryData) {
             final SummaryStatistics dataSummaryStatistics = new SummaryStatistics();
             categoryDataSummaryStatistics.add(dataSummaryStatistics);
             for (final double val : data) {
                 dataSummaryStatistics.addValue(val);
             }
         }
 
         return anovaStats(categoryDataSummaryStatistics, false);
     }
 
     /**
      * Performs an ANOVA test, evaluating the null hypothesis that there
      * is no difference among the means of the data categories.
      *
      * <p><strong>Preconditions</strong>: <ul>
      * <li>The categoryData <code>Collection</code> must contain
      * <code>double[]</code> arrays.</li>
      * <li> There must be at least two <code>double[]</code> arrays in the
      * <code>categoryData</code> collection and each of these arrays must
      * contain at least two values.</li>
      * <li>alpha must be strictly greater than 0 and less than or equal to 0.5.
      * </li></ul><p>
      * This implementation uses the
      * {@link org.apache.commons.statistics.distribution.FDistribution
      * commons-math F Distribution implementation} to estimate the exact
      * p-value, using the formula<pre>
      *   p = survivalProbability(F)</pre>
      * where <code>F</code> is the F value and <code>survivalProbability = 1 - cumulativeProbability</code>
      * is the commons-statistics implementation of the F distribution.
      * <p>True is returned iff the estimated p-value is less than alpha.</p>
      *
      * @param categoryData <code>Collection</code> of <code>double[]</code>
      * arrays each containing data for one category
      * @param alpha significance level of the test
      * @return true if the null hypothesis can be rejected with
      * confidence 1 - alpha
      * @throws NullArgumentException if <code>categoryData</code> is <code>null</code>
      * @throws DimensionMismatchException if the length of the <code>categoryData</code>
      * array is less than 2 or a contained <code>double[]</code> array does not have
      * at least two values
      * @throws OutOfRangeException if <code>alpha</code> is not in the range (0, 0.5]
      * @throws ConvergenceException if the p-value can not be computed due to a convergence error
      * @throws MaxCountExceededException if the maximum number of iterations is exceeded
      */
     public boolean anovaTest(final Collection<double[]> categoryData,
                              final double alpha)
         throws NullArgumentException, DimensionMismatchException,
         OutOfRangeException, ConvergenceException, MaxCountExceededException {
 
         if (alpha <= 0 || alpha > 0.5) {
             throw new OutOfRangeException(
                     LocalizedFormats.OUT_OF_BOUND_SIGNIFICANCE_LEVEL,
                     alpha, 0, 0.5);
         }
         return anovaPValue(categoryData) < alpha;
     }
 
     /**
      * This method actually does the calculations (except P-value).
      *
      * @param categoryData <code>Collection</code> of <code>double[]</code>
      * arrays each containing data for one category
      * @param allowOneElementData if true, allow computation for one catagory
      * only or for one data element per category
      * @return computed AnovaStats
      * @throws NullArgumentException if <code>categoryData</code> is <code>null</code>
      * @throws DimensionMismatchException if <code>allowOneElementData</code> is false and the number of
      * categories is less than 2 or a contained SummaryStatistics does not contain
      * at least two values
      */
     private AnovaStats anovaStats(final Collection<SummaryStatistics> categoryData,
                                   final boolean allowOneElementData)
         throws NullArgumentException, DimensionMismatchException {
 
         NullArgumentException.check(categoryData);
 
         if (!allowOneElementData) {
             // check if we have enough categories
             if (categoryData.size() < 2) {
                 throw new DimensionMismatchException(LocalizedFormats.TWO_OR_MORE_CATEGORIES_REQUIRED,
                                                      categoryData.size(), 2);
             }
 
             // check if each category has enough data
             for (final SummaryStatistics array : categoryData) {
                 if (array.getN() <= 1) {
                     throw new DimensionMismatchException(LocalizedFormats.TWO_OR_MORE_VALUES_IN_CATEGORY_REQUIRED,
                                                          (int) array.getN(), 2);
                 }
             }
         }
 
         int dfwg = 0;
         double sswg = 0;
         double totsum = 0;
         double totsumsq = 0;
         int totnum = 0;
 
         for (final SummaryStatistics data : categoryData) {
 
             final double sum = data.getSum();
             final double sumsq = data.getSumsq();
             final int num = (int) data.getN();
             totnum += num;
             totsum += sum;
             totsumsq += sumsq;
 
             dfwg += num - 1;
             final double ss = sumsq - ((sum * sum) / num);
             sswg += ss;
         }
 
         final double sst = totsumsq - ((totsum * totsum) / totnum);
         final double ssbg = sst - sswg;
         final int dfbg = categoryData.size() - 1;
         final double msbg = ssbg / dfbg;
         final double mswg = sswg / dfwg;
         final double f = msbg / mswg;
 
         return new AnovaStats(dfbg, dfwg, f);
     }
 
     /**
         Convenience class to pass dfbg,dfwg,F values around within OneWayAnova.
         No get/set methods provided.
     */
     private static final class AnovaStats {
 
         /** Degrees of freedom in numerator (between groups). */
         private final int dfbg;
 
         /** Degrees of freedom in denominator (within groups). */
         private final int dfwg;
 
         /** Statistic. */
         private final double f;
 
         /**
          * Constructor.
          * @param dfbg degrees of freedom in numerator (between groups)
          * @param dfwg degrees of freedom in denominator (within groups)
          * @param f statistic
          */
         private AnovaStats(int dfbg, int dfwg, double f) {
             this.dfbg = dfbg;
             this.dfwg = dfwg;
             this.f = f;
         }
     }
 }