How to avoid NaN in optim?

Question

Suppose, I have a function and want to optimize it. But if I use optim() which gives warnings(). How can I avoid these warnings of NaN?

    myfun<-function(par, x){

    f<- sum(x)*length(x)+sum(log(gamma(par))*x)+1
    return(-f)
    }
    optim(0.1, myfun, x=c(1,5,4,7,8,5,6,5,45,8))

$par
[1] 4.203895e-46

$value
    [1] -10762.39
    $counts
function gradient 
     502       NA 
$convergence
    [1] 1
    $message
NULL

There were 50 or more warnings (use warnings() to see the first 50)


warnings()

2: In log(gamma(par)) : NaNs produced
3: In log(gamma(par)) : NaNs produced
4: In log(gamma(par)) : NaNs produced

Answer 1

As AdamO said, you have an issue with non-positive values. Purely for optimization purposes, exponentiating to make sure things are non-negative and adding some salt to avoid zero generally does the trick. So... something like:

myfun<-function(par, x){
 par <- exp(par) + 10^-10
 f<- sum(x)*length(x)+sum(log(gamma(par))*x)+1
 return(-f)
}

Then none of the two optimizers will give you any trouble.

optim( log(0.1), myfun, x=c(1,5,4,7,8,5,6,5,45,8), method="BFGS")
optim( log(0.1), myfun, x=c(1,5,4,7,8,5,6,5,45,8), method="CG")
# I logged because I exponentiate in the function.

Basically you have a constrained optimization problem and you want to express it as an unconstrained one. Therefore you exponentiate your solution space to make sure it is non-negative one.

Just a word of caution: Exponentiation can occasionally lead to problems if you test for very large values, eg. par=113; would lead to evaluating "gamma(exp(113))" which is actually quite a large number. :)

Answer 2

The problem is that your optimal solution is a boundary value (at 0, and negative values of par aren't allowed). You can verify this by plotting myfun with the given values of $x$. I don't think general optimizers are well configured to handle such issues.

Answer 3

Use lgamma(par) instead of log(gamma(par)).