}
/**** Computes Hessian and covariance matrix ***/
-void hesscov(double **matcov, double **hess, double p[], int npar, double delti[], double ftolhess, double (*func)(double []))
+void hesscov(double **matcov, double **hess, double p[], int npar, double delti[], double ftolhess, char label[][16], double (*func)(double []))
{
double **a,**y,*x,pd;
/* double **hess; */
int *indx;
double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar);
- double hessij(double p[], double **hess, double delti[], int i, int j,double (*func)(double []),int npar);
+ double hessij(double p[], double **hess, double delti[], int i, int j,char label[][16],double (*func)(double []),int npar);
void lubksb(double **a, int npar, int *indx, double b[]) ;
void ludcmp(double **a, int npar, int *indx, double *d) ;
double gompertz(double p[]);
if (j>i) {
printf(".%d-%d",i,j);fflush(stdout);
fprintf(ficlog,".%d-%d",i,j);fflush(ficlog);
- hess[i][j]=hessij(p,hess, delti,i,j,func,npar);
+ hess[i][j]=hessij(p,hess, delti,i,j,label,func,npar);
hess[j][i]=hess[i][j];
/*printf(" %lf ",hess[i][j]);*/
}
-double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,double (*func)(double []),int npar)
+double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,char label[][16],double (*func)(double []),int npar)
{
int i;
- int l=1, lmax=20;
+ int l=1, lmax=10;
double k1,k2,k3,k4,res,fx;
double p2[MAXPARM+1];
int k, kmax=1;
int firstime=0;
fx=func(x);
- for (k=1; k<=kmax; k=k+10) {
+ /* for (k=1; k<=kmax; k=k+10) { /\* Widening around the supposed maximum *\/ */
+ for (l=0; l<=lmax; l++) { /* Widening of a factor 10, 20 around the supposed maximum up to 100 times */
for (i=1;i<=npar;i++) p2[i]=x[i];
+ k=1+10*l;
p2[thetai]=x[thetai]+delti[thetai]*k;
p2[thetaj]=x[thetaj]+delti[thetaj]*k;
k1=func(p2)-fx;
k4=func(p2)-fx;
res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */
if(k1*k2*k3*k4 <0.){
- firstime=1;
- kmax=kmax+10;
- }
- if(kmax >=10 || firstime ==1){
+ /* firstime=1; */
+ /* kmax=kmax+10; /\* Trying to enlarge of a factor 10 on all directions in order to get a maximum on each direction *\/ */
+ /* } */
+ /* if(kmax >=10 || firstime ==1){ */
/* What are the thetai and thetaj? thetai/ncovmodel thetai=(thetai-thetai%ncovmodel)/ncovmodel +thetai%ncovmodel=(line,pos) */
- printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
+ printf("Warning: directions %s->%s, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",label[thetai],label[thetaj], ftol);
fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
if ((lc2 <0) || (lc1 <0) ){
- printf("Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj);
- fprintf(ficlog, "Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj);
- printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
- fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
+ printf("Warning: sub Hessian matrix '%d%d' does not have positive eigen values lc1=%.12e lc2=%.12e \n",thetai,thetaj,lc1,lc2);
+ fprintf(ficlog, "Warning: sub Hessian matrix '%d%d' does not have positive eigen values lc1=%.12e lc2=%.12e \n",thetai,thetaj,lc1,lc2);
}
#endif
}
double *p, *pstart; /* p=param[1][1] pstart is for starting values guessed by freqsummary */
double **matcov; /* Matrix of covariance */
double **hess; /* Hessian matrix */
+
+ char labelVar[100][8] ; /* 100 Variables of max length 8=age*age */
+ char label[100][16]; /* 100 Variables of max length 16=V10*V12*age*age */
+
+
double ***delti3; /* Scale */
double *delti; /* Scale */
double ***eij, ***vareij;
free_ivector(flatdir,1,npar);
#endif /* LINMINORIGINAL*/
#endif /* POWELL */
- hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz);
+ hesscov(matcov, hess, p, NDIM, delti, 1e-4, label, gompertz); /* Careful */
for(i=1; i <=NDIM; i++)
for(j=i+1;j<=NDIM;j++)
fprintf(ficlog,"\nCovariance matrix\n ");
for(i=1; i <=NDIM; i++) {
for(j=1;j<=NDIM;j++){
- printf("%f ",matcov[i][j]);
- fprintf(ficlog,"%f ",matcov[i][j]);
+ printf("%f ",matcov[i][j]);
+ fprintf(ficlog,"%f ",matcov[i][j]);
}
printf("\n "); fprintf(ficlog,"\n ");
}
if(mle != 0){
/* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */
- ftolhess=ftol; /* Usually correct */
- hesscov(matcov, hess, p, npar, delti, ftolhess, func);
printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
fprintf(fichtm, "\n<p>The Wald test results are output only if the maximimzation of the Likelihood is performed (mle=1)\n</br>Parameters, Wald tests and Wald-based confidence intervals\n</br> W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n</br> And Wald-based confidence intervals plus and minus 1.96 * W \n </br> It might be better to visualize the covariance matrix. See the page '<a href=\"%s\">Matrix of variance-covariance of one-step probabilities and its graphs</a>'.\n</br>",optionfilehtmcov);
fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">");
fprintf(fichtm, "\n<tr><th>Model=</th><th>1</th><th>+ age</th>");
+ /* ncovmodel=2+nbocc(model,'+')+1; */
+ /* nforce= (nlstate+ndeath-1)*nlstate; /\* Number of forces ij from state i to j *\/ */
+ /* npar= nforce*ncovmodel; /\* Number of parameters like aij*\/ */
+ /* char labelVar[100][16] ; /\* 100 Variables of max length 16 V10*V12*age*age *\/ */
+ /* char label[100][16]; */
+ /* printf("HELLLOOOOO\n"); */
+ /* labelVar[1]="1"; */
+ jk=1;
+ strcpy(labelVar[jk],"1");
+ /* printf("%s\n",labelVar[1]); */
+ /* printf("HHHH %s HH\n",labelVar[1]); */
+ jk++;
+ strcpy(labelVar[jk],"age");
+ /* printf("HELLLOOOOO22222\n"); */
if(nagesqr==1){
+ jk++;
+ sprintf(labelVar[jk],"age*age");
printf(" + age*age ");
fprintf(ficres," + age*age ");
fprintf(ficlog," + age*age ");
fprintf(fichtm, "<th>+ age*age</th>");
}
+ /* printf("HELLLOOOOO3\n"); */
for(j=1;j <=ncovmodel-2-nagesqr;j++){
if(Typevar[j]==0) {
+ jk++;
+ sprintf(labelVar[jk],"V%d",Tvar[j]);
printf(" + V%d ",Tvar[j]);
+ printf(" + V%d jk=%d ncovmodel=%d j=%d",Tvar[j], jk, ncovmodel, j);
fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);
}else if(Typevar[j]==1) {
+ jk++;
+ sprintf(labelVar[jk],"V%d*age",Tvar[j]);
printf(" + V%d*age ",Tvar[j]);
fprintf(fichtm, "<th>+ V%d*age</th>",Tvar[j]);
}else if(Typevar[j]==2) {
+ jk++;
+ sprintf(labelVar[jk],"V%d*V%d",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
fprintf(fichtm, "<th>+ V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
}else if(Typevar[j]==3) { /* TO VERIFY */
+ jk++;
+ sprintf(labelVar[jk],"V%d*V%d*age</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
fprintf(fichtm, "<th>+ V%d*V%d*age</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
}
}
fprintf(fichtm, "</tr>\n");
-
+ jk=0;
+ for(i=1;i<=nlstate;i++){
+ /* printf("i=%d\n",i); */
+ for(j=1;j<=nlstate+ndeath;j++){
+ if(j !=i ){
+ /* printf("i=%d, j=%d\n",i,j); */
+ for(k=1;k<=ncovmodel;k++){
+ jk++;
+ sprintf(label[jk],"%d%d-%s",i,j,labelVar[k]);
+ /* printf("i=%d,j=%d, k=%d,jk=%d, label[%d]=%s, labelVar[%d]=%s\n",i,j,k,jk,jk,label[jk],k,labelVar[k]); */
+ /* printf("i=%d,j=%d, k=%d,jk=%d,labelVar[%d]=%s, %d%d%s\n",i,j,k,jk,k,labelVar[k],i,j,labelVar[k]); */
+ } /* end k */
+ }
+ } /* end j */
+ } /* end i */
+ /* for(jk=1;jk<=npar;jk++) */
+ /* printf("\n\nATTENTION \n\n %d label=%s\n",jk,label[jk]); */
+
+ ftolhess=ftol; /* Usually correct */
+ hesscov(matcov, hess, p, npar, delti, ftolhess, label, func);
+
for(i=1,jk=1; i <=nlstate; i++){
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
git://henry.ined.fr / .git / commitdiff