CODEWITHSUNDEEP
r
Miyazaki
Miyazaki

Reputation: 99

function return NaN

I dont understand why the below function return NaN even the answer shouldn't be NaN. I have tried to search similar question but most of them in other programming language that i couldn't understand.

Code: 

a0 = 1.5; a1 = 0.4; b1 = 0.3; g1= 0.7
nu = rep(0,1)
h.new = rep(0,1)
ddp = rep(0,1)

nu[1]=0

h.new[1] = a0/(1-a1-b1)
ddp[1] = 0.5*log(g1)- g1*h.new[1] + 
  nu[1]*(log(nu[1])-1) - log(factorial(nu[1]))+(g1)*nu[1]* 
  (1+log(h.new[1]/nu[1]))

This is the output:

###> ddp
###[1] NaN

Then I do manual calculation, clearly show that ddp is not NaN:

###h.new = 1.5/(1-0.4-0.3) = 5

###ddp 
###= 0.5*log(0.7)- (0.7)*(5) + 
###(0)*(log(0)-1) - log(0!)+(0.7)*(0)* 
###(1+log(5/0))

###= 0.5*log(0.7)- (0.7)*(5)

###= -3.6783

I did know for sum() function, na.rm=TRUE could ommit all the missing value, but in this case, how to modify the code to make it return to the correct answer?

Thanks in advance.

Upvotes: 0

Views: 2528

Answers (2)

Mark
Mark

Reputation: 4537

nu[1] is 0 and the log of 0 is undefined. R evaluates log(0) as -Inf Then you have log(5/0) - what is 5/0? It's undefined. So you have lots of undefined values which will create an overall undefined result. As noted in Ben's answer, the NaN is created when you multiply 0 by Inf or -Inf

However, even if these two elements were to resolve to 0, the code you've written would evaluate to -3.68, not -0.52. So, you're doing something wrong by hand.

Upvotes: 2

Ben Bolker
Ben Bolker

Reputation: 226182

In general 0*Inf gives NaN in R: the computer doesn't/can't know how you want to take the limit. If you know that you want these limits to be zero, you can specify how your computation works in the special case nu==0:

term3 <- if (nu==0) 0 else nu*(log(nu)-1)
term5 <- if (nu==0) 0 else g1*nu*(1+log(h.new/nu))
ddp = 0.5*log(g1)-
    g1*h.new + 
    term3 - 
    log(factorial(nu))+
    term5

You could do this slightly more generally with a zprod function:

zprod <- function(x,y) {
     if (x==0) 0 else x*y
}
ddp = 0.5*log(g1)- g1*h.new + 
      zprod(nu,log(nu)-1) - log(factorial(nu))+
      g1*zprod(nu,1+log(h.new/nu))

Upvotes: 2

Related Questions