stat.h

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

/* 
 * 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
 *
 */

/*!
 * \file stat.h
 * \brief Definitions for functions that perform common
 * statistical operations on Scythe Matrix objects.
 *
 * \note As is the case throughout the library, we provide both
 * general and default template definitions of the Matrix-returning
 * functions in this file, explicitly providing documentation for only
 * the general template versions.
 */

#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 */
  
  /*! 
	* \brief Calculate the sum of a Matrix
	*
	* This function calculates the sum of a matrix by adding each element
	* in turn.
	*
	* \param A The matrix to be summed.
	*
	* \see prod(const Matrix<T,PO,PS> &A)
	* \see sumc(const Matrix<T,PO,PS> &A)
	* \see prodc(const Matrix<T,PO,PS> &A)
	*/


  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 */
   
  /*! 
	* \brief Calculate the sum of each column in a Matrix
	*
	* This function calculates the sum of each column in a matrix by 
	* consecutively adding elements in a single column, looping through all 
	* columns, and returning the results.
	*
	* \param A The matrix to be summed.
	*
	* \see prod(const Matrix<T,PO,PS> &A)
	* \see sum(const Matrix<T,PO,PS> &A)
	* \see prodc(const Matrix<T,PO,PS> &A)
	*/

  SCYTHE_STATMETH_COL(sum)
  
  /* Calculate the product of a Matrix */
  
   /*! 
	* \brief Calculate the product of a Matrix
	*
	* This function calculates the product of a matrix by beginning with the 
	* first element of a matrix, and consecutively multiplying each entry.
	*
	* \param A The matrix to be multiplied.
	*
	* \see sumc(const Matrix<T,PO,PS> &A)
	* \see sum(const Matrix<T,PO,PS> &A)
	* \see prodc(const Matrix<T,PO,PS> &A)
	*/
	
  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 */
  
   /*! 
	* \brief Calculate the product of each column of a Matrix
	*
	* This function calculates the product of each column of a matrix by 
	* multiplying all elements of a single column, looping through all columns,
	* and returning the results.
	*
	* \param A The matrix to be multiplied.
	*
	* \see sumc(const Matrix<T,PO,PS> &A)
	* \see sum(const Matrix<T,PO,PS> &A)
	* \see prod(const Matrix<T,PO,PS> &A)
	*/
	
  SCYTHE_STATMETH_COL(prod)
  
  /* Calculate the mean of a Matrix */
    
   /*! 
	* \brief Calculate the mean of a Matrix
	*
	* This function calculates the mean of a matrix by summing all elements of 
	* the matrix, and dividing by the total number of elements in the matrix.
	*
	* \param A The matrix to be averaged.
	*
	* \see sum(const Matrix<T,PO,PS> &A)
	* \see meanc(const Matrix<T,PO,PS> &A)
	* \see median(const Matrix<T,PO,PS> &A)
	* \see mode(const Matrix<T,PO,PS> &A)
	* \see variance(const Matrix<T,PO,PS> &A)
	*/

  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 */
  
   /*! 
	* \brief Calculate the mean of each column of a Matrix
	*
	* This function calculates the mean of each column of a matrix by summing 
	* all elements of a column in the matrix, divding by the total number of 
	* elements in the column, and looping over every column in the matrix.
	*
	* \param A The matrix to be averaged.
	*
	* \see sumc(const Matrix<T,PO,PS> &A)
	* \see mean(const Matrix<T,PO,PS> &A)
	* \see medianc(const Matrix<T,PO,PS> &A)
	* \see modec(const Matrix<T,PO,PS> &A)
	* \see variancec(const Matrix<T,PO,PS> &A)
	*/

  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
   */
   
   /*! 
	* \brief Calculate the median of a Matrix
	*
	* This function calculates the median of a matrix by first sorting the elements
	* of the matrix, and then finding the middle element.
	*
	* \param A The matrix whose median is of interest.
	*
	* \see medianc(const Matrix<T,PO,PS> &A)
	* \see mean(const Matrix<T,PO,PS> &A)
	* \see mode(const Matrix<T,PO,PS> &A)
	*/

  template <typename T, matrix_order PO, matrix_style PS>
  T
  median (const Matrix<T,PO,PS> &A)
  {
    Matrix<T, PO, PS> 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 */
  
   /*! 
	* \brief Calculate the median of each column a Matrix
	*
	* This function calculates the median of each column of a matrix by first 
	* sorting the elements and locating the middle in a single column, and then
	* looping over all columns.
	*
	* \param A The matrix whose medians are of interest.
	*
	* \see median(const Matrix<T,PO,PS> &A)
	* \see meanc(const Matrix<T,PO,PS> &A)
	* \see modec(const Matrix<T,PO,PS> &A)
	*/

  SCYTHE_STATMETH_COL(median)

  /* Calculate the mode of a matrix */
  
   /*! 
	* \brief Calculate the mode of a Matrix
	*
	* This function calculates the mode of a matrix by determining which value of
	* the matrix occurs with the highest frequency.
	*
	* \param A The matrix whose mode is of interest.
	*
	* \see modec(const Matrix<T,PO,PS> &A)
	* \see mean(const Matrix<T,PO,PS> &A)
	* \see median(const Matrix<T,PO,PS> &A)
	*/
	
  template <typename T, matrix_order PO, matrix_style PS>
  T
  mode (const Matrix<T,PO,PS> &A)
  {
    Matrix<T, PO, PS> 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;
  }

   /*! 
	* \brief Calculate the mode of the columns of a Matrix
	*
	* This function calculates the mode of the columns of a matrix by 
	* determining which value in a single column of the matrix occurs 
	* most frequently, and then looping over all columns.
	*
	* \param A The matrix whose modes are of interest.
	*
	* \see mode(const Matrix<T,PO,PS> &A)
	* \see meanc(const Matrix<T,PO,PS> &A)
	* \see medianc(const Matrix<T,PO,PS> &A)
	*/

  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_);
      }