4
$\begingroup$

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
$\endgroup$
1
  • 4
    $\begingroup$ would this not be a better question to ask stackoverflow? $\endgroup$
    – ACD
    Jan 16 '13 at 4:43
3
$\begingroup$

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. :)

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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

$\endgroup$
1
  • 2
    $\begingroup$ Generally speaking on stack exchange, longer answers are preferred. To improve your answer you could add things such as more detail as to what the underlying problem with log(gamma(par)) is, why lgamma(par) should fix (or work around) the problem, and how to avoid this sort of problem in the future. $\endgroup$
    – Mark Booth
    Jan 21 '13 at 16:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy