Оценка функции 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