finite_diff_gradient.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUNCTOR_FINITE_DIFF_GRADIENT_HPP
#define STAN_MATH_PRIM_FUNCTOR_FINITE_DIFF_GRADIENT_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
 
namespace stan {
namespace math {
 
template <typename F>
void finite_diff_gradient(const F& f, const Eigen::VectorXd& x, double& fx,
                          Eigen::VectorXd& grad_fx, double epsilon = 1e-03) {
  Eigen::VectorXd x_temp(x);
  int d = x.size();
  grad_fx.resize(d);
 
  fx = f(x);
 
  for (int i = 0; i < d; ++i) {
    double delta_f = 0.0;
 
    x_temp(i) = x(i) + 3.0 * epsilon;
    delta_f = f(x_temp);
 
    x_temp(i) = x(i) + 2.0 * epsilon;
    delta_f -= 9.0 * f(x_temp);
 
    x_temp(i) = x(i) + epsilon;
    delta_f += 45.0 * f(x_temp);
 
    x_temp(i) = x(i) + -3.0 * epsilon;
    delta_f -= f(x_temp);
 
    x_temp(i) = x(i) + -2.0 * epsilon;
    delta_f += 9.0 * f(x_temp);
 
    x_temp(i) = x(i) + -epsilon;
    delta_f -= 45.0 * f(x_temp);
 
    delta_f /= 60 * epsilon;
 
    x_temp(i) = x(i);
    grad_fx(i) = delta_f;
  }
}
}  // namespace math
}  // namespace stan
#endif