Go to the documentation of this file.
3
5
6 namespace stan {
7 namespace math {
8
9 /**
10  * Calculate the value and the gradient of the specified function
11  * at the specified argument using finite difference.
12  *
13  * <p>The functor must implement
14  *
15  * <code>
16  * double
17  * operator()(const
18  * Eigen::Matrix<double, Eigen::Dynamic, 1>&)
19  * </code>
20  *
21  * Error should be on order of epsilon ^ 6.
22  * The reference for this algorithm is:
23  *
24  * De Levie: An improved numerical approximation
25  * for the first derivative, page 3
26  *
27  * This function involves 6 calls to f.
28  *
29  * @tparam F Type of function
30  * @param[in] f Function
31  * @param[in] x Argument to function
32  * @param[out] fx Function applied to argument
34  * @param[in] epsilon perturbation size
35  */
36 template <typename F>
37 void finite_diff_gradient(const F& f, const Eigen::Matrix<double, -1, 1>& x,
38  double& fx, Eigen::Matrix<double, -1, 1>& grad_fx,
39  double epsilon = 1e-03) {
40  using Eigen::Dynamic;
41  using Eigen::Matrix;
42  Matrix<double, Dynamic, 1> x_temp(x);
43
44  int d = x.size();
46
47  fx = f(x);
48
49  for (int i = 0; i < d; ++i) {
50  double delta_f = 0.0;
51
52  x_temp(i) = x(i) + 3.0 * epsilon;
53  delta_f = f(x_temp);
54
55  x_temp(i) = x(i) + 2.0 * epsilon;
56  delta_f -= 9.0 * f(x_temp);
57
58  x_temp(i) = x(i) + epsilon;
59  delta_f += 45.0 * f(x_temp);
60
61  x_temp(i) = x(i) + -3.0 * epsilon;
62  delta_f -= f(x_temp);
63
64  x_temp(i) = x(i) + -2.0 * epsilon;
65  delta_f += 9.0 * f(x_temp);
66
67  x_temp(i) = x(i) + -epsilon;
68  delta_f -= 45.0 * f(x_temp);
69
70  delta_f /= 60 * epsilon;
71
72  x_temp(i) = x(i);