[R] always NaN after some running in R, but all fine in S-plus

[Yingfu Xie]

Hello, S-plus and R helpers,(sorry for cross-post)

I wrote some simple C code for one likelihood to be optimized (using
optim(MASS)). I use same function, same data, same starting points and same
DLL in R and S-plus for comparison. (I compiled it with 'Rcmd SHLIB
likelihood.c' and the header files of it include only R.h and math.h). While
it works quite fine in S-plus, it forever returns NaN for the value of
log-likelihood after several normal running in R. What is worse, it can
never return a normal value once NaN appears. I just cannot understand why.
Has anyone had similar experience? I searched R and S-plus archive for
nothing. I attach more details below and really appreciate any help from
you!

R:1.9.0  S-plus:6.2 OS: Windows 2000

Best Regards,
Yingfu

The output in R are:

> symbol.C("like")
[1] "like"
> likelihood(data=data03) # The likelihood function calling .C("like",...)
[1] 5850.12
> myoptim(data=data03) # Optimization using optim with fn=likelihood

Initial parameters are
1 10 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0.2 0.3 0.5 
1 10 0.2 0.2 0.2[1] 5850.154     #I print first 5 parameters and the value 
1.001 10 0.2 0.2 0.2[1] 5850.308 #of the minus log-likelihood every run
0.999 10 0.2 0.2 0.2[1] 5849.997
...........                      #similar outputs are omitted
1 10 0.2 0.2 0.2[1] 5848.578
0 0 0 0 0[1] NaN                 #It seems that here it jumps to the edge of
0.001 0 0 0 0[1] NaN             #box constrain, all zeros in my case.
.............                    #And it returns NaN. From now on, it is  
0 0 0 0 0[1] NaN                 # forever NaN,
0 0 0 0 0[1] NaN                 # even if it jumps back near to the initial

0.995734 9.95734 0.1991468 0.1991468 0.1991468[1] NaN    #parameters.
0.996734 9.95734 0.1991468 0.1991468 0.1991468[1] NaN
.............
1 10 0.2 0.2 0.201[1] NaN
1 10 0.2 0.2 0.199[1] NaN
$par
 [1]  1.0 10.0  0.2  0.2  0.2  0.1  0.1  0.1  0.1  0.2  0.3  0.5

$value          # Since NaN is not allowed with method L-BFGS-B, I set it 
[1] 1e+05       # returns a bigger 10000 if the value is NaN or NA.

$counts
function gradient 
       7        7 

$convergence
[1] 0

What is more, it became worse now:
> likelihood(data=data03) #It returns a NaN also!!
[1] NaN
> myoptim(data=data03)
Initial parameters are
1 10 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0.2 0.3 0.5 
1 10 0.2 0.2 0.2[1] NaN
1.001 10 0.2 0.2 0.2[1] NaN
................
1 10 0.2 0.2 0.2[1] NaN
$par
 [1]  1.0 10.0  0.2  0.2  0.2  0.1  0.1  0.1  0.1  0.2  0.3  0.5

$value
[1] 1e+05

$counts
function gradient 
       1        1 

$convergence
[1] 0

$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

Now, the problem seems lie in my C code and algorithm. However, in S-plus,
it works fine:

> is.loaded("like")
[1] T
> myoptim(data=data03)
Initial parameters are
 1 10 0.2 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0.3 0.5 
 [1]  1.0 10.0  0.2  0.2  0.2 [1] 5834.419
 [1]  1.001 10.000  0.200  0.200 [1] 5834.579
..............                                  #it never jumps to zeros
[1] 0.7701213 4.3442599 0.1424875 0.3072238 0.1274840 [1] 5106.236
$par:
 [1] 0.77012134 4.34425988 0.14248754 0.30722383 0.12748400
 [6] 0.48420116 0.00000000 0.02095689 0.00000000 0.61156935
[11] 0.77179635 0.00000000

$value:
[1] 5106.049

$counts:
 function gradient 
       21       21

$convergence:
[1] 0

$message:
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

I tried the 'bad' parameters in S-plus:
> likelihood(par=rep(0,12),data=data03) #It is NA.
[1] NA
> likelihood(data=data03)               #It is normal again.
[1] 5834.421

Thank you very much for your time!
###########################################

This message has been scanned by F-Secure\ Anti-Virus for Mi...{{dropped}}