CODEWITHSUNDEEP
roptimization
Bazman
Bazman

Reputation: 2150

Finding the Maximum of a Function with numerical derivatives in R

I wish to numerically find the maximum of the function multiplied by Beta 3 shown on p346 of the following link when tau=30:

http://www.ssc.upenn.edu/~fdiebold/papers/paper49/Diebold-Li.pdf

They give the answer on p347 as 0.0609.

I would like to confirm this numerically in R. I.e. to take the derivative and find the value where it reaches zero.

library(numDeriv)

x <- 30
testh <- function(lambda){ ((1-exp(-lambda*30))/(lambda*30)) - exp(-lambda*30) }
grad_h <- function(lambda){
  val <- grad(testh, lambda)
  return(val^2)
}

OptLam <- optimize(f=grad_h, interval=c(0.0001,120), tol=0.0000000000001)

I take the square of the gradient as I want the minimum to be at zero.

Unfortunately, the answer comes back as Lambda=120!! With lambda at 120 the value of the objective function is 5.36e-12.

By working by hand I can func a lower value of the numerical derivative that is closer to zero (it is also close to the analytical value given above):

grad_h(0.05977604)

## [1] 4.24494e-12

Why is the function above not finding this lower value? I have set the tolerance very high so it should be able to find such this optimal value?

Is it possible to correct the existing method so that it gives the correct answer?

Is there a better way to find the maximum gradient of a function numerically in R? For example is there an optimizer that looks for zero rather than trying to find a minimum of maximum?

Upvotes: 2

Views: 3265

Answers (1)

Rorschach
Rorschach

Reputation: 32426

You can use uniroot to find where the derivative is 0. This might work for you,

grad_h <- function(lambda){
    val=grad(testh,lambda)
    return(val)
}
## The root
res <- uniroot(grad_h, c(0,120), tol=1e-10)

## see it
ls <- seq(0.001, 1, length=1000)
plot(ls, testh(ls), col="salmon")
abline(v=res$root, col="steelblue", lwd=2, lty=2)
text(x=res$root, y=testh(res$root),
     labels=sprintf("(%f, %s)", res$root,
         format(testh(res$root), scientific = T)), adj=-0.1)

enter image description here

Upvotes: 3

Related Questions