optimize.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/optimize.h
 *
 */

/*!
 * \file optimize.h
 * \brief Definitions of functions for doing numerical optimization
 * and related operations.
 *
 * This file contains a number of functions that are useful for
 * numerical optimization and maximum likelihood estimation.  In
 * addition, it contains some basic facilities for evaluating definite
 * integrals.
 *
 * As is the case across Scythe, we provide both general and default
 * template definitions for the functions in this file that return
 * Matrix objects.  The general definitions allow the user to
 * customize the matrix_order and matrix_style of the returned Matrix,
 * while the default versions return concrete matrices of the same
 * matrix_order as the first (or only) Matrix argument to the
 * function.  In cases where we supply these two types of definitions,
 * we explicitly document only the general version, although the
 * default definition will typically appear in the function list
 * below.
 *
 * \note
 * Doxygen has some difficulty dealing with overloaded templates.
 * Under certain circumstances it does not correctly process the
 * definitions of default templates.  In these cases, the definition
 * for the default template will not even appear in the function list.
 * We provide default templates for all of the Matrix-returning
 * functions in this file.
 *
 */

#ifndef SCYTHE_OPTIMIZE_H
#define SCYTHE_OPTIMIZE_H

#ifdef SCYTHE_COMPILE_DIRECT
#include "matrix.h"
#include "algorithm.h"
#include "error.h"
#include "rng.h"
#include "distributions.h"
#include "la.h"
#include "ide.h"
#include "smath.h"
#include "stat.h"
#else
#include "scythestat/matrix.h"
#include "scythestat/algorithm.h"
#include "scythestat/error.h"
#include "scythestat/rng.h"
#include "scythestat/distributions.h"
#include "scythestat/la.h"
#include "scythestat/ide.h"
#include "scythestat/smath.h"
#include "scythestat/stat.h"
#endif

/* We want to use an anonymous namespace to make the following consts
 * and functions local to this file, but mingw doesn't play nice with
 * anonymous namespaces so we do things differently when using the
 * cross-compiler.
 */
#ifdef __MINGW32__
#define SCYTHE_MINGW32_STATIC static
#else
#define SCYTHE_MINGW32_STATIC
#endif

namespace scythe {
#ifndef __MINGW32__
  namespace {
#endif

    /* Functions (private to this file) that do very little... */
    template <typename T, matrix_order O, matrix_style S>
    SCYTHE_MINGW32_STATIC T donothing (const Matrix<T,O,S>& x)
    {
      return (T) 0.0;
    }

    template <typename T>
    SCYTHE_MINGW32_STATIC T donothing (T& x)
    {
      return (T) 0.0;
    }
#ifndef __MINGW32__
  }
#endif


  /* Return the machine epsilon 
   * Notes: Algorithm taken from Sedgewick, Robert. 1992. Algorithms
   * in C++. Addison Wesley. pg. 561
   */
   /*! \brief Compute the machine epsilon.
    *
    * The epsilon function returns the machine epsilon: the smallest
    * number that, when summed with 1, produces a value greater than
    * one.
    */
  template <typename T>
  T
  epsilon()
  {
    T eps, del, neweps;
    del    = (T) 0.5;
    eps    = (T) 0.0;
    neweps = (T) 1.0;
  
    while ( del > 0 ) {
      if ( 1 + neweps > 1 ) {  /* Then the value might be too large */
        eps = neweps;    /* ...save the current value... */
        neweps -= del;    /* ...and decrement a bit */
      } else {      /* Then the value is too small */
        neweps += del;    /* ...so increment it */
      }
      del *= 0.5;      /* Reduce the adjustment by half */
    }

    return eps;
  }
  
   /*! \brief Calculate the definite integral of a function from a to b.
    *
    * This function calculates the definite integral of a univariate
    * function on the interval \f$[a,b]\f$.
    *
    * \param fun The function (or functor) whose definite integral is
    * to be calculated.  This function should both take and return a
    * single argument of type T.
    * \param a The starting value of the interval.
    * \param b The ending value of the interval.
    * \param N The number of subintervals to calculate.  Increasing
    * this number will improve the accuracy of the estimate but will
    * also increase run-time.
    *
    * \throw scythe_invalid_arg (Level 1)
    *
    * \see adaptsimp(FUNCTOR fun, T a, T b, unsigned int& N, double tol = 1e-5)
    * \note
    * Users will typically wish to implement \a fun in terms of a
    * functor.  Using a functor provides a generic way in which to
    * evaluate functions with more than one parameter.  Furthermore,
    * although one can pass a function pointer to this routine,
    * the compiler cannot inline and fully optimize code
    * referenced by function pointers.
    */
  template <typename T, typename FUNCTOR>
  T  intsimp (FUNCTOR fun, T a, T b, unsigned int N)
  {
    SCYTHE_CHECK_10(a > b, scythe_invalid_arg,
        "Lower limit larger than upper");
    
    T I = (T) 0;
    T w = (b - a) / N;
    for (unsigned int i = 1; i <= N; ++i)
      I += w * (fun(a +(i - 1) *w) + 4 * fun(a - w / 2 + i * w) +
          fun(a + i * w)) / 6;
   
    return I;
  }
  
   /*! \brief Calculate the definite integral of a function from a to b.
    *
    * This function calculates the definite integral of a univariate
    * function on the interval \f$[a,b]\f$.
    *
    * \param fun The function (or functor) whose definite integral is 
    * to be calculated.  This function should both take and return a
    * single argument of type T.
    * \param a The starting value of the interval.
    * \param b The ending value of the interval.
    * \param N The number of subintervals to calculate.  Increasing
    * this number will improve the accuracy of the estimate but will
    * also increase run-time.
    * \param tol The accuracy required.  Both accuracy and run-time
    * decrease as this number increases.
    *
    * \throw scythe_invalid_arg (Level 1)
    *
    * \see intsimp(FUNCTOR fun, T a, T b, unsigned int& N)
    *
    * \note
    * Users will typically wish to implement \a fun in terms of a
    * functor.  Using a functor provides a generic way in which to
    * evaluate functions with more than one parameter.  Furthermore,
    * although one can pass a function pointer to this routine,
    * the compiler cannot inline and fully optimize code
    * referenced by function pointers.
    */
  template <typename T, typename FUNCTOR>
  T adaptsimp(FUNCTOR fun, T a, T b, unsigned int N, double tol = 1e-5)
  {
    SCYTHE_CHECK_10(a > b, scythe_invalid_arg,
        "Lower limit larger than upper");

    T I = intsimp(fun, a, b, N);
    if (std::fabs(I - intsimp(fun, a, b, N / 2)) > tol)
      return adaptsimp(fun, a, (a + b) / 2, N, tol)
        + adaptsimp(fun, (a + b) / 2, b, N, tol);

    return I;
  }

   /*! \brief Calculate gradient of a function using a forward
    * difference formula.
    *
    * This function numerically calculates the gradient of a
    * vector-valued function at \a theta using a forward difference
    * formula.
    *
    * \param fun The function to calculate the gradient of.  This
    * function should both take and return a single Matrix (vector) of 
    * type T.
    * \param theta The column vector of values at which to calculate 
    * the gradient of the function.
    *
    * \see gradfdifls(FUNCTOR fun, T alpha, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
    * \see jacfdif(FUNCTOR fun, const Matrix<T,PO,PS>& theta)
    * \see hesscdif(FUNCTOR fun, const Matrix<T,PO,PS>& theta)
    *
    * \throw scythe_dimension_error (Level 1)
    *
    * \note
    * Users will typically wish to implement \a fun in terms of a
    * functor.  Using a functor provides a generic way in which to
    * evaluate functions with more than one parameter.  Furthermore,
    * although one can pass a function pointer to this routine,
    * the compiler cannot inline and fully optimize code
    * referenced by function pointers.
    */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS, typename FUNCTOR>
  Matrix<T, RO, RS>
  gradfdif (FUNCTOR fun, const Matrix<T,PO,PS>& theta)
      
  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");

    unsigned int k = theta.size();
    T h = std::sqrt(epsilon<T>());
    h = std::sqrt(h);

    Matrix<T,RO,RS> grad(k, 1);
    Matrix<T,RO> e;
    Matrix<T,RO> temp;
    for (unsigned int i = 0; i < k; ++i) {
      e = Matrix<T,RO>(k, 1);
      e[i] = h;
      temp = theta + e;
      donothing(temp); // XXX I don't understand this
      e = temp - theta;
      grad[i] = (fun(theta + e) - fun(theta)) / e[i];
    }

    return grad;
  }

  // Default template version
  template <typename T, matrix_order O, matrix_style S, 
            typename FUNCTOR>
  Matrix<T, O, Concrete>
  gradfdif (FUNCTOR fun, const Matrix<T,O,S>& theta)
  {
    return gradfdif<O,Concrete>(fun, theta);
  }

   /*! \brief Calculate the first derivative of the function using
    * a forward difference formula.
    *
    * This function numerically calculates the first derivative of a
    * function with respect to \a alpha at \f$theta + alpha \cdot p\f$
    * using a forward difference formula.  This function is primarily
    * useful for linesearches.
    *
    * \param fun The function to calculate the first derivative of.
    * This function should take a single Matrix<T> argument and return
    * a value of type T.
    * \param alpha Double the step length.
    * \param theta A Matrix (vector) of parameter values at which to
    * calculate the gradient.
    * \param p A direction vector.
    *
    * \see gradfdif(FUNCTOR fun, const Matrix<T,PO,PS>& theta)
    * \see jacfdif(FUNCTOR fun, const Matrix<T,PO,PS>& theta)
    * \see hesscdif(FUNCTOR fun, const Matrix<T,PO,PS>& theta)