>  # 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()