Оценка функции NA/NaN и Ньютону не удалось найти минимальные сообщения об ошибках от TMB

C++, R пользователи.

Я застрял в сообщении об ошибке Na/NaN в TMB, несмотря на помощь других пользователей. Я понял, что сообщение об ошибке появляется, когда функции журнала получают отрицательные значения в большинстве случаев. Я хотел бы знать, где и как сначала атаковать в моем коде следующие сообщения об ошибках: NA/NaN function evaluation, Newton failed to find minimum

Еще один вопрос: эта проблема может быть для неправильных начальных значений?

Я ценю ваши комментарии, советы, идеи заранее!

файл cpp

#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 >> size of n+1*/


      tmpP(t) = P(t) + r*P(t)*(1-P(t))-C(t)/k;
      P(t+1) = posfun(tmpP(t), Type(0.01), fpen);
      //B(t+1) = posfun(tmpB,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;
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

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:\\ResSC.csv")
names(fit.sc)<-c("logB","B","K","P(B/K)","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))

Сообщения об ошибках от 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 для начального значения

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

0 ответов

Другие вопросы по тегам