CODEWITHSUNDEEP
rfunctionoptimizationparameterslog-likelihood
WiseRohin
WiseRohin

Reputation: 25

Optim() function in R

I am trying to optimize the parameters of a function with a bound using R's optim() function. Here is the function: Function, because mathjax does not work The code snippet below is my implementation of its Log-Likelihood function, which is here

likelihood_levy <- function(params, x){
  mu <- params[1]
  sigma <- params[2]
  n <- length(x)

  ans <- ((1.5 * n) * log(sigma)) +
      sum(log(1 / (x - mu))) -
      (n * log(sqrt(2*pi) * sigma)) - 
      (0.5 * sigma * sum(1 / (x - mu)))
  }

  return(ans)
}

optim(c(-10, 2), likelihood_levy, x = sample1, method="L-BFGS-B",
      lower = c(min(sample1) - 0.001, 0))

How do I add parameter bounds to the function definition? For example, the first parameter must be lower than x and the second shouldn't be lower than 0.

Note: If I use the "L-BFGS-B" method of optim() function, I get an error that says "L-BFGS-B needs finite values of 'fn'" Any solutions to this problem would be great too!

Thanks!

Upvotes: 1

Views: 2532

Answers (2)

G. Grothendieck
G. Grothendieck

Reputation: 269371

There are multiple problems:

  1. There is an extraneous right brace bracket just before the return statement.
  2. Even if lower ensures that x - mu is positive we can still have problems when the numeric gradient is calculated so use a derivative free method (which we do below) or provide a gradient function to optim.
  3. optim calculates the minimum but we want the maximum likelihood, not the minimum likelihood. Tell optim to maximize (which we do below) or use the negative log likelihood instead of the log likelihood.
  4. The starting value should be feasible, i.e. it should be such that the log likelihood can be calculated at it and produce a finite value while satisfying the constraint that mu < min(x).
  5. We need to handle the situation where mu < min(x) does not hold.
  6. sample1 is missing so we can't reproduce the problem in the question.

Fixing all these we have

likelihood_levy <- function(params, x){
  mu <- params[1]
  sigma <- params[2]
  n <- length(x)
  ans <- if (mu >= min(x)) -Inf
    else ((1.5 * n) * log(sigma)) +
      sum(log(1 / (x - mu))) -
      (n * log(sqrt(2*pi) * sigma)) - 
      (0.5 * sigma * sum(1 / (x - mu)))
  return(ans)
}

set.seed(123)
sample1 <- rnorm(25)
st <- c(min(sample1) - 0.001, 1)
res <- optim(st, likelihood_levy, x = sample1, control = list(fnscale = -1))

str(res)

giving

List of 5
 $ par        : num [1:2] -2.25 1.61
 $ value      : num -46.6
 $ counts     : Named int [1:2] 45 NA
  ..- attr(*, "names")= chr [1:2] "function" "gradient"
 $ convergence: int 0
 $ message    : NULL

Upvotes: 2

Robert Long
Robert Long

Reputation: 6802

The problem is that the function likelihood_levy is either returning a value greater than .Machine$double.xmax or within the function, a division by zero is occurring. This line may be the culprit:

sum(log(1 / (x - mu)))

if x is equal to, or just very close to mu, the problem would occur. You could add some defensive code to check for this. However:

the first parameter must be lower than x

I don't think there is a way to do that with optim

Upvotes: 0

Related Questions