Hi users,<br><br>I played with ADMB to test it before my real application. I first simulated data resulting from the model I would like to test. It is a weird non-linear model with 24 parameters (so 25 with sigma). JFY, the general expression for this model (in R): Y[i]<-A[i]+(w0+w1*X+w2*X^2)*exp(-(s0+s1*X+s2*X^2)*(B[i]-443))+a1[i]*V^a2[i]+b*Z*(B[i]/443)^n<br>
(don't worry it's an ocean optic model...)<br>"i" is the wavelength of light and there are 8 wavelengths<br>A[i], B[i] are constant depending on i<br>X, V and Z are my independent variables<br>and w0, w1, w2, s0, s1, s2, a1[i], a2[i], b, n are the parameters<br>
<br>I fixed the parameters and simulate Y to test the ADMB script. The retrieved parameters have exactly (and amazingly!) the same values as I fixed (which is not the case when using optim function in R or fminsearch in matlab and it's much longer).<br>
The problem is that I have the "error" message: "Hessian does not appear to be positive definite" and I cannot get a "sdreport" for the parameters or the variable Y.<br><br>I copy paste the last part of the ADMB run as well as the tpl file hereafter.<br>
<br>Thanks for your help,<br><br>Sylvain<br><br>25 variables; iteration 740; function evaluation 907<br>Function value -2.7878808e+05; maximum gradient component mag 5.7352e+18<br>Var Value Gradient |Var Value Gradient |Var Value Gradient <br>
1-36.23879 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460 2.70468e+16<br> 4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288 1.15859e+18<br> 7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865 1.04332e+15<br>
10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213 -2.80633e+16<br> 13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428 7.12243e+15<br> 16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483 3.84233e+13<br>
19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250 1.08502e+15<br> 22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177 2.82176e+09<br> 25 227.8695 -9.76063e+06 |<br> ic > imax in fminim is answer attained ?<br>
Var Value Gradient |Var Value Gradient |Var Value Gradient <br> 1-36.23882 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460 2.70468e+16<br> 4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288 1.15859e+18<br>
7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865 1.04332e+15<br> 10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213 -2.80633e+16<br> 13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428 7.12243e+15<br>
16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483 3.84233e+13<br> 19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250 1.08502e+15<br> 22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177 2.82176e+09<br>
25 227.8695 -9.76063e+06 |<br>Function minimizer not making progress ... is minimum attained?<br>Minimprove criterion = 0.0000e+00<br><br> - final statistics:<br>25 variables; iteration 741; function evaluation 937<br>
Function value -2.7879e+05; maximum gradient component mag 5.7352e+18<br>Exit code = 1; converg criter 1.0000e-12<br>Var Value Gradient |Var Value Gradient |Var Value Gradient <br> 1-36.23882 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460 2.70468e+16<br>
4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288 1.15859e+18<br> 7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865 1.04332e+15<br> 10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213 -2.80633e+16<br>
13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428 7.12243e+15<br> 16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483 3.84233e+13<br> 19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250 1.08502e+15<br>
22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177 2.82176e+09<br> 25 227.8695 -9.76063e+06 |<br>Estimating row 1 out of 25 for hessian<br>Estimating row 2 out of 25 for hessian<br>Estimating row 3 out of 25 for hessian<br>
Estimating row 4 out of 25 for hessian<br>Estimating row 5 out of 25 for hessian<br>Estimating row 6 out of 25 for hessian<br>Estimating row 7 out of 25 for hessian<br>Estimating row 8 out of 25 for hessian<br>Estimating row 9 out of 25 for hessian<br>
Estimating row 10 out of 25 for hessian<br>Estimating row 11 out of 25 for hessian<br>Estimating row 12 out of 25 for hessian<br>Estimating row 13 out of 25 for hessian<br>Estimating row 14 out of 25 for hessian<br>Estimating row 15 out of 25 for hessian<br>
Estimating row 16 out of 25 for hessian<br>Estimating row 17 out of 25 for hessian<br>Estimating row 18 out of 25 for hessian<br>Estimating row 19 out of 25 for hessian<br>Estimating row 20 out of 25 for hessian<br>Estimating row 21 out of 25 for hessian<br>
Estimating row 22 out of 25 for hessian<br>Estimating row 23 out of 25 for hessian<br>Estimating row 24 out of 25 for hessian<br>Estimating row 25 out of 25 for hessian<br>Warning -- Hessian does not appear to be positive definite<br>
Hessian does not appear to be positive definite<br><br>///////////////////////////////////////////////////////////////////////////////////////////////////////////<br>TPLfile<br>///////////////////////////////////////////////////////////////////////////////////////////////////////////<br>
DATA_SECTION<br> init_int Nobs_tot //Total number of observations<br> init_vector Kw(1,8) //Value of Water absorption for each wavelength<br> init_matrix Data(1,Nobs_tot,1,6) // Data arranged in column such as :Kd, Chl, Acdm, Bbp, lambda, convert_lambda<br>
vector Kd(1,Nobs_tot) // Kd data as a vector<br> vector Chl(1,Nobs_tot) // Chlorophyll data as a vector<br> vector Acdm(1,Nobs_tot) // Acdm data as a vector<br> vector Bbp(1,Nobs_tot) // Bbp data as a vector<br>
vector lambda(1,Nobs_tot) // Lambda as vector<br> vector convert_lambda(1,Nobs_tot) // Lambda as an indice 1=320, 2=340, 3=380, 4=412, 5=443, 6=490, 7=510, 8=555<br><br>PARAMETER_SECTION<br> init_number logSigma;<br>
number sigma;<br> init_vector w(1,3) //3 common parameters for the quadratic function of acdm (equivalent to CC)<br> init_vector s(1,3) //3 common parameters for the quadratic function of acdm in the exponent term (equivalent to S)<br>
init_bounded_vector loga1(1,8,-7,-1) //8 wavelength-specific parameters that HAVE to be positive<br> init_bounded_vector loga2(1,8,-1,-0.0001) //8 wavelength-specific parameters (as an exponent) that Have to be positive<br>
init_bounded_number logbeta0(-20,-10) //1 common parameter for the Bbp that HAS to be positive<br> init_bounded_number logeta(-0.7,1.1) //1 common parameter for the power of wavelength ratio that HAS to be positive<br>
vector pred_Kd(1,Nobs_tot) //one vector of the predicted value of Kd<br> vector a1(1,8)<br> sdreport_vector a3(1,8)<br> vector a2(1,8)<br> number beta0<br> number eta<br> objective_function_value f<br><br>
PRELIMINARY_CALCS_SECTION<br> Kd=column(Data,1);<br> Chl=column(Data,2);<br> Acdm=column(Data,3);<br> Bbp=column(Data,4);<br> lambda=column(Data,5); <br> convert_lambda=column(Data,6);<br><br>PROCEDURE_SECTION<br> f=0;<br>
sigma=exp(logSigma);<br> a1=exp(loga1);<br> a2=exp(loga2);<br> a3=a1;<br> beta0=exp(logbeta0);<br> eta=-exp(logeta);<br> for (int i=1;i<=Nobs_tot;i++)<br> {<br> pred_Kd(i)=Kw(convert_lambda(i))+(w(1)+w(2)*Acdm(i)+w(3)*Acdm(i)*Acdm(i))*exp(-(s(1)+s(2)*Acdm(i)+s(3)*Acdm(i)*Acdm(i))*(lambda(i)-443.0))+a1(convert_lambda(i))*pow(Chl(i),a2(convert_lambda(i)))+beta0*Bbp(i)*pow((lambda(i)/443.0),eta);<br>
}<br>// cout<<sum(square(pred_Kd-Kd))<<endl;<br> f=0.5*(Nobs_tot*log(2*M_PI*square(sigma))+sum(square(pred_Kd-Kd)/square(sigma)));<br><br>RUNTIME_SECTION<br> maximum_function_evaluations 200000000<br> convergence_criteria 1.e-12<br>
<br>TOP_OF_MAIN_SECTION<br> arrmblsize = 100000000;<br> gradient_structure::set_MAX_NVAR_OFFSET(3000);<br><br>