* src/imach.c: improving warning on Hessian inversion and nan

diff --git a/src/imach.c b/src/imach.c
index ef8910fc9250873e636af2269bcd82d830ed542e..d70ee5f15f7340221b213cc843163d4aaf5ecbb5 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1427,7 +1427,7 @@ double gnuplotversion=GNUPLOTVERSION;
 /* $State$ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 /* $State$ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;

-char copyright[]="November 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024";
+char copyright[]="December 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024";

 char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];


@@ -6785,6 +6785,18 @@ void hesscov(double **matcov, double **hess, double p[], int npar, double delti[
   printf("\n");
   fprintf(ficlog,"\n");
 
   printf("\n");
   fprintf(ficlog,"\n");
 

+ 
+  printf("\n#Hessian matrix#\n");
+  fprintf(ficlog,"\n#Hessian matrix#\n");
+  for (i=1;i<=npar;i++) { 
+    for (j=1;j<=npar;j++) { 
+      printf("%.6e ",hess[i][j]);
+      fprintf(ficlog,"%.6e ",hess[i][j]);
+    }
+    printf("\n");
+    fprintf(ficlog,"\n");
+  }
+

   printf("\nInverting the hessian to get the covariance matrix. Wait...\n");
   fprintf(ficlog,"\nInverting the hessian to get the covariance matrix. Wait...\n");
   
   printf("\nInverting the hessian to get the covariance matrix. Wait...\n");
   fprintf(ficlog,"\nInverting the hessian to get the covariance matrix. Wait...\n");
   


@@ -6805,17 +6817,6 @@ void hesscov(double **matcov, double **hess, double p[], int npar, double delti[
     }
   }
 
     }
   }
 

-  printf("\n#Hessian matrix#\n");
-  fprintf(ficlog,"\n#Hessian matrix#\n");
-  for (i=1;i<=npar;i++) { 
-    for (j=1;j<=npar;j++) { 
-      printf("%.6e ",hess[i][j]);
-      fprintf(ficlog,"%.6e ",hess[i][j]);
-    }
-    printf("\n");
-    fprintf(ficlog,"\n");
-  }
-

   /* printf("\n#Covariance matrix#\n"); */
   /* fprintf(ficlog,"\n#Covariance matrix#\n"); */
   /* for (i=1;i<=npar;i++) { */
   /* printf("\n#Covariance matrix#\n"); */
   /* fprintf(ficlog,"\n#Covariance matrix#\n"); */
   /* for (i=1;i<=npar;i++) { */


@@ -6829,7 +6830,7 @@ void hesscov(double **matcov, double **hess, double p[], int npar, double delti[
 
   /* Recompute Inverse */
   /* for (i=1;i<=npar;i++) */
 
   /* Recompute Inverse */
   /* for (i=1;i<=npar;i++) */

-  /*   for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; */
+  /*   for (j=1;j<=npar;j+<+) a[i][j]=matcov[i][j]; */

   /* ludcmp(a,npar,indx,&pd); */
 
   /*  printf("\n#Hessian matrix recomputed#\n"); */
   /* ludcmp(a,npar,indx,&pd); */
 
   /*  printf("\n#Hessian matrix recomputed#\n"); */


@@ -16244,9 +16245,9 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
 
     if(mle != 0){
       /* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */
 
     if(mle != 0){
       /* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */

-      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.\nPlease notice that if the LogLikelihood is not maximised in all directions, the Hessian should not be inverted; but sometimes it can be successfully inverted but producing wrong results like negative variances. In those case the square roots of a negative number produce 'nan' which means 'not a number'.\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.\nPlease notice that if the LogLikelihood is not maximised in all directions, the Hessian should not be inverted; but sometimes it can be successfully inverted but producing wrong results like negative variances. In those case the square roots of a negative number produce 'nan' which means 'not a number'.\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>Please notice that if the LogLikelihood is not maximised in all directions, the Hessian should not be inverted; but sometimes it can be successfully inverted but producing wrong results like negative variances. In those case the square roots of a negative number produce 'nan' which means 'not a number'.\n</br>",optionfilehtmcov);
+      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.\nPlease notice that if the LogLikelihood is not maximised in all directions, the Hessian might not be invertable; but sometimes it can be successfully inverted but producing wrong results like negative variances. In those case the square roots of a negative number produce 'nan' which means 'not a number' on Unix or -1.#IND on Windows.\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.\nPlease notice that if the LogLikelihood is not maximised in all directions, the Hessian might not be invertable; but sometimes it can be successfully inverted but producing wrong results like negative variances. In those case the square roots of a negative number produce 'nan' which means 'not a number' on Unix or -1.#IND on Windows.\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>Please notice that if the LogLikelihood is not maximised in all directions, the Hessian might not be invertable; but sometimes it can be successfully inverted but producing wrong results like negative variances. In those case the square roots of a negative number produce 'nan' which means 'not a number' on Unix or -1.#IND on Windows.\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; */
       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; */