stat.h

Go to the documentation of this file.

/* 
 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin
 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn,
 * and Daniel Pemstein.  All Rights Reserved.
 *
 * This program is free software; you can redistribute it and/or
 * modify under the terms of the GNU General Public License as
 * published by Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.  See the text files
 * COPYING and LICENSE, distributed with this source code, for further
 * information.
 * --------------------------------------------------------------------
 *  scythestat/stat.h
 *
 */

#ifndef SCYTHE_STAT_H
#define SCYTHE_STAT_H

#ifdef SCYTHE_COMPILE_DIRECT
#include "matrix.h"
#include "algorithm.h"
#include "error.h"
#else
#include "scythestat/matrix.h"
#include "scythestat/algorithm.h"
#include "scythestat/error.h"
#endif

#include <numeric>
#include <set>


namespace scythe {

  namespace {
    typedef unsigned int uint;
  }

/* A macro for defining column versions of a function.  That is,
 * when expanded, this macro produces general and default template
 * functions that compute function NAME on each column in a matrix and
 * return a row vector with the results.  We use this to generate
 * column versions of every function in this header file.
 */
#define SCYTHE_STATMETH_COL(NAME)                                     \
  template <matrix_order RO, matrix_style RS, typename T,             \
            matrix_order PO, matrix_style PS>                         \
  Matrix<T,RO,RS>                                                     \
  NAME ## c (const Matrix<T,PO,PS>& A)                                \
  {                                                                   \
    Matrix<T,RO,RS> res (1, A.cols(), false);                         \
                                                                      \
    for (uint j = 0; j < A.cols(); ++j)                               \
      res[j] = NAME(A(_, j));                                         \
                                                                      \
    return res;                                                       \
  }                                                                   \
                                                                      \
  template <typename T, matrix_order O, matrix_style S>               \
  Matrix<T,O,Concrete>                                                \
  NAME ## c (const Matrix<T,O,S>& A)                                  \
  {                                                                   \
    return NAME ## c<O,Concrete>(A);                                  \
  }         


  /* Calculate the sum of a Matrix */
  
  template <typename T, matrix_order PO, matrix_style PS>
  T
  sum (const Matrix<T,PO,PS> &A)
  {
    return (std::accumulate(A.begin_f(), A.end_f(), (T) 0));
  }

  /* Calculate the sum of each column in a Matrix */
   
  SCYTHE_STATMETH_COL(sum)
  
  /* Calculate the product of a Matrix */
  
   
  template <typename T, matrix_order PO, matrix_style PS>
  T
  prod (const Matrix<T,PO,PS> &A)
  {
    return std::accumulate(A.begin_f(), A.end_f(), (T) 1, 
                           std::multiplies<T> ());
  }

  /* Calculate the product of each column of a matrix */
  
  SCYTHE_STATMETH_COL(prod)
  
  /* Calculate the mean of a Matrix */
    
   
  template <typename T, matrix_order PO, matrix_style PS>
  T
  mean (const Matrix<T,PO,PS> &A)
  {
    return (std::accumulate(A.begin_f(), A.end_f(), (T) 0) / A.size());
  }

  /* Calculate the mean of each column of a Matrix */
  
  SCYTHE_STATMETH_COL(mean)
  
  /* Calculate the median of a matrix.  Uses a sort but I'll implement
   * the randomized alg when I figure out how to generalize it to
   * even-length lists
   */
   
   
  template <typename T, matrix_order PO, matrix_style PS>
  T
  median (const Matrix<T,PO,PS> &A)
  {
    Matrix<T, PO, Concrete> temp(A);
    uint n = temp.size();

    sort(temp.begin(), temp.end());
    if (n % 2 == 0)
      return ((temp[n / 2] + temp[n / 2 - 1]) / 2);
    else
      return temp[(uint) ::floor(n / 2)];
  }

  /* Calculate the median of each column of a matrix */
  
  SCYTHE_STATMETH_COL(median)

  /* Calculate the mode of a matrix */
  
   
  template <typename T, matrix_order PO, matrix_style PS>
  T
  mode (const Matrix<T,PO,PS> &A)
  {
    Matrix<T, PO, Concrete> temp(A);
    
    sort(temp.begin(), temp.end());

    T last = temp[0];
    uint cnt = 1;
    T cur_max = temp[0];
    uint max_cnt = 1;
    
    for (uint i = 1; i < temp.size(); ++i) {
      if (last == temp[i]) {
        ++cnt;
      } else {
        last = temp[i];
        cnt = 1;
      }
      if (cnt > max_cnt) {
        max_cnt = cnt;
        cur_max = temp[i];
      }
    }

    return cur_max;
  }

  SCYTHE_STATMETH_COL(mode)

  /* Calculate the variance of a Matrix */

  /* A functor that encapsulates a single variance calculation step.
   * Also used by skew and kurtosis. */
  namespace {
    template <typename T, typename T2>
    struct var_step : std::binary_function<T, T, T>
    {
      T constant_;
      T2 divisor_;
      T exponent_;
      var_step (T c, T2 d, T e) : constant_ (c), divisor_ (d),
                                    exponent_ (e) {}
      T operator() (T last, T x) const
      {
        return (last + std::pow(constant_ - x, exponent_) / divisor_);
      }
    };
  }
  
  template <typename T, matrix_order PO, matrix_style PS>
  T
  var (const Matrix<T,PO,PS> &A)
  {
    return var(A, mean(A));
  }

  /* Calculate the variances of each column of a Matrix. */
  
  SCYTHE_STATMETH_COL(var)

 
  template <typename T, matrix_order PO, matrix_style PS>
  T
  var (const Matrix<T,PO,PS> &A, T mu)
  {
    return std::accumulate(A.begin_f(), A.end_f(), (T) 0, 
                         var_step<T, uint> (mu, A.size() - 1, 2));
  }
  
  /* Calculate the standard deviation of a Matrix (not std cause of namespace std:: */
  
  template <typename T, matrix_order PO, matrix_style PS>
  T
  sd (const Matrix<T,PO,PS> &A)
  {
    return std::sqrt(var(A));
  }
  
  /* Calculate the standard deviation of each column of a Matrix */
  SCYTHE_STATMETH_COL(sd)

  
  template <typename T, matrix_order PO, matrix_style PS>
  T
  sd (const Matrix<T,PO,PS> &A, T mu)
  {
    return std::sqrt(var(A, mu));
  }

  /* Calculate the skew of a Matrix */

  template <typename T, matrix_order PO, matrix_style PS>
  T
  skew (const Matrix<T,PO,PS> &A)
  {
    T mu = mean(A);
    T sde = sd(A, mu);
    return std::accumulate(A.begin_f(), A.end_f(), (T) 0, 
              var_step<T, T> (mu, A.size() * std::pow(sde, 3), 3));
  }

  /* Calculate the skew of each column of a Matrix. */
  
  SCYTHE_STATMETH_COL(skew)
  
  /* Calculate the kurtosis of a Matrix */
    
   
  template <typename T, matrix_order PO, matrix_style PS>
  T
  kurtosis (const Matrix<T,PO,PS> &A)
  {
    T mu = mean(A);
    T sde = sd(A, mu);
    return (std::accumulate(A.begin_f(), A.end_f(), (T) 0, 
              var_step<T, T> (mu, A.size() * std::pow(sde, 4), 4))
            - 3);
  }
  
  /* Calculate the kurtosis of each column of a Matrix. */
  
  SCYTHE_STATMETH_COL(kurtosis)

  /* Calculates the maximum element in a Matrix */
  template <typename T, matrix_order PO, matrix_style PS>
  T
  max (const Matrix<T,PO,PS> &A)
  {
    return *(max_element(A.begin_f(), A.end_f()));
  }

  SCYTHE_STATMETH_COL(max)

  /* Calculates the minimum element in a Matrix */
  
  
  template <typename T, matrix_order PO, matrix_style PS>
  T
  min (const Matrix<T,PO,PS> &A)
  {
    return *(min_element(A.begin_f(), A.end_f()));
  }
  
  SCYTHE_STATMETH_COL(min)

  /* Find the index of the max element */  
  
  template <typename T, matrix_order PO, matrix_style PS>
  unsigned int
  maxind (const Matrix<T,PO,PS> &A)
  {
    return (max_element(A.begin_f(), A.end_f())).get_index();
  }
  
  SCYTHE_STATMETH_COL(maxind)
  
  /* Find the index of the min element */
  
  
  template <typename T, matrix_order PO, matrix_style PS>
  unsigned int
  minind (const Matrix<T,PO,PS> &A)
  {
    return (min_element(A.begin_f(), A.end_f())).get_index();
  }

  SCYTHE_STATMETH_COL(minind)

} // end namespace scythe


#endif /* SCYTHE_STAT_H */