Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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
share|improve this question

migrated from stats.stackexchange.com Jan 16 '13 at 9:36

This question came from our site for people interested in statistics, machine learning, data analysis, data mining, and data visualization.

3  
would this not be a better question to ask stackoverflow? –  ACD Jan 16 '13 at 4:43

3 Answers 3

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

share|improve this answer

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.

share|improve this answer

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

share|improve this answer
1  
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. –  Mark Booth Jan 21 '13 at 16:20

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.