mathematical optimization - R optimize multiple parameters

Question

I am using R optim() function to estimate set of parameters which optimize user defined function shown below. But optim() out put is:

Error in optim(pstart, llAgedepfn, method = "L-BFGS-B", upper = up, lower = lo) : L-BFGS-B needs finite values of 'fn'

Please help. The complete script is shown below:

dataM<-cbind(c(1.91,0.29,0.08,0.02,0.01,0.28,0.45,0.36,0.42,0.17,0.16,0.06,0.17,0.17,0.12),
               c(0.27,4.54,0.59,0.05,0.04,0.13,0.48,0.68,0.66,0.18,0.11,0.06,0.08,0.08,0.08),
               c(0.07,0.57,4.48,0.48,0.02,0.05,0.09,0.43,0.78,0.52,0.17,0.10,0.05,0.05,0.14),
               c(0.02,0.04,0.44,4.34,0.36,0.09,0.07,0.11,0.41,0.77,0.43,0.10,0.03,0.04,0.14),
               c(0.01,0.04,0.01,0.36,2.20,0.46,0.19,0.15,0.19,0.34,0.62,0.30,0.09,0.03,0.22),
               c(0.22,0.11,0.05,0.09,0.45,0.91,0.61,0.43,0.37,0.26,0.41,0.63,0.29,0.16,0.15),
               c(0.31,0.35,0.07,0.05,0.16,0.54,0.81,0.59,0.48,0.36,0.33,0.43,0.47,0.26,0.20),
               c(0.22,0.45,0.29,0.08,0.11,0.34,0.53,0.85,0.71,0.39,0.27,0.26,0.26,0.28,0.38),
               c(0.22,0.36,0.44,0.26,0.12,0.24,0.36,0.59,0.91,0.61,0.35,0.28,0.20,0.22,0.29),
               c(0.09,0.10,0.30,0.49,0.22,0.17,0.28,0.33,0.62,0.80,0.52,0.29,0.20,0.11,0.46),
               c(0.10,0.07,0.12,0.32,0.48,0.32,0.30,0.27,0.42,0.61,0.78,0.47,0.33,0.23,0.49),
               c(0.04,0.04,0.06,0.08,0.24,0.53,0.41,0.28,0.36,0.36,0.50,0.67,0.51,0.19,0.47),
               c(0.10,0.05,0.04,0.02,0.07,0.23,0.43,0.26,0.23,0.23,0.33,0.48,0.75,0.51,0.49),
               c(0.05,0.04,0.03,0.05,0.02,0.10,0.19,0.22,0.21,0.10,0.18,0.14,0.40,0.79,0.82),
               c(0.03,0.02,0.03,0.03,0.06,0.04,0.06,0.12,0.11,0.18,0.16,0.14,0.16,0.34,1.26)
)

NormCM <- dataM/eigen(CMWkday)$values[1] #Normalizing the contact mtrix - divide by the largest eigen value

w <- c(495,528,548,603,617,634,720,801,957,937,798,755,795,1016,2469) 

g2 <- c(770,622,726,559,410,547,564,472,399,397,340,308,337,91,84) 

h2 <- c(269,426,556,430,271,284,303,207,194,181,126,106,74,24,23) 

z2 <- h2/g2

g1 <- c(774,527,665,508,459,539,543,492,402,412,365,342,213,146,152) 

h1 <- c(56,31,84,173,103,85,123,70,71,80,55,25,18,12,26) 
z1 <- h1/g1

#### Normal loglikelihood #########

llnormfn <- function(q) {  

  tol <- 1e-9
  final.size.start <- 0.8
  zeta <- rep(final.size.start, nrow(NormCM))
  last.zeta <- rep(0, nrow(NormCM))
  first.run <- T
  current.diff <- tol+1
  loglik <- 0

  while (current.diff > tol) {

    zeta <- 1-exp(-(q*(zeta%*%NormCM)))
    current.diff <- sum(abs(last.zeta-zeta))
    last.zeta <-zeta

  }
  mu <- c(zeta)

  zigma <- z1*(1-z1)/g1 + (z1+mu)*(1-(z1+mu))/g2

  logliknorm <- -sum((((z2-z1)-mu)**2)/2*zigma + 0.5*log(2*pi*zigma))

  return(logliknorm)

} 

pstart <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
up <- c(5,5,5,5,5,5,5,5,5,5,5,5,5,5,5)
lo <- c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1)
estm <- optim(pstart, llnormfn, method = "L-BFGS-B", upper = up, lower = lo )

Answer 1

Your llnormfn doesn't return a finite value for all values of its parameters within the range. For example at the upper limit:

> llnormfn(up)
[1] NaN
Warning message:
In log(2 * pi * zigma) : NaNs produced

Because zigma must be less than zero here.

If you restrict the range a bit you can eventually find a spot where it does work...

> llnormfn(up-2)
[1] NaN
Warning message:
In log(2 * pi * zigma) : NaNs produced
> llnormfn(up-3)
[1] 42.96818

Let's check it works at the lower range:

> llnormfn(lo)
[1] 41.92578

that looks fine. So either you've set that upper limit outside the computationally valid range of your function, or you've got a bug in your llnormfn function, or both, or something else.

If you do run the optimisation with a reduced upper bound you do get convergence:

> estm <- optim(pstart, llnormfn, method = "L-BFGS-B", upper = up-3, lower = lo )
> estm
$par
 [1] 1.9042672 1.0891264 0.9916916 0.6208685 1.2413983 1.4822433 1.1243878
 [8] 1.5224263 1.3686933 1.4876350 1.6231518 2.0000000 2.0000000 2.0000000
[15] 2.0000000

$value
[1] 38.32182

$counts
function gradient 
      23       23 

$convergence
[1] 0

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

Although you might notice some of those parameters are at the upper value (2.0) which is an alarm bell.

Check your function behaves sensibly for its input values - try fixing all-but-one and plotting how llnormfn behaves while varying one. I just had a quick look and the function does not look smooth at all, with lots of discontinuities, so I doubt BFGS is a good method for optimising.

e.g varying the fifth parameter between 0.1 and 2:

> s = seq(0.1,2,len=300)
> ss = sapply(1:length(s),function(i){ll=lo;ll[5]=s[i];llnormfn(ll)})
> plot(s,ss)

gives: