--- imach096d/src/imach.c 2002/04/05 15:45:00 1.39 +++ imach096d/src/imach.c 2002/05/21 18:44:41 1.42 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.39 2002/04/05 15:45:00 lievre Exp $ +/* $Id: imach.c,v 1.42 2002/05/21 18:44:41 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -56,7 +56,7 @@ #include #define MAXLINE 256 -#define GNUPLOTPROGRAM "wgnuplot" +#define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ #define FILENAMELENGTH 80 /*#define DEBUG*/ @@ -1330,10 +1330,10 @@ void prevalence(int agemin, float agemax j=cptcoveff; if (cptcovn<1) {j=1;ncodemax[1]=1;} - for(k1=1; k1<=j;k1++){ + for(k1=1; k1<=j;k1++){ for(i1=1; i1<=ncodemax[k1];i1++){ j1++; - + for (i=-1; i<=nlstate+ndeath; i++) for (jk=-1; jk<=nlstate+ndeath; jk++) for(m=agemin; m <= agemax+3; m++) @@ -1352,41 +1352,46 @@ void prevalence(int agemin, float agemax if ((k2>=dateprev1) && (k2<=dateprev2)) { if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; - if (m0) + freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i]; + else + freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; + freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i]; + } } } } } - for(i=(int)agemin; i <= (int)agemax+3; i++){ - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][i]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) + for(i=(int)agemin; i <= (int)agemax+3; i++){ + for(jk=1; jk <=nlstate ; jk++){ + for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) + pp[jk] += freq[jk][m][i]; + } + for(jk=1; jk <=nlstate ; jk++){ + for(m=-1, pos=0; m <=0 ; m++) pos += freq[jk][m][i]; } - for(jk=1; jk <=nlstate ; jk++){ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) - pp[jk] += freq[jk][m][i]; - } - - for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk]; - - for(jk=1; jk <=nlstate ; jk++){ - if( i <= (int) agemax){ - if(pos>=1.e-5){ - probs[i][jk][j1]= pp[jk]/pos; - } - } - } - + for(jk=1; jk <=nlstate ; jk++){ + for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) + pp[jk] += freq[jk][m][i]; } + + for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk]; + + for(jk=1; jk <=nlstate ; jk++){ + if( i <= (int) agemax){ + if(pos>=1.e-5){ + probs[i][jk][j1]= pp[jk]/pos; + } + } + } + + } } } - + free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); @@ -1531,18 +1536,29 @@ void tricode(int *Tvar, int **nbcode, in /*********** Health Expectancies ****************/ -void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm) +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 ) + { /* Health expectancies */ - int i, j, nhstepm, hstepm, h, nstepm; + int i, j, nhstepm, hstepm, h, nstepm, k, cptj; double age, agelim, hf; - double ***p3mat; + double ***p3mat,***varhe; + double **dnewm,**doldm; + double *xp; + double **gp, **gm; + double ***gradg, ***trgradg; + int theta; + + varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage); + xp=vector(1,npar); + dnewm=matrix(1,nlstate*2,1,npar); + doldm=matrix(1,nlstate*2,1,nlstate*2); fprintf(ficreseij,"# Health expectancies\n"); fprintf(ficreseij,"# Age"); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) - fprintf(ficreseij," %1d-%1d",i,j); + fprintf(ficreseij," %1d-%1d (SE)",i,j); fprintf(ficreseij,"\n"); if(estepm < stepm){ @@ -1582,27 +1598,114 @@ void evsij(char fileres[], double ***eij /* if (stepm >= YEARM) hstepm=1;*/ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2); + gp=matrix(0,nhstepm,1,nlstate*2); + gm=matrix(0,nhstepm,1,nlstate*2); + /* Computed by stepm unit matrices, product of hstepm matrices, stored in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij); - - /*for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++) printf("%f %.5f\n", age*12+h, p3mat[1][1][h]);*/ + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ + + /* Computing Variances of health expectancies */ + + for(theta=1; theta <=npar; theta++){ + for(i=1; i<=npar; i++){ + xp[i] = x[i] + (i==theta ?delti[theta]:0); + } + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + + cptj=0; + for(j=1; j<= nlstate; j++){ + for(i=1; i<=nlstate; i++){ + cptj=cptj+1; + for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){ + gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; + } + } + } + + + for(i=1; i<=npar; i++) + xp[i] = x[i] - (i==theta ?delti[theta]:0); + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + + cptj=0; + for(j=1; j<= nlstate; j++){ + for(i=1;i<=nlstate;i++){ + cptj=cptj+1; + for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){ + gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; + } + } + } + + + + for(j=1; j<= nlstate*2; j++) + for(h=0; h<=nhstepm-1; h++){ + gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; + } + + } + +/* End theta */ + + trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar); + + for(h=0; h<=nhstepm-1; h++) + for(j=1; j<=nlstate*2;j++) + for(theta=1; theta <=npar; theta++) + trgradg[h][j][theta]=gradg[h][theta][j]; + + + for(i=1;i<=nlstate*2;i++) + for(j=1;j<=nlstate*2;j++) + varhe[i][j][(int)age] =0.; + + printf("%d||",(int)age);fflush(stdout); + for(h=0;h<=nhstepm-1;h++){ + for(k=0;k<=nhstepm-1;k++){ + matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]); + for(i=1;i<=nlstate*2;i++) + for(j=1;j<=nlstate*2;j++) + varhe[i][j][(int)age] += doldm[i][j]*hf*hf; + } + } + + + /* Computing expectancies */ for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; - /* 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]);*/ + +/* 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]);*/ + } + fprintf(ficreseij,"%3.0f",age ); + cptj=0; for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ - fprintf(ficreseij," %9.4f", eij[i][j][(int)age]); + cptj++; + fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) ); } fprintf(ficreseij,"\n"); + + free_matrix(gm,0,nhstepm,1,nlstate*2); + free_matrix(gp,0,nhstepm,1,nlstate*2); + free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2); + free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar); free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } + free_vector(xp,1,npar); + free_matrix(dnewm,1,nlstate*2,1,npar); + free_matrix(doldm,1,nlstate*2,1,nlstate*2); + free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage); } /************ Variance ******************/ @@ -1842,8 +1945,17 @@ void varprob(char fileres[], double **ma if((ficresprob=fopen(fileresprob,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprob); } - printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob); + printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); +fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n"); + fprintf(ficresprob,"# Age"); + for(i=1; i<=nlstate;i++) + for(j=1; j<=(nlstate+ndeath);j++) + fprintf(ficresprob," p%1d-%1d (SE)",i,j); + + + fprintf(ficresprob,"\n"); + xp=vector(1,npar); dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); @@ -1933,7 +2045,7 @@ void varprob(char fileres[], double **ma fprintf(ficresprob,"\n%d ",(int)age); for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++) - fprintf(ficresprob,"%.3e (%.3e) ",gm[i],doldm[i][i]); + fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i])); } } @@ -1979,9 +2091,10 @@ Interval (in months) between two waves: fprintf(fichtm,"\n - Parameter file with estimated parameters and the covariance matrix: %s
\n - - Variances of life expectancies by age and initial health status (estepm=%d months): v%s
\n + - Variance of one-step probabilities: prob%s
\n + - Variances of life expectancies by age and initial health status (estepm=%d months): v%s
\n - Health expectancies with their variances: t%s
\n - - Standard deviation of stationary prevalences: vpl%s
\n",rfileres,rfileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres); + - Standard deviation of stationary prevalences: vpl%s
\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres); if(popforecast==1) fprintf(fichtm,"\n - Prevalences forecasting: f%s
\n @@ -2004,24 +2117,24 @@ fprintf(fichtm,"
  • Graphs
  • "); fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); fprintf(fichtm," ************\n


    "); } - fprintf(fichtm,"
    - Probabilities: pe%s%d.gif
    -",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); + fprintf(fichtm,"
    - Probabilities: pe%s%d.png
    +",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); for(cpt=1; cpt- Prevalence of disability : p%s%d%d.gif
    -",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); + fprintf(fichtm,"
    - Prevalence of disability : p%s%d%d.png
    +",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"
    - Observed and stationary prevalence (with confident -interval) in state (%d): v%s%d%d.gif
    -",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); +interval) in state (%d): v%s%d%d.png
    +",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"\n
    - Health life expectancies by age and initial health state (%d): exp%s%d%d.gif
    -",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); + fprintf(fichtm,"\n
    - Health life expectancies by age and initial health state (%d): exp%s%d%d.png
    +",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } fprintf(fichtm,"\n
    - Total life expectancy by age and -health expectancies in states (1) and (2): e%s%d.gif
    -",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); +health expectancies in states (1) and (2): e%s%d.png
    +",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); fprintf(fichtm,"\n"); } } @@ -2049,9 +2162,11 @@ m=pow(2,cptcoveff); for (k1=1; k1<= m ; k1 ++) { #ifdef windows - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1); + fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1); #endif #ifdef unix +fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres); #endif @@ -2071,15 +2186,15 @@ for (i=1; i<= nlstate ; i ++) { } fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1)); #ifdef unix -fprintf(ficgp,"\nset ter gif small size 400,300"); +fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\n"); #endif -fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1); } } /*2 eme*/ for (k1=1; k1<= m ; k1 ++) { - fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage); + fprintf(ficgp,"\nset out \"e%s%d.png\" \n\n",strtok(optionfile, "."),k1); + fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage); for (i=1; i<= nlstate+1 ; i ++) { k=2*i; @@ -2104,27 +2219,36 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0"); else fprintf(ficgp,"\" t\"\" w l 0,"); } - fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1); } /*3eme*/ for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt<= nlstate ; cpt ++) { - k=2+nlstate*(cpt-1); - 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); + k=2+nlstate*(2*cpt-2); + fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1); + fprintf(ficgp,"set ter png small\nset size 0.65,0.65\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); + /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); + for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); +fprintf(ficgp,"\" t \"e%d1\" w l",cpt); +fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); + for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); +fprintf(ficgp,"\" t \"e%d1\" w l",cpt); + +*/ for (i=1; i< nlstate ; i ++) { - 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); + 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); + } - fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1); - } } + } /* CV preval stat */ for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt