> # numdif.r October 2007, 2009 > # > # numerical differentiation demonstration > # > # compute first and second derivatives of log(gamma(x)) at x=1/2 > # > h <- 2**(-(2:32)) > x <- 1/2 > fx <- lgamma(x) > fxph <- lgamma(x+h) > fxmh <- lgamma(x-h) > # now numerical derivatives > fp1 <- (fxph - fx)/h # one sided -- first difference > fp2 <- (fxph - fxmh)/(2*h) # central difference > # second derivatives -- how are these different? > fpp1 <- (fxph - 2*fx + fxmh)/(h*h) > fpp2 <- ((fxph - fx) - (fx - fxmh))/(h*h) > abs(fpp1-fpp2) [1] 8.881784e-16 3.552714e-15 0.000000e+00 0.000000e+00 0.000000e+00 [6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [26] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [31] 0.000000e+00 > # > # targets > digamma(x) # first derivative [1] -1.96351 > trigamma(x) # second derivative [1] 4.934802 > # table > cbind(h,fp1,fp2,fpp1,fpp2) h fp1 fp2 fpp1 fpp2 [1,] 2.500000e-01 -1.476336 -2.169483 5.545177 5.545177 [2,] 1.250000e-01 -1.692284 -2.008978 5.067110 5.067110 [3,] 6.250000e-02 -1.819352 -1.974565 4.966841 4.966841 [4,] 3.125000e-02 -1.889025 -1.966255 4.942750 4.942750 [5,] 1.562500e-02 -1.925627 -1.964195 4.936785 4.936785 [6,] 7.812500e-03 -1.944403 -1.963681 4.935298 4.935298 [7,] 3.906250e-03 -1.953914 -1.963553 4.934926 4.934926 [8,] 1.953125e-03 -1.958702 -1.963521 4.934833 4.934833 [9,] 9.765625e-04 -1.961103 -1.963513 4.934810 4.934810 [10,] 4.882812e-04 -1.962306 -1.963511 4.934804 4.934804 [11,] 2.441406e-04 -1.962908 -1.963510 4.934803 4.934803 [12,] 1.220703e-04 -1.963209 -1.963510 4.934802 4.934802 [13,] 6.103516e-05 -1.963359 -1.963510 4.934802 4.934802 [14,] 3.051758e-05 -1.963435 -1.963510 4.934802 4.934802 [15,] 1.525879e-05 -1.963472 -1.963510 4.934803 4.934803 [16,] 7.629395e-06 -1.963491 -1.963510 4.934801 4.934801 [17,] 3.814697e-06 -1.963501 -1.963510 4.934807 4.934807 [18,] 1.907349e-06 -1.963505 -1.963510 4.934814 4.934814 [19,] 9.536743e-07 -1.963508 -1.963510 4.934937 4.934937 [20,] 4.768372e-07 -1.963509 -1.963510 4.935059 4.935059 [21,] 2.384186e-07 -1.963509 -1.963510 4.937500 4.937500 [22,] 1.192093e-07 -1.963510 -1.963510 4.937500 4.937500 [23,] 5.960464e-08 -1.963510 -1.963510 4.968750 4.968750 [24,] 2.980232e-08 -1.963510 -1.963510 5.125000 5.125000 [25,] 1.490116e-08 -1.963510 -1.963510 6.000000 6.000000 [26,] 7.450581e-09 -1.963510 -1.963510 6.000000 6.000000 [27,] 3.725290e-09 -1.963510 -1.963510 8.000000 8.000000 [28,] 1.862645e-09 -1.963510 -1.963510 0.000000 0.000000 [29,] 9.313226e-10 -1.963510 -1.963510 128.000000 128.000000 [30,] 4.656613e-10 -1.963510 -1.963510 512.000000 512.000000 [31,] 2.328306e-10 -1.963510 -1.963510 0.000000 0.000000 > # done > rm(list=ls()) > q() > # spendley.r October 2007 > # > # numerical differentiation using Spendley, et al method > # W Spendley, GR Hext, and FR Himsworth (1962) Technometrics, 4:441 > # > spend <- function(f,x,typx) { + d <- length(x) # dimension + dall <- ((d+1)*(d+2))/2 # all points + h <- typx*1.e-5 # step sizes + X <- matrix(c(rep(0,d),diag(h,d)),d+1,d,byrow=T) # 'simplex' pts + I <- c(matrix(0:d,d+1,d+1,byrow=T)) # make a vector + J <- c(matrix(0:d,d+1,d+1,byrow=F)) + Is <- I[!(I>J)] # indices where I < J + Js <- J[!(I>J)] # length(Is) = dall + BigX <- (X[Is+1,]+X[Js+1,])/2 # design points + y <- apply(matrix(x,length(Is),d,byrow=T)+BigX,1,f) + a0 <- y[1] + y <- y - a0 # centering simplifies + y0i <- y[1+(1:d)] # 0 and i values + yii <- (y[Is==Js])[1+(1:d)] # i and i values + a <- 2*( 2*y0i - yii/2 )/h # gradient + bii <- ( yii - 2*y0i )/(h*h) # diagonal of hessian + Hv <- rep(0,d*d) # fill with zeroes + indx <- (d*(Is-1) + Js )[(d+2):dall] # indices for hessian + Hv[indx] <- y[(d+2):dall] # fill the ones we have + H <- matrix(Hv,d,d) - matrix(y0i,d,d) # lower triangular part + H <- H + t(H) # symmetrize and subtract + H <- diag(1/h,d) %*% H %*% diag(1/h,d) # rescale + diag(H) <- bii # fix diagonal + H <- 4*H # everything off by 4 + spend <- list(f0=a0,grad=a,hess=H) } # end function > # > # define test function > fx <- function(x) (log((x[1]*x[2])/(x[3]*x[4])))**2 > # > # test point > v <- spend(fx,c(1,2,3,4),rep(1,4)) > v $f0 [1] 3.210402 $grad [1] -3.5835189 -1.7917595 1.1945063 0.8958797 $hess [,1] [,2] [,3] [,4] [1,] 5.5834626 0.9999646 -0.6666667 -0.49999116 [2,] 0.9999646 1.3958257 -0.3333334 -0.25002223 [3,] -0.6666667 -0.3333334 -0.1759481 0.16667556 [4,] -0.4999912 -0.2500222 0.1666756 -0.09899637 > # done > rm(list=ls()) > q()