I'm trying to fit data in R with least squares and negative log likelihood.
However, I'm finding NLL won't find the global minimum that LS does find:
G <- function(pars, data) {
pars['A'] * (1 - exp(-pars['age']*data$age))^pars['theta']
}
NLL = function(pars, data) {
if (pars['age'] < 0){
return(-Inf)
}
if (pars['age'] > 1){
return(-Inf)
}
if (pars['theta'] > 10){
return(-Inf)
}
if (pars['sd'] < 0){
return(-Inf)
}
result = -sum(dnorm(x = data$agbd - G(pars, data), mean = 0, sd = pars['sd'], log = TRUE), na.rm = TRUE)
ifelse(result == 0, -Inf, result)
}
nls <- function(pars, data) {
result <- sum((G(pars, data) - data$agbd)^2)
ifelse(result == 0, -Inf, result)
}
data_example <- data.frame(age = 1:32, agbd = c(19.0, 0.0, 5.0, 5.5, 14.0, 12.5, 15.0, 21.0, 28.0, 31.0, 26.0, 34.5, 51.5, 43.5, 52.5, 47.0, 42.5, 57.0, 48.0, 44.0, 52.5, 76.0, 59.0,
79.0, 62.0, 77.5, 49.5, 58.5, 77.0, 55.5, 77.0, 89.5))
pars <- c(A = 100, age = 0.5, theta = 1.5)
o_nls <- optim(pars, nls, data=data_example, method = 'BFGS')
o_nls # o_nls$value = 2612.017
#####################
# note that
o_nls <- optim(pars, nls, data=data_example)
o_nls # o_nls$value = 2628.707
# so the default won't do as well as BFGS.
####################
pars <- c(A = 100, age = 0.5, theta = 1.5, sd = 0.5)
o_nll <- optim(pars, NLL, data=data_example)
o_nll
sum((data_example$agbd - G(o_nll$par, data_example))^2) #2638.79
I tried different data with the same functional form, and in some cases, for nonlinear least squares, the default optim method indeed did better than BFGS.
Any idea why NLL is not returning the global minimum in any case I try it? It shouldn't be outperformed by nls here.