optimize.h

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

    *
    * \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 <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2, typename FUNCTOR>
  T
  gradfdifls (FUNCTOR fun, T alpha, const Matrix<T,PO1,PS1>& theta, 
              const Matrix<T,PO2,PS2>& p)

  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");
    SCYTHE_CHECK_10(! p.isColVector(), scythe_dimension_error,
        "p not column vector");

    unsigned int k = theta.size();
    T h = std::sqrt(epsilon<T>()); 
    h = std::sqrt(h);
    //T h = std::sqrt(2.2e-16);

    T deriv;

    for (unsigned int i = 0; i < k; ++i) {
      T temp = alpha + h;
      donothing(temp);
      T e = temp - alpha;
      deriv = (fun(theta + (alpha + e) * p) - fun(theta + alpha * p)) 
              / e;
    }
    
    return deriv;
  }

   /*! \brief Calculate the Jacobian of a function using a forward
    * difference formula.
    *
    * This function numerically calculates the Jacobian of a
    * vector-valued function using a forward difference formula.
    *
    * \param fun The function to calculate the Jacobian of.  This
    * function should both take and return a Matrix (vector) of type
    * T.
    * \param theta The column vector of parameter values at which to 
    * take the Jacobian of \a fun.
    *
    * \see gradfdif(FUNCTOR fun, const Matrix<T,PO,PS>& theta)
    * \see gradfdifls(FUNCTOR fun, T alpha, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
    * \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>
  jacfdif (FUNCTOR fun, const Matrix<T,PO,PS>& theta)
  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");

    Matrix<T,RO> fval = fun(theta);
    unsigned int k = theta.rows();
    unsigned int n = fval.rows();

    T h = std::sqrt(epsilon<T>()); //2.2e-16
    h = std::sqrt(h);

    Matrix<T,RO,RS> J(n,k);
    Matrix<T,RO> e;
    Matrix<T,RO> temp;
    Matrix<T,RO> fthetae;
    Matrix<T,RO> ftheta;
    
    for (int i = 0; i < k; ++i) {
      e = Matrix<T,RO>(k,1);
      e[i] = h;
      temp = theta + e;
      donothing(temp); /// XXX ??
      e = temp - theta;
      fthetae = fun(theta + e);
      ftheta = fun(theta);
      for (unsigned int j = 0; j < n; ++j) {
        J(j,i) = (fthetae[j] - ftheta[j]) / e[i];
      }
    }
   
    return J;
  }

  // default template
  template <typename T, matrix_order PO, matrix_style PS,
            typename FUNCTOR>
  Matrix<T,PO,PS>
  jacfdif (FUNCTOR fun, const Matrix<T,PO,PS>& theta)
  {
    return jacfdif<PO,Concrete>(fun, theta);
  }


   /*! \brief Calculate the Hessian of a function using a central
    * difference formula.
    *
    * This function numerically calculates the Hessian of a
    * vector-valued function using a central difference formula.
    *
    * \param fun The function to calculate the Hessian of.  This
    * function should take a Matrix (vector) of type T and return a
    * single value of type T.
    * \param theta The column vector of parameter values at which to 
    * calculate the Hessian.
    *
    * \see gradfdif(FUNCTOR fun, const Matrix<T,PO,PS>& theta)
    * \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)
    *
    * \throw scythe_dimension_error
    *
    * \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>
  hesscdif (FUNCTOR fun, const Matrix<T,PO,PS>& theta)
  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");
    
    T fval = fun(theta);

    //std::cout << std::endl;
    //std::cout << "hesscdif theta = " << theta << "\n";
    //std::cout << "hesscdif fun(theta) = " << fval << std::endl;

    unsigned int k = theta.rows();

    // stepsize CAREFUL -- THIS IS MACHINE SPECIFIC !!!!
    T h2 = std::sqrt(epsilon<T>());
    //T h2 = (T) 1e-10;
    T h = std::sqrt(h2); 

    Matrix<T, RO, RS> H(k,k);

    //std::cout << "h2 = " << h2 << "  h = " << h << std::endl;

    Matrix<T,RO> ei;
    Matrix<T,RO> ej;
    Matrix<T,RO> temp;

    for (unsigned int i = 0; i < k; ++i) {
      ei = Matrix<T,RO>(k, 1);
      ei[i] = h;
      temp = theta + ei;
      donothing(temp); // XXX Again, I'm baffled
      ei = temp - theta;
      for (unsigned int j = 0; j < k; ++j){
        ej = Matrix<T,RO>(k,1);
        ej[j] = h;
        temp = theta + ej;
        donothing(temp); // XXX and again
        ej = temp - theta;
        
        if (i == j) {
          H(i,i) = ( -fun(theta + 2.0 * ei) + 16.0 *
              fun(theta + ei) - 30.0 * fval + 16.0 *
              fun(theta - ei) -
              fun(theta - 2.0 * ei)) / (12.0 * h2);
        } else {
          H(i,j) = ( fun(theta + ei + ej) - fun(theta + ei - ej)
              - fun(theta - ei + ej) + fun(theta - ei - ej))
            / (4.0 * h2);
        }
      }
    }
       
    //std::cout << "end of hesscdif, H = " << H << "\n";
    return H;
  }

  // default template
  template <typename T, matrix_order PO, matrix_style PS,
            typename FUNCTOR>
  Matrix<T,PO,PS>
  hesscdif (FUNCTOR fun, const Matrix<T,PO,PS>& theta)
  {
    return hesscdif<PO,Concrete>(fun, theta);
  }

   /*! \brief Find the step length that minimizes an implied 1-dimensional function.
    *
    * This function performs a line search to find the step length
    * that approximately minimizes an implied one dimensional
    * function.
    *
    * \param fun The function to minimize.  This function should take
    * one Matrix (vector) argument of type T and return a single value
    * of type T.
    * \param theta A column vector of parameter values that anchor the
    * 1-dimensional function.
    * \param p A direction vector that creates the 1-dimensional
    * function.
    *
    * \see linesearch2(FUNCTOR fun, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p, rng<RNGTYPE>& runif)
    * \see zoom(FUNCTOR fun, T alpha_lo, T alpha_hi, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
    * \see BFGS(FUNCTOR fun, const Matrix<T,PO,PS>& theta, rng<RNGTYPE>& runif, unsigned int maxit, T tolerance, bool trace = false)
    *
    * \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 <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2, typename FUNCTOR>
  T linesearch1 (FUNCTOR fun, const Matrix<T,PO1,PS1>& theta,
                 const Matrix<T,PO2,PS2>& p)
  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");
    SCYTHE_CHECK_10(! p.isColVector(), scythe_dimension_error,
        "p not column vector");

    T alpha_bar = (T) 1.0;
    T rho = (T) 0.9;
    T c   = (T) 0.5;
    T alpha = alpha_bar;
    Matrix<T,PO1> fgrad = gradfdif(fun, theta);

    while (fun(theta + alpha * p) > 
           (fun(theta) + c * alpha * t(fgrad) * p)[0]) {
      alpha = rho * alpha;
    }

    return alpha;
  }

   /*! \brief Find the step length that minimizes an implied 1-dimensional function.
    *
    * This function performs a line search to find the step length
    * that approximately minimizes an implied one dimensional
    * function.
    *
    * \param fun The function to minimize.  This function should take
    * one Matrix (vector) argument of type T and return a single value
    * of type T.
    * \param theta A column vector of parameter values that anchor the
    * 1-dimensional function.
    * \param p A direction vector that creates the 1-dimensional
    * function.
    * \param runif A random uniform number generator function object
    * (an object that returns a random uniform variate on (0,1) when
    * its () operator is invoked).
    *
    * \see linesearch1(FUNCTOR fun, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
    * \see zoom(FUNCTOR fun, T alpha_lo, T alpha_hi, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
    * \see BFGS(FUNCTOR fun, const Matrix<T,PO,PS>& theta, rng<RNGTYPE>& runif, unsigned int maxit, T tolerance, bool trace = false)
    *
    * \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 <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2, typename FUNCTOR,
            typename RNGTYPE>
  T linesearch2 (FUNCTOR fun, const Matrix<T,PO1,PS1>& theta,
                 const Matrix<T,PO2,PS2>& p, rng<RNGTYPE>& runif)
  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");
    SCYTHE_CHECK_10(! p.isColVector(), scythe_dimension_error,
        "p not column vector");