* src/imach.c: 'Not a number' more clear if Hessian not invertable

When the Likelihood can't be maximised in all directions, the Hessian
can't be estimated correctly and can be inverted, or, if it is
inverted, some computed variances are 'negative' and their square
roots produce, in C language, what is called a 'nan' which means 'Not a number'.
A new sentence is printed when the Wald statistics are produced.

diff --git a/src/imach.c b/src/imach.c
index 58b36d84be5eff95b06d56bdb01adc22f00e6d36..ef8910fc9250873e636af2269bcd82d830ed542e 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -6818,8 +6818,8 @@ void hesscov(double **matcov, double **hess, double p[], int npar, double delti[
 
   /* printf("\n#Covariance matrix#\n"); */
   /* fprintf(ficlog,"\n#Covariance matrix#\n"); */
-  /* for (i=1;i<=npar;i++) {  */
-  /*   for (j=1;j<=npar;j++) {  */
+  /* for (i=1;i<=npar;i++) { */
+  /*   for (j=1;j<=npar;j++) { */
   /*     printf("%.6e ",matcov[i][j]); */
   /*     fprintf(ficlog,"%.6e ",matcov[i][j]); */
   /*   } */
@@ -6970,8 +6970,8 @@ double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,
     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 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);
+      printf("Warning: sub Hessian matrix '%d%d' (%s->%s) does not have positive eigen values lc1=%.12e lc2=%.12e \n",thetai,thetaj,label[thetai],label[thetaj],lc1,lc2);
+      fprintf(ficlog, "Warning: sub Hessian matrix '%d%d' (%s->%s) does not have positive eigen values  lc1=%.12e lc2=%.12e \n",thetai,thetaj,label[thetai],label[thetaj],lc1,lc2);
     }
 #endif
   }
@@ -7034,7 +7034,13 @@ double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,
 
 /************** Inverse of matrix **************/
 void ludcmp(double **a, int n, int *indx, double *d) 
-{ 
+{
+/*   Given a matrix a[1..n][1..n], this routine replaces it by the LU (lower upper) decomposition of a rowwise */
+/* permutation of itself. a and n are input. a is output, arranged as in equation  L.U=a; */
+/* indx[1..n] is an output vector that records the row permutation effected by the partial */
+/* pivoting; d is output as ±1 depending on whether the number of row interchanges was even */
+/* or odd, respectively. This routine is used in combination with lubksb to solve linear equations */
+/* or invert a matrix. */
   int i,imax,j,k; 
   double big,dum,sum,temp; 
   double *vv; 
@@ -7047,10 +7053,13 @@ void ludcmp(double **a, int n, int *indx, double *d)
       if ((temp=fabs(a[i][j])) > big) big=temp; 
     if (big == 0.0){
       printf(" Singular Hessian matrix at row %d:\n",i);
+      fprintf(ficlog," Singular Hessian matrix at row %d:\n",i);
       for (j=1;j<=n;j++) {
        printf(" a[%d][%d]=%f,",i,j,a[i][j]);
        fprintf(ficlog," a[%d][%d]=%f,",i,j,a[i][j]);
       }
+      printf(" ERROR Exiting The Likelihood is too flat in some directions \n(probably because of lack of cases transiting) and the Hessian can be inverted)\n");
+      fprintf(ficlog," ERROR Exiting The Likelihood is too flat in some directions \n (probably because of lack of cases transiting) and the Hessian can be inverted\n");
       fflush(ficlog);
       fclose(ficlog);
       nrerror("Singular matrix in routine ludcmp"); 
@@ -16235,9 +16244,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 */
-      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);
+      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);
       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; */