R help - Jun 2008 - L-BFGS-B needs finite values of 'fn'

Hi,

When I run the following code,

r <- c(3,4,4,3,5,4,5,9,8,11,12,13)
n <- rep(15,12)
x <- c(0, 1.1, 1.3, 2.0, 2.2, 2.8, 3.7, 3.9, 4.4, 4.8, 5.9, 6.8)
x <- log10(x)
fr <- function(c, alpha, beta) {
  P <- c + (1-c) * pnorm(alpha + beta * x)
  P <- pmax(pmin(P,1),0)
  -(sum(log(choose(n,r))) + sum(r * log(P)) + sum((n -r)* log(1-P)))
 }
fit <- mle((fr), start = list(c =0.2, alpha = 0, beta =0.1), method
"L-BFGS-B")

I got:

Error in optim(start, f, method = method, hessian = TRUE, ...) :
  L-BFGS-B needs finite values of 'fn'

When I change "L-BFGS-B" to "BFGS" or anything else that
optim()
support, there's no any error message.

What's wrong?

Any suggestions will be appreciated.

Jinsong

On Tue, 24 Jun 2008, Jinsong Zhao wrote:
> Hi,
>
> When I run the following code,
>
> r <- c(3,4,4,3,5,4,5,9,8,11,12,13)
> n <- rep(15,12)
> x <- c(0, 1.1, 1.3, 2.0, 2.2, 2.8, 3.7, 3.9, 4.4, 4.8, 5.9, 6.8)
> x <- log10(x)
So x[1] = -Inf
> fr <- function(c, alpha, beta) {
>  P <- c + (1-c) * pnorm(alpha + beta * x)
>  P <- pmax(pmin(P,1),0)
>  -(sum(log(choose(n,r))) + sum(r * log(P)) + sum((n -r)* log(1-P)))
> }
> fit <- mle((fr), start = list(c =0.2, alpha = 0, beta =0.1), method >
"L-BFGS-B")
>
> I got:
>
> Error in optim(start, f, method = method, hessian = TRUE, ...) :
>  L-BFGS-B needs finite values of 'fn'
>
> When I change "L-BFGS-B" to "BFGS" or anything else
that optim()
> support, there's no any error message.
>
> What's wrong?
What the message says.  Your fr() function is returning non-finite values. 
Look at what happens if P is 0 or 1 for example.  You seem to need 0 < c <
1, and you have not imposed that constraint.  Adding a trace tp print 
c(c, alpha, beta) I got

[1] 0.2 0.0 0.1
[1] 0.201 0.000 0.100
[1] 0.199 0.000 0.100
[1] 0.200 0.001 0.100
[1]  0.200 -0.001  0.100
[1] 0.200 0.000 0.101
[1] 0.200 0.000 0.099
[1] -0.63692943 -0.54093515  0.01670358
Error in optim(start, f, method = method, hessian = TRUE, ...) :
   L-BFGS-B needs finite values of 'fn'

and other methods also give negative 'c' (and backtrack).  There are
also
problems if beta < 0 (as c + (1-c) * pnorm(alpha + beta * x) = 1 for x = 
-Inf).  So

fit <- mle((fr), start = list(c =0.2, alpha = 0, beta =0.1),
            method = "L-BFGS-B", lower=c(0, -Inf, 0), upper = c(1,
Inf, Inf))

works.

There is no point in using "L-BFGS-B" in a 3-parameter problem unless
you
do impose constraints.

Also, dbinom() will give a more stable way to compute a binomial 
log-likelihood.
> Any suggestions will be appreciated.
>
> Jinsong
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

Another question related with the same code:
> r <- c(3,4,4,3,5,4,5,9,8,11,12,13)
> n <- rep(15,12)
> x <- c(0, 1.1, 1.3, 2.0, 2.2, 2.8, 3.7, 3.9, 4.4, 4.8, 5.9, 6.8)
> x <- log10(x)
> fr <- function(c, alpha, beta) {
>   P <- c + (1-c) * pnorm(alpha + beta * x)
>   P <- pmax(pmin(P,1),0)
>   -(sum(log(choose(n,r))) + sum(r * log(P)) + sum((n -r)* log(1-P)))
>  }
> fit <- mle((fr), start = list(c =0.2, alpha = 0, beta =0.1), method >
"BFGS")
> 
The above code could terminate normally. However, if ``beta'' set to 0, 
then the following error message will appear:

Error in optim(start, f, method = method, hessian = TRUE, ...) :
   initial value in 'vmmin' is not finite

I have learn that "drc" package could deal with this situation, it set
the first observation to control, however, I don't find the way to 
implement it in the source code, for the code is too complex to me.

Any suggestion? Thanks in advance!

Jinsong