(Module): Add correlation matrix of one-step

probabilities and covariance matrix

diff --git a/src/imach.c b/src/imach.c
index 0c01f952ef1bfe364470cf4e9a0e77cffb65a725..d1c74555e36515da874fc8cac64dd290ffc47f16 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -96,7 +96,7 @@ double jmean; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */\r
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */\r
 FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;\r
-FILE *ficgp,*ficresprob,*ficpop;\r
+FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;\r
 FILE *ficreseij;\r
   char filerese[FILENAMELENGTH];\r
  FILE  *ficresvij;\r
@@ -1641,14 +1641,10 @@ void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, in
          }\r
        }\r
       }\r
-     \r
-   \r
-\r
       for(j=1; j<= nlstate*2; j++)\r
        for(h=0; h<=nhstepm-1; h++){\r
          gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];\r
        }\r
-\r
      } \r
    \r
 /* End theta */\r
@@ -1658,15 +1654,15 @@ void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, in
      for(h=0; h<=nhstepm-1; h++)\r
       for(j=1; j<=nlstate*2;j++)\r
        for(theta=1; theta <=npar; theta++)\r
-       trgradg[h][j][theta]=gradg[h][theta][j];\r
-\r
+         trgradg[h][j][theta]=gradg[h][theta][j];\r
+     \r
 \r
      for(i=1;i<=nlstate*2;i++)\r
       for(j=1;j<=nlstate*2;j++)\r
        varhe[i][j][(int)age] =0.;\r
 \r
      printf("%d|",(int)age);fflush(stdout);\r
-    for(h=0;h<=nhstepm-1;h++){\r
+     for(h=0;h<=nhstepm-1;h++){\r
       for(k=0;k<=nhstepm-1;k++){\r
        matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);\r
        matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);\r
@@ -1931,7 +1927,7 @@ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], dou
 void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)\r
 {\r
   int i, j, i1, k1, j1, z1;\r
-  int k=0, cptcode;\r
+  int k=0,l, cptcode;\r
   double **dnewm,**doldm;\r
   double *xp;\r
   double *gp, *gm;\r
@@ -1939,24 +1935,43 @@ void varprob(char fileres[], double **matcov, double x[], double delti[], int nl
   double age,agelim, cov[NCOVMAX];\r
   int theta;\r
   char fileresprob[FILENAMELENGTH];\r
+  char fileresprobcov[FILENAMELENGTH];\r
+  char fileresprobcor[FILENAMELENGTH];\r
 \r
   strcpy(fileresprob,"prob"); \r
   strcat(fileresprob,fileres);\r
   if((ficresprob=fopen(fileresprob,"w"))==NULL) {\r
     printf("Problem with resultfile: %s\n", fileresprob);\r
   }\r
+  strcpy(fileresprobcov,"probcov"); \r
+  strcat(fileresprobcov,fileres);\r
+  if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {\r
+    printf("Problem with resultfile: %s\n", fileresprobcov);\r
+  }\r
+  strcpy(fileresprobcor,"probcor"); \r
+  strcat(fileresprobcor,fileres);\r
+  if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {\r
+    printf("Problem with resultfile: %s\n", fileresprobcor);\r
+  }\r
   printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);\r
+  printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);\r
+  printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);\r
   \r
-fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");\r
+  fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");\r
   fprintf(ficresprob,"# Age");\r
+  fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");\r
+  fprintf(ficresprobcov,"# Age");\r
+  fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");\r
+  fprintf(ficresprobcov,"# Age");\r
   for(i=1; i<=nlstate;i++)\r
-    for(j=1; j<=(nlstate+ndeath);j++)\r
+    for(j=1; j<=(nlstate+ndeath);j++){\r
       fprintf(ficresprob," p%1d-%1d (SE)",i,j);\r
-\r
-\r
+      fprintf(ficresprobcov," p%1d-%1d ",i,j);\r
+      fprintf(ficresprobcor," p%1d-%1d ",i,j);\r
+    }  \r
   fprintf(ficresprob,"\n");\r
-\r
-\r
+  fprintf(ficresprobcov,"\n");\r
+  fprintf(ficresprobcor,"\n");\r
   xp=vector(1,npar);\r
   dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);\r
   doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));\r
@@ -1971,15 +1986,20 @@ fprintf(ficresprob,"#One-step probabilities and standard deviation in parenthese
 \r
     if  (cptcovn>0) {\r
       fprintf(ficresprob, "\n#********** Variable "); \r
+      fprintf(ficresprobcov, "\n#********** Variable "); \r
+      fprintf(ficresprobcor, "\n#********** Variable "); \r
       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);\r
       fprintf(ficresprob, "**********\n#");\r
+      for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);\r
+      fprintf(ficresprobcov, "**********\n#");\r
+      for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);\r
+      fprintf(ficresprobcor, "**********\n#");\r
     }\r
     \r
       for (age=bage; age<=fage; age ++){ \r
        cov[2]=age;\r
        for (k=1; k<=cptcovn;k++) {\r
          cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];\r
-         \r
        }\r
        for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];\r
        for (k=1; k<=cptcovprod;k++)\r
@@ -2037,16 +2057,33 @@ fprintf(ficresprob,"#One-step probabilities and standard deviation in parenthese
          }\r
        }\r
      \r
-     /*printf("\n%d ",(int)age);\r
+       /*printf("\n%d ",(int)age);\r
      for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){\r
        printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));\r
      }*/\r
 \r
        fprintf(ficresprob,"\n%d ",(int)age);\r
+       fprintf(ficresprobcov,"\n%d ",(int)age);\r
+       fprintf(ficresprobcor,"\n%d ",(int)age);\r
 \r
        for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)\r
-         fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));\r
-  \r
+         fprintf(ficresprob,"%12.3e (%12.3e) ",gm[i],sqrt(doldm[i][j]));\r
+       for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){\r
+         fprintf(ficresprobcov,"%12.3e ",gm[i]);\r
+         fprintf(ficresprobcor,"%12.3e ",gm[i]);\r
+       }\r
+       i=0;\r
+       for (k=1; k<=(nlstate);k++){\r
+         for (l=1; l<=(nlstate+ndeath);l++){ \r
+           i=i++;\r
+           fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);\r
+           fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);\r
+           for (j=1; j<=i;j++){\r
+             fprintf(ficresprobcov," %12.3e",doldm[i][j]);\r
+             fprintf(ficresprobcor," %12.3e",doldm[i][j]/sqrt(doldm[i][i])/sqrt(doldm[j][j]));\r
+           }\r
+         }\r
+       }\r
       }\r
     }\r
     free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));\r
@@ -2056,7 +2093,8 @@ fprintf(ficresprob,"#One-step probabilities and standard deviation in parenthese
   }\r
   free_vector(xp,1,npar);\r
   fclose(ficresprob);\r
-  \r
+  fclose(ficresprobcov);\r
+  fclose(ficresprobcor);\r
 }\r
 \r
 \r
@@ -2099,9 +2137,11 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
  fprintf(fichtm,"\n<li> Result files (second order: variances)<br>\n\r
  - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n\r
  - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n\r
+ - Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n\r
+ - Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%s</a> <br>\n\r
  - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n \r
  - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n\r
- - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);\r
+ - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);\r
 \r
  if(popforecast==1) fprintf(fichtm,"\n\r
  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n\r
@@ -2159,7 +2199,7 @@ void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],ch
   int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;\r
   int ng;\r
   strcpy(optionfilegnuplot,optionfilefiname);\r
-  strcat(optionfilegnuplot,".gp.txt");\r
+  strcat(optionfilegnuplot,".gp");\r
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {\r
     printf("Problem with file %s",optionfilegnuplot);\r
   }\r
@@ -2304,7 +2344,7 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
         for(k=1; k<=(nlstate+ndeath); k++) {\r
           if (k != k2){\r
             if(ng==2)\r
-              fprintf(ficgp," %f*exp(p%d+p%d*x",stepm/YEARM,i,i+1);\r
+              fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);\r
             else\r
               fprintf(ficgp," exp(p%d+p%d*x",i,i+1);\r
             ij=1;\r
@@ -2675,7 +2715,7 @@ int main(int argc, char *argv[])
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;\r
   \r
 \r
-  char version[80]="Imach version 0.8f, May 2002, INED-EUROREVES ";\r
+  char version[80]="Imach version 0.8g, May 2002, INED-EUROREVES ";\r
   char *alph[]={"a","a","b","c","d","e"}, str[4];\r
 \r
 \r