Version 0.8a. Standard deviation of health expectancies eij.

diff --git a/src/imach.c b/src/imach.c
index bc64e0b530c52b88bf193e7db008823abd8b3209..0b94cf682b0dcf16b8a93bf4c7eec8a44ed75231 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1352,8 +1352,11 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
            if ((k2>=dateprev1) && (k2<=dateprev2)) {\r
              if(agev[m][i]==0) agev[m][i]=agemax+1;\r
              if(agev[m][i]==1) agev[m][i]=agemax+2;\r
-             if (m<lastpass) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];\r
-             /* freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];  */ \r
+             if (m<lastpass) \r
+               if (calagedate>0) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];\r
+             else\r
+              freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];\r
+              freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i]; \r
            }\r
          }\r
        }\r
@@ -1386,7 +1389,7 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
        }\r
     }\r
   }\r
-  \r
+\r
   \r
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);\r
   free_vector(pp,1,nlstate);\r
@@ -1531,18 +1534,29 @@ void tricode(int *Tvar, int **nbcode, int imx)
 \r
 /*********** Health Expectancies ****************/\r
 \r
-void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm)\r
+void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov )\r
+\r
 {\r
   /* Health expectancies */\r
-  int i, j, nhstepm, hstepm, h, nstepm;\r
+  int i, j, nhstepm, hstepm, h, nstepm, k, cptj;\r
   double age, agelim, hf;\r
-  double ***p3mat;\r
+  double ***p3mat,***varhe;\r
+  double **dnewm,**doldm;\r
+  double *xp;\r
+  double **gp, **gm;\r
+  double ***gradg, ***trgradg;\r
+  int theta;\r
+\r
+  varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage);\r
+  xp=vector(1,npar);\r
+  dnewm=matrix(1,nlstate*2,1,npar);\r
+  doldm=matrix(1,nlstate*2,1,nlstate*2);\r
   \r
   fprintf(ficreseij,"# Health expectancies\n");\r
   fprintf(ficreseij,"# Age");\r
   for(i=1; i<=nlstate;i++)\r
     for(j=1; j<=nlstate;j++)\r
-      fprintf(ficreseij," %1d-%1d",i,j);\r
+      fprintf(ficreseij," %1d-%1d (SE)",i,j);\r
   fprintf(ficreseij,"\n");\r
 \r
   if(estepm < stepm){\r
@@ -1582,27 +1596,113 @@ void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, in
     /* if (stepm >= YEARM) hstepm=1;*/\r
     nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */\r
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+    gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2);\r
+    gp=matrix(0,nhstepm,1,nlstate*2);\r
+    gm=matrix(0,nhstepm,1,nlstate*2);\r
+\r
     /* Computed by stepm unit matrices, product of hstepm matrices, stored\r
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */\r
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);  \r
-  \r
-    /*for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++) printf("%f %.5f\n", age*12+h, p3mat[1][1][h]);*/\r
+ \r
 \r
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */\r
+\r
+    /* Computing Variances of health expectancies */\r
+\r
+     for(theta=1; theta <=npar; theta++){\r
+      for(i=1; i<=npar; i++){ \r
+       xp[i] = x[i] + (i==theta ?delti[theta]:0);\r
+      }\r
+      hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  \r
+  \r
+      cptj=0;\r
+      for(j=1; j<= nlstate; j++){\r
+       for(i=1; i<=nlstate; i++){\r
+         cptj=cptj+1;\r
+         for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){\r
+           gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;\r
+         }\r
+       }\r
+      }\r
+     \r
+     \r
+      for(i=1; i<=npar; i++) \r
+       xp[i] = x[i] - (i==theta ?delti[theta]:0);\r
+      hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  \r
+      \r
+      cptj=0;\r
+      for(j=1; j<= nlstate; j++){\r
+       for(i=1;i<=nlstate;i++){\r
+         cptj=cptj+1;\r
+         for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){\r
+           gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;\r
+         }\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
+\r
+     trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar);\r
+\r
+     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
+\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
+    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
+       for(i=1;i<=nlstate*2;i++)\r
+         for(j=1;j<=nlstate*2;j++)\r
+           varhe[i][j][(int)age] += doldm[i][j]*hf*hf;\r
+      }\r
+    }\r
+\r
+     \r
+    /* Computing expectancies */\r
     for(i=1; i<=nlstate;i++)\r
       for(j=1; j<=nlstate;j++)\r
        for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){\r
          eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;\r
-         /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/\r
+         \r
+/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/\r
+\r
        }\r
+\r
     fprintf(ficreseij,"%3.0f",age );\r
+    cptj=0;\r
     for(i=1; i<=nlstate;i++)\r
       for(j=1; j<=nlstate;j++){\r
-       fprintf(ficreseij," %9.4f", eij[i][j][(int)age]);\r
+       cptj++;\r
+       fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );\r
       }\r
     fprintf(ficreseij,"\n");\r
+   \r
+    free_matrix(gm,0,nhstepm,1,nlstate*2);\r
+    free_matrix(gp,0,nhstepm,1,nlstate*2);\r
+    free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2);\r
+    free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);\r
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
   }\r
+  free_vector(xp,1,npar);\r
+  free_matrix(dnewm,1,nlstate*2,1,npar);\r
+  free_matrix(doldm,1,nlstate*2,1,nlstate*2);\r
+  free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage);\r
 }\r
 \r
 /************ Variance ******************/\r
@@ -1842,8 +1942,17 @@ void varprob(char fileres[], double **matcov, double x[], double delti[], int nl
   if((ficresprob=fopen(fileresprob,"w"))==NULL) {\r
     printf("Problem with resultfile: %s\n", fileresprob);\r
   }\r
-  printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);\r
+  printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);\r
   \r
+fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");\r
+  fprintf(ficresprob,"# Age");\r
+  for(i=1; i<=nlstate;i++)\r
+    for(j=1; j<=(nlstate+ndeath);j++)\r
+      fprintf(ficresprob," p%1d-%1d (SE)",i,j);\r
+\r
+\r
+  fprintf(ficresprob,"\n");\r
+\r
 \r
   xp=vector(1,npar);\r
   dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);\r
@@ -1933,7 +2042,7 @@ void varprob(char fileres[], double **matcov, double x[], double delti[], int nl
        fprintf(ficresprob,"\n%d ",(int)age);\r
 \r
        for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)\r
-         fprintf(ficresprob,"%.3e (%.3e) ",gm[i],doldm[i][i]);\r
+         fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));\r
   \r
       }\r
     }\r
@@ -1979,9 +2088,10 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
 \r
  fprintf(fichtm,"\n\r
  - Parameter file with estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\n\r
- - Variances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n\r
+  - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n\r
+ - Variances 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: <a href=\"t%s\">t%s</a> <br>\n\r
- - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);\r
+ - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,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
@@ -2111,10 +2221,19 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),c
 \r
   for (k1=1; k1<= m ; k1 ++) { \r
     for (cpt=1; cpt<= nlstate ; cpt ++) {\r
-      k=2+nlstate*(cpt-1);\r
+      k=2+nlstate*(2*cpt-2);\r
       fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);\r
+      /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);\r
+ for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");\r
+fprintf(ficgp,"\" t \"e%d1\" w l",cpt);\r
+fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);\r
+ for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");\r
+fprintf(ficgp,"\" t \"e%d1\" w l",cpt);\r
+\r
+*/\r
       for (i=1; i< nlstate ; i ++) {\r
-       fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);\r
+       fprintf(ficgp," ,\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+2*i,cpt,i+1);\r
+\r
       } \r
       fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);\r
     }\r
@@ -2514,7 +2633,7 @@ int main(int argc, char *argv[])
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; \r
   int mobilav=0,popforecast=0;\r
   int hstepm, nhstepm;\r
-  double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1;\r
+  double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate;\r
 \r
   double bage, fage, age, agelim, agebase;\r
   double ftolpl=FTOL;\r
@@ -2532,7 +2651,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.8b, March 2002, INED-EUROREVES ";\r
+  char version[80]="Imach version 0.8a, May 2002, INED-EUROREVES ";\r
   char *alph[]={"a","a","b","c","d","e"}, str[4];\r
 \r
 \r
@@ -3249,9 +3368,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   if((stepm == 1) && (strcmp(model,".")==0)){\r
     prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);\r
     if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);\r
-    free_matrix(mint,1,maxwav,1,n);\r
-    free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);\r
-    free_vector(weight,1,n);}\r
+  } \r
   else{\r
     erreur=108;\r
     printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);\r
@@ -3281,6 +3398,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){
     printf("Problem with variance resultfile: %s\n", fileresv);exit(0);\r
   }\r
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);\r
+  calagedate=-1;\r
+prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);\r
 \r
   k=0;\r
   for(cptcov=1;cptcov<=i1;cptcov++){\r
@@ -3303,7 +3422,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){
 \r
       eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);\r
       oldm=oldms;savm=savms;\r
-      evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm);  \r
+      evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);  \r
+ \r
       vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);\r
       oldm=oldms;savm=savms;\r
        varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);\r
@@ -3326,9 +3446,11 @@ while((c=getc(ficpar))=='#' && c!= EOF){
        for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){\r
          for(i=1, epj[j]=0.;i <=nlstate;i++) {\r
            epj[j] += prlim[i][i]*eij[i][j][(int)age];\r
+           /*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/\r
          }\r
          epj[nlstate+1] +=epj[j];\r
        }\r
+\r
        for(i=1, vepp=0.;i <=nlstate;i++)\r
          for(j=1;j <=nlstate;j++)\r
            vepp += vareij[i][j][(int)age];\r
@@ -3340,7 +3462,9 @@ while((c=getc(ficpar))=='#' && c!= EOF){
       }\r
     }\r
   }\r
-\r
+free_matrix(mint,1,maxwav,1,n);\r
+    free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);\r
+    free_vector(weight,1,n);\r
   fclose(ficreseij);\r
   fclose(ficresvij);\r
   fclose(ficrest);\r