nlminb(): NA/NaN function evaluation & Newton failed to find minimum

1,706 views
Skip to first unread message

Ellie

unread,
May 28, 2018, 9:39:53 PM5/28/18
to TMB Users
Dear TMB users

I've been stuck in Na/NaN error message with my TMB code.

I figured out that the error message pops up when log functions get negative values in most cases (i.e., log(-something)) 

I would like to know where and how to check my code with the following error messages:

nlminb(model$par, model$fn, model$gr) : NA/NaN function evaluation

Newton failed to find minimum.

One more question: is this issue can be for wrong initial values?

I appreciate your comments, feedback, advice, ideas in advance.

Attached files are same as below
***sessionInfo***
R version 3.5.0 (2018-04-23)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows >= 8 (build 9200)


attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base    


other attached packages:
[1] TMB_1.7.13


loaded via a namespace (and not attached):
[1] compiler_3.5.0  Matrix_1.2-14   tools_3.5.0     grid_3.5.0      lattice_0.20-35



***CPP file***
#include <TMB.hpp>
template<class Type>
  Type posfun(Type x, Type eps, Type &pen)
{
  if ( x >= eps ){
    return x;
  } else {
    pen += Type(0.01) * pow(x-eps,2);
    return eps/(Type(2.0)-x/eps);
  }
  }

template<class Type>
  Type objective_function<Type>::operator() ()
{
  //data
  DATA_VECTOR(C);
  DATA_VECTOR(I);
  int n = C.size();

  //std::cout << "C "<< C << std::endl;
  //std::cout << "I "<< I << std::endl;

  PARAMETER(logR);
  PARAMETER(logK);
  PARAMETER(logQ);
  PARAMETER(logsdproc);   //log(sd) in the process error;
  PARAMETER(logSigma);
  PARAMETER_VECTOR(P);

  //procedures:

  Type r = exp(logR);
  Type k = exp(logK);
  Type q = exp(logQ);
  Type sdproc = exp(logsdproc);
  Type sigma = exp(logSigma);

  //derived parameters

  vector<Type> Ihat(n);
  vector<Type> fvec1(n);
  vector<Type> tmpP(n); 
  Type f = 0.0;
  Type fpen = 0.0;


for(int t=0; t<(n-1); t++)   { /* from 0 to (n-1) >> size of n*/


      tmpP(t) = P(t) + r*P(t)*(1-P(t))-C(t)/k;
      P(t+1) = posfun(tmpP(t), Type(0.01), fpen);

      f += fpen;

      f -= dnorm(log(P(t+1)), log(tmpP(t)), sdproc, true);
      fvec1(t)=f;
  }

for(int t=0; t<n; t++)  {

  Ihat(t)=q*P(t)*k;
}

f -= sum(dnorm(log(I), log(Ihat), sigma, true));

std::cout << "fvec: "<< fvec1 << std::endl; // posted the result below
std::cout << "temP: "<< tmpP << std::endl;
std::cout << "Ihat: "<< Ihat << std::endl;
std::cout << "r: " << r << std::endl;
std::cout << "K: " << k << std::endl;
std::cout<< "q: "<< q <<std::endl;
std::cout<< "sig_p: "<< sdproc << std::endl;
std::cout<< "sig_o: " << sigma << std::endl;

ADREPORT(P);

  return f;
}


***r code ***
hake <- read.table("schaefer.dat", header=TRUE)
names(hake) <- c("t", "C", "I")
n=c(dim(hake)[1]);  #the number of Ps

fit.sc<-read.csv("C:\\Desktop\\ResSC.csv")
names(fit.sc)<-c("logB","B","K","P(B/K)","P") # for initial P
Pval<-fit.sc$P 
n

length(Pval)
Pval<-round(Pval,1) # for initial values for P
parameters<-list(logR=-1.0, logK=8.0, logQ=-7.75, logsdproc = -5, logSigma = -3, P=Pval)
require(TMB);
compile("schakev4.cpp", "-O1 -g", DLLFLAGS="")
gdbsource("schakev4.R", interactive=TRUE) #debug first!

dyn.load(dynlib("schakev4"))

model<- MakeADFun(hake, parameters, random="P", DLL="schakev4"); 
model$par
fit <- nlminb(model$par, model$fn, model$gr, lower=c(-2, 7, -8, -6, -4, (Pval-0.3)), upper=c(-0.5, 9, -6, -3, -1, Pval))
rep <- sdreport(model)
print(summary(rep))

***Error messages from gdbsource***
    Optimizing tape... Done
iter: 1  value: 36.48854 mgc: 312.3532 ustep: 0.0006262857
iter: 2  value: 36.4721 mgc: 250.0312 ustep: 0.01264024
iter: 3  value: 30.03316 mgc: 460.8754 ustep: 0.1125176
iter: 4  value: 29.44054 mgc: 109.9746 ustep: 0.3355028
iter: 5  value: 29.43263 mgc: 11.68191 ustep: 0.5792682
iter: 6  value: 29.43263 mgc: 0.1750893 ustep: 0.7611206
iter: 7  value: 29.43263 mgc: 4.715644e-05 ustep: 0.872435
Error in newton(par = c(P = 1, P = 1, P = 0.9, P = 0.9, P = 0.8, P = 0.7,  :
  Newton failed to find minimum.
iter: 1  value: 36.48854 mgc: 312.3532 ustep: 0.0006262857
iter: 2  value: 36.4721 mgc: 250.0312 ustep: 0.01264024
iter: 3  value: 30.03316 mgc: 460.8754 ustep: 0.1125176
iter: 4  value: 29.44054 mgc: 109.9746 ustep: 0.3355028
iter: 5  value: 29.43263 mgc: 11.68191 ustep: 0.5792682
iter: 6  value: 29.43263 mgc: 0.1750893 ustep: 0.7611206
iter: 7  value: 29.43263 mgc: 4.715644e-05 ustep: 0.872435
Error in newton(par = c(P = 1, P = 1, P = 0.9, P = 0.9, P = 0.8, P = 0.7,  :
  Newton failed to find minimum.
In addition: Warning message:
In nlminb(model$par, model$fn, model$gr) : NA/NaN function evaluation
outer mgc:  NaN
Error in nlminb(model$par, model$fn, model$gr) :
  gradient function must return a numeric vector of length 5

***DATA: Schaefer ***
year    catch   cpue
1967    15.9    61.89
1968    25.7    78.98
1969    28.5    55.59
1970    23.7    44.61
1971    25.0    56.89
1972    33.3    38.27
1973    28.2    33.84
1974    19.7    36.13
1975    17.5    41.95
1976    19.3    36.63
1977    21.6    36.33
1978    23.1    38.82
1979    22.5    34.32
1980    22.5    37.64
1981    23.6    34.01
1982    29.1    32.16
1983    14.4    26.88
1984    13.2    36.61
1985    28.4    30.07
1986    34.6    30.75
1987    37.5    23.36
1988    25.9    22.36
1989    25.3    21.91

***Pval for initial value***
1.000
0.967
0.903
0.867
0.774
0.725
0.656
0.610
0.483
0.442
0.420
0.401
0.352
0.339
0.331
0.353
0.404
0.461
0.490
0.506
0.518
0.535
0.545
0.589

***vector fvec1 (from gdbsource)***
fvec:
-4.08106
-8.16212
-12.2432
-16.3242
-20.4053
-24.4864
-28.5674
-32.6485
-36.7296
-40.8106
-44.8917
-48.9727
-53.0538
-57.1349
-61.2159
 -65.297
 -69.378
-73.4591
-77.5402
-81.6212
-85.7023
-89.7834
-93.8644
       0
schakev4.r
schakev4.cpp
ResSC.csv
schaefer.dat

Kasper Kristensen

unread,
May 29, 2018, 3:15:52 AM5/29/18
to TMB Users

One thing you may find useful for debugging (without gdb) is to set

obj$env$tracepar <- TRUE

This will print every parameter passed to obj$fn.


If this is not sufficient you may want to do what Brad suggests, i.e. to use CppAD::PrintFor rather than std::cout.
To use PrintFor with TMB a few extra steps are required.

1. Add a dummy version of PrintFor to your cpp file to handle the case Type=double:

namespace CppAD {
  void PrintFor(const char* before, const double& var) {  }
}


2. Disable the tape optimizer on the R side. (otherwise the print statements are optimized out of the tape):

config(optimize.instantly=0, DLL="schakeV4")

You can now use PrintFor from within the objective function like this:

Type a;
CppAD::printFor("a=", a);

Hans Skaug

unread,
May 29, 2018, 8:27:54 AM5/29/18
to TMB Users
An observation on your code:

You have P as the "latent variable", but it is log(P) that has a gaussian distribution:

      f -= dnorm(log(P(t+1)), log(tmpP(t)), sdproc, true);

I recommend using a logP parameterization of the model: random="logP"

Ellie

unread,
Jul 1, 2018, 2:22:31 AM7/1/18
to TMB Users
Dear Mr. Skaug

I've changed my code following your suggestion (i.e., use logP not P), but I've got same warning and error messages which are implying that the optimization failed.
Also, I figured out the warning and error messages are come up dependently with the condition on the first logP value in .cpp (i.e., logP(0)=0.0; at the line #48).
It would be appreciated if you give me any clues or suggestions on this code.

I attached code and data files for reproducible example.

1. under condition logP=0

Error in ev(obj$env$par) : Wrong range component.
In addition: Warning messages:
1: In nlminb(model$par, model$fn, model$gr, lower = c(lower, rep(min(parameters$logP),  :
  NA/NaN function evaluation
2: In he(par) : restarting interrupted promise evaluation
outer mgc:  NaN 
Error in nlminb(model$par, model$fn, model$gr, lower = c(lower, rep(min(parameters$logP),  : 
  gradient function must return a numeric vector of length 4

2. without the condition logP=0

iter: 1  value: 151.7913 mgc: 312.358 ustep: 0.001442378 
iter: 2  value: 83.5787 mgc: 1898.813 ustep: 0.03807486 
iter: 3  value: 62.84069 mgc: 675.6235 ustep: 0.1952083 
iter: 4  value: 59.60101 mgc: 194.6069 ustep: 0.4418796 
iter: 5  value: 59.4731 mgc: 33.23974 ustep: 0.6647738 
iter: 6  value: 59.47282 mgc: 1.516256 ustep: 0.8153551 
iter: 7  value: 59.47282 mgc: 0.003623545 ustep: 0.9029799 
iter: 8  value: 59.47282 mgc: 1.168406e-07 ustep: 0.9502575 
Error in newton(par = c(logP = 0, logP = -0.0337184160322439, logP = -0.101591860412746,  : 
  Newton failed to find minimum.
In addition: Warning message:
In nlminb(model$par, model$fn, model$gr, lower = c(lower, rep(min(parameters$logP),  :
  NA/NaN function evaluation


<cpp file>
#include <TMB.hpp>

template<class Type>
Type posfun(Type x, Type eps, Type &pen)
{
  if ( x >= eps ){
    return x;
  } else {
    pen += Type(0.01)*pow(x-eps,2);
    return eps/(Type(2.0)-x/eps);
  }
};

template<class Type>
Type objective_function<Type>::operator() ()
{
  //data
  DATA_VECTOR(C);
  DATA_VECTOR(I);
  int n = C.size();
  
  std::cout<<"C: "<<C<<std::endl;
  std::cout<<"n: "<<n<<std::endl; 
 
  //free parameters
  PARAMETER(logR);
  PARAMETER(logK);
  PARAMETER(logQ);
  PARAMETER(logSdObs);      //SdObs: SD of the observed error;                 
  PARAMETER_VECTOR(logP);   //skaug's recommendation  //the size = n;
    
  //procedures:
  Type r = exp(logR);
  Type k = exp(logK);
  Type q = exp(logQ);
  Type SdObs = exp(logSdObs);
  Type SdProc = Type(0.5)*SdObs;  //assumption; 
  
  //derived parameters
  vector<Type> P(n);    
  vector<Type> Ihat(n);
  
  Type f;
  Type fpen=0.0;
  Type tmpP;   
  
  //logP(0)=0.0;
  for(int t=0; t<(n-1); t++)   { 
    P(t)=exp(logP(t));
    tmpP = P(t) + r*P(t)*(1-P(t))-C(t)/k;
    P(t+1) = posfun(tmpP, Type(0.01), fpen);
    f+=fpen;
    
    logP(t+1)=log(P(t+1)); 
    f -= dnorm((logP(t+1)+1.e-10), log(tmpP+1.e-10), SdProc, true);
  };
  std::cout << "f: "<< f << std::endl;
  std::cout << "fpen: "<< fpen << std::endl;
  
  for(int t=0; t<n; t++)  {
    Ihat(t)=q*P(t)*k;   //Ihat(t)=q*B(t);
  };
  
  f -= sum(dnorm(log(I), log(Ihat+1.e-10), SdObs, true));
  std::cout << "temP: "<< tmpP << std::endl;
  std::cout << "I "<< I << std::endl;
  std::cout << "Ihat: "<< Ihat << std::endl;
  std::cout << "r: " << r << std::endl;
  std::cout << "K: " << k << std::endl;
  std::cout<< "q: "<< q <<std::endl;
  std::cout<< "SD_p: "<< SdProc << std::endl;
  std::cout<< "SD_o: " << SdObs << std::endl;

  REPORT(logP);        
  //ADREPORT(logP);

  return f;
}



<R file>
####For hake data
getwd()
hake <- read.table("schaefer.dat", header=TRUE)
names(hake) <- c("t", "C", "I")
n=c(dim(hake)[1]);  #the number of Bs

fit.sc<-read.csv("ResSC.csv") #for desktop
names(fit.sc)<-c("logB","B","K","P(B/K)","P")
Pval<-fit.sc$P

require(TMB)
compile("schakev6_play.cpp", "-O1 -g", DLLFLAGS="")
dyn.load(dynlib("schakev6_play"))
config(optimize.instantly=0, DLL="schakev6_play")
parameters<-list(logR=-1.00001, logK=8.0, logQ=-7.75, logSdObs = -3, logP=log(Pval))
parvec=c()

#to give boundaries
lower=c()
upper=c()
for (i in 1:length(parvec)) {
  lower[i]=log(parvec[i])-1
  upper[i]=log(parvec[i])+1
}
lower
upper

model<- MakeADFun(hake, parameters, random="logP", DLL="schakev6_play")
fit <- nlminb(model$par, model$fn, model$gr, 
              lower=c(lower, rep(min(parameters$logP),n)), 
              upper=c(upper, rep(min(parameters$logP),n))) #for logP

rep <- sdreport(model)
print(summary(rep))






2018년 5월 29일 화요일 오후 9시 27분 54초 UTC+9, Hans Skaug 님의 말:
schakev6_play.cpp
schakev6_play.R
schaefer.dat
ResSC.csv
Reply all
Reply to author
Forward
0 new messages