[ADMB Users] "Hessian does not appear to be positive definite"
Sylvain Bonhommeau
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
(don't worry it's an ocean optic model...)
"i" is the wavelength of light and there are 8 wavelengths
A[i], B[i] are constant depending on i
X, V and Z are my independent variables
and w0, w1, w2, s0, s1, s2, a1[i], a2[i], b, n are the parameters
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).
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.
I copy paste the last part of the ADMB run as well as the tpl file hereafter.
Thanks for your help,
Sylvain
25 variables; iteration 740; function evaluation 907
Function value -2.7878808e+05; maximum gradient component mag 5.7352e+18
Var Value Gradient |Var Value Gradient |Var Value Gradient
1-36.23879 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460 2.70468e+16
4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288 1.15859e+18
7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865 1.04332e+15
10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213 -2.80633e+16
13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428 7.12243e+15
16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483 3.84233e+13
19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250 1.08502e+15
22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177 2.82176e+09
25 227.8695 -9.76063e+06 |
ic > imax in fminim is answer attained ?
Var Value Gradient |Var Value Gradient |Var Value Gradient
1-36.23882 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460 2.70468e+16
4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288 1.15859e+18
7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865 1.04332e+15
10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213 -2.80633e+16
13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428 7.12243e+15
16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483 3.84233e+13
19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250 1.08502e+15
22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177 2.82176e+09
25 227.8695 -9.76063e+06 |
Function minimizer not making progress ... is minimum attained?
Minimprove criterion = 0.0000e+00
- final statistics:
25 variables; iteration 741; function evaluation 937
Function value -2.7879e+05; maximum gradient component mag 5.7352e+18
Exit code = 1; converg criter 1.0000e-12
Var Value Gradient |Var Value Gradient |Var Value Gradient
1-36.23882 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460 2.70468e+16
4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288 1.15859e+18
7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865 1.04332e+15
10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213 -2.80633e+16
13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428 7.12243e+15
16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483 3.84233e+13
19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250 1.08502e+15
22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177 2.82176e+09
25 227.8695 -9.76063e+06 |
Estimating row 1 out of 25 for hessian
Estimating row 2 out of 25 for hessian
Estimating row 3 out of 25 for hessian
Estimating row 4 out of 25 for hessian
Estimating row 5 out of 25 for hessian
Estimating row 6 out of 25 for hessian
Estimating row 7 out of 25 for hessian
Estimating row 8 out of 25 for hessian
Estimating row 9 out of 25 for hessian
Estimating row 10 out of 25 for hessian
Estimating row 11 out of 25 for hessian
Estimating row 12 out of 25 for hessian
Estimating row 13 out of 25 for hessian
Estimating row 14 out of 25 for hessian
Estimating row 15 out of 25 for hessian
Estimating row 16 out of 25 for hessian
Estimating row 17 out of 25 for hessian
Estimating row 18 out of 25 for hessian
Estimating row 19 out of 25 for hessian
Estimating row 20 out of 25 for hessian
Estimating row 21 out of 25 for hessian
Estimating row 22 out of 25 for hessian
Estimating row 23 out of 25 for hessian
Estimating row 24 out of 25 for hessian
Estimating row 25 out of 25 for hessian
Warning -- Hessian does not appear to be positive definite
Hessian does not appear to be positive definite
///////////////////////////////////////////////////////////////////////////////////////////////////////////
TPLfile
///////////////////////////////////////////////////////////////////////////////////////////////////////////
DATA_SECTION
init_int Nobs_tot //Total number of observations
init_vector Kw(1,8) //Value of Water absorption for each wavelength
init_matrix Data(1,Nobs_tot,1,6) // Data arranged in column such as :Kd, Chl, Acdm, Bbp, lambda, convert_lambda
vector Kd(1,Nobs_tot) // Kd data as a vector
vector Chl(1,Nobs_tot) // Chlorophyll data as a vector
vector Acdm(1,Nobs_tot) // Acdm data as a vector
vector Bbp(1,Nobs_tot) // Bbp data as a vector
vector lambda(1,Nobs_tot) // Lambda as vector
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
PARAMETER_SECTION
init_number logSigma;
number sigma;
init_vector w(1,3) //3 common parameters for the quadratic function of acdm (equivalent to CC)
init_vector s(1,3) //3 common parameters for the quadratic function of acdm in the exponent term (equivalent to S)
init_bounded_vector loga1(1,8,-7,-1) //8 wavelength-specific parameters that HAVE to be positive
init_bounded_vector loga2(1,8,-1,-0.0001) //8 wavelength-specific parameters (as an exponent) that Have to be positive
init_bounded_number logbeta0(-20,-10) //1 common parameter for the Bbp that HAS to be positive
init_bounded_number logeta(-0.7,1.1) //1 common parameter for the power of wavelength ratio that HAS to be positive
vector pred_Kd(1,Nobs_tot) //one vector of the predicted value of Kd
vector a1(1,8)
sdreport_vector a3(1,8)
vector a2(1,8)
number beta0
number eta
objective_function_value f
PRELIMINARY_CALCS_SECTION
Kd=column(Data,1);
Chl=column(Data,2);
Acdm=column(Data,3);
Bbp=column(Data,4);
lambda=column(Data,5);
convert_lambda=column(Data,6);
PROCEDURE_SECTION
f=0;
sigma=exp(logSigma);
a1=exp(loga1);
a2=exp(loga2);
a3=a1;
beta0=exp(logbeta0);
eta=-exp(logeta);
for (int i=1;i<=Nobs_tot;i++)
{
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);
}
// cout<<sum(square(pred_Kd-Kd))<<endl;
f=0.5*(Nobs_tot*log(2*M_PI*square(sigma))+sum(square(pred_Kd-Kd)/square(sigma)));
RUNTIME_SECTION
maximum_function_evaluations 200000000
convergence_criteria 1.e-12
TOP_OF_MAIN_SECTION
arrmblsize = 100000000;
gradient_structure::set_MAX_NVAR_OFFSET(3000);
dave fournier
I think there is a bit lost in translation here.
>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).
Do you mean the you used known parameters to simulate some data and then
fit the data and got the parameters back. Is that the same data set
for which the results you posted here are from? anyway if you send me
the tpl file and every thing needed to run the model I'll take a look at it.
Dave
--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com
_______________________________________________
Users mailing list
Us...@admb-project.org
http://lists.admb-project.org/mailman/listinfo/users
dave fournier
Your problme is that the data are too good so the model fits perfectly
with an estimated sigma of almost 0. I added a tiny bit of noise to the
data and it converged fine.
Kd=column(Data,1);
random_number_generator rng(301);
dvector tmp(1,Nobs_tot);
Kd=elem_prod(Kd,exp(.01*tmp)); // add some lognormal noise to Kd
However with real data you may need to civilize the code a bit.
dave fournier
interest other users. If the data are perfect then the residuals are
amost zero and the estimate for the standard deviation is almost zero so
you get a
0/0
situation. very unlikely that this will occur with real data.
If you are moving on to random effects you should first civilize the
model a bit. In somethng like
exp( a(1) + a(2)*v + a(3)*v^2 )
you should at least center v as in
w = v -mean(v)
mfexp( a(1) + a(2)*w + a(3)*w^2 )
The use of mfexp is to avoid overflow.
also you might fit a(1) first as in
PARAMETER_SECTION
init_number a1 // fit in phase 1
init_vector a(2,3,2) .. fit in phase 2
mfexp( a1 + a(2)*w + a(3)*w^2 )
Thanks guys. It runs well now with both methods you mentioned. I seems
like it is not good to have perfect data! I gonna try the random effect
module now...
Cheers,
Sylvain
2009/4/9 Mark Maunder <mmau...@iattc.org>
I see that Dave found that it was related to having no error in the
simulated data. I guess if you fixed the standard deviation (use a
negative phase) it might work without adding error.
Mark
From: Sylvain Bonhommeau [mailto:sylvain.b...@gmail.com]
Sent: Thursday, April 09, 2009 8:39 AM
To: Mark Maunder
Subject: Re: [ADMB Users] "Hessian does not appear to be positive
definite"
Hi Mark,
Thanks for your answer.
I simulated the data with R (code below). There is no error added to
data for now but that is the aim of today. I specified a pin file and
the initial conditions are not (but close to) the "known" parameters.
Thanks,
Sylvain
Sylvain Bonhommeau
Thanks for your answer. That's what I did or tried to do and that's the results for the dataset with known parameters. Here are the .tpl, .pin and .dat files.
Thanks a bunch,
Sylvain
Sylvain Bonhommeau
I want to test the ability of my procedure to estimate parameters so I simulated data with known parameters. As you mentioned previously, when there is no errors, it is not possible to get the standard deviation. However in the .par file, the value of parameters are almost perfect. The problem is when I add a small noise to my simulated data, i.e. multiplying by a random variable normally distributed with mean=1 and sd=0.01, the estimated parameters are not the one I chose...
I attach the .tpl, .pin and the two .dat file (without and with error) as well as the R script to simulate the data.
(this is a simplified version of the model I sent you before)
Thanks for your help,
Sylvain
I think we should keep this on the list as these are issures whihc may
dave fournier
What you are coming up against here are the limits of estimability.
If you remove the std dev from the model and just fit the least squares
and then
calculate the mean-squared residuals you will get 0.000821022 so the model
fits the data almost perfectly. the fact that they are not the values
you used
to generate the data simply reflects the fact that there is not a lot of
information in the data
i.e. thee are different combinations of the model parameters which fit
the model (almost) equally well.
Dave
> Hi,
>
> I want to test the ability of my procedure to estimate parameters so I
> simulated data with known parameters. As you mentioned previously,
> when there is no errors, it is not possible to get the standard
> deviation. However in the .par file, the value of parameters are
> almost perfect. The problem is when I add a small noise to my
> simulated data, i.e. multiplying by a random variable normally
> distributed with mean=1 and sd=0.01, the estimated parameters are not
> the one I chose...
>
> I attach the .tpl, .pin and the two .dat file (without and with error)
> as well as the R script to simulate the data.
>
> (this is a simplified version of the model I sent you before)
>
>
> Thanks for your help,
>
> Sylvain
>
>
>
> 2009/4/9 dave fournier <ot...@otter-rsch.com
> <mailto:ot...@otter-rsch.com>>
> 2009/4/9 Mark Maunder <mmau...@iattc.org <mailto:mmau...@iattc.org>>
> Us...@admb-project.org <mailto:Us...@admb-project.org>
> http://lists.admb-project.org/mailman/listinfo/users