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