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!

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