From c3bf2d2728b60c455972466550533cd0b08e1152 Mon Sep 17 00:00:00 2001 From: nbrouard Date: Wed, 25 Dec 2024 12:36:15 +0100 Subject: [PATCH] * 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. --- src/imach.c | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/src/imach.c b/src/imach.c index 58b36d8..ef8910f 100644 --- 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

The Wald test results are output only if the maximimzation of the Likelihood is performed (mle=1)\n
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
",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

The Wald test results are output only if the maximimzation of the Likelihood is performed (mle=1)\n
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
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
",optionfilehtmcov); fprintf(fichtm,"\n"); fprintf(fichtm, "\n"); /* ncovmodel=2+nbocc(model,'+')+1; */ -- 2.43.0
Model=1+ age