1# Copyright (C) 2011 Jelmer Ypma. All Rights Reserved.
2# This code is published under the L-GPL.
3#
4# File:   finite.diff.R
5# Author: Jelmer Ypma
6# Date:   24 July 2011
7#
8# Approximate the gradient of a function using finite differences.
9#
10# Input:
11#        func : calculate finite difference approximation for the gradient of this function
12#        .x : evaluate at this point
13#        indices : calculate finite difference approcximation only for .x-indices in this vector (optional)
14#        stepSize : relative size of the step between x and x+h (optional)
15#       ... : arguments that are passed to the user-defined function (func)
16#
17# Output: matrix with finite difference approximations
18
19
20# http://en.wikipedia.org/wiki/Numerical_differentiation
21finite.diff <- function( func, .x, indices=1:length(.x), stepSize=sqrt( .Machine$double.eps ), ... ) {
22
23    # if we evaluate at values close to 0, we need a different step size
24    stepSizeVec <- pmax( abs(.x), 1 ) * stepSize
25
26    fx <- func( .x, ... )
27    approx.gradf.index <- function(i, .x, func, fx, stepSizeVec, ...) {
28        x_prime <- .x
29        x_prime[i] <- .x[i] + stepSizeVec[i]
30        stepSizeVec[i] <- x_prime[i] - .x[i]
31        fx_prime <- func( x_prime, ... )
32        return( ( fx_prime - fx )/stepSizeVec[i] )
33    }
34    grad_fx <- sapply (indices, approx.gradf.index, .x=.x, func=func, fx=fx, stepSizeVec=stepSizeVec, ... )
35
36    return( grad_fx )
37}
38