--- imach/src/imach.c 2002/04/19 13:46:19 1.40 +++ imach/src/imach.c 2002/05/30 17:44:35 1.46 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.40 2002/04/19 13:46:19 lievre Exp $ +/* $Id: imach.c,v 1.46 2002/05/30 17:44:35 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*/ @@ -96,7 +96,7 @@ double jmean; /* Mean space between 2 wa double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; -FILE *ficgp,*ficresprob,*ficpop; +FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor; FILE *ficreseij; char filerese[FILENAMELENGTH]; FILE *ficresvij; @@ -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,110 @@ 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 ******************/ @@ -1621,7 +1720,7 @@ void varevsij(char fileres[], double *** double age,agelim, hf; int theta; - fprintf(ficresvij,"# Covariances of life expectancies\n"); + fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n"); fprintf(ficresvij,"# Age"); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) @@ -1756,7 +1855,7 @@ void varprevlim(char fileres[], double * double age,agelim; int theta; - fprintf(ficresvpl,"# Standard deviation of prevalences limit\n"); + fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n"); fprintf(ficresvpl,"# Age"); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %1d-%1d",i,i); @@ -1828,7 +1927,7 @@ void varprevlim(char fileres[], double * void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax) { int i, j, i1, k1, j1, z1; - int k=0, cptcode; + int k=0,l, cptcode; double **dnewm,**doldm; double *xp; double *gp, *gm; @@ -1836,15 +1935,43 @@ void varprob(char fileres[], double **ma double age,agelim, cov[NCOVMAX]; int theta; char fileresprob[FILENAMELENGTH]; + char fileresprobcov[FILENAMELENGTH]; + char fileresprobcor[FILENAMELENGTH]; strcpy(fileresprob,"prob"); strcat(fileresprob,fileres); 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); - - + strcpy(fileresprobcov,"probcov"); + strcat(fileresprobcov,fileres); + if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { + printf("Problem with resultfile: %s\n", fileresprobcov); + } + strcpy(fileresprobcor,"probcor"); + strcat(fileresprobcor,fileres); + if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { + printf("Problem with resultfile: %s\n", fileresprobcor); + } + printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); + printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); + printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); + + fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); + fprintf(ficresprob,"# Age"); + fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); + fprintf(ficresprobcov,"# Age"); + fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); + fprintf(ficresprobcov,"# Age"); + for(i=1; i<=nlstate;i++) + for(j=1; j<=(nlstate+ndeath);j++){ + fprintf(ficresprob," p%1d-%1d (SE)",i,j); + fprintf(ficresprobcov," p%1d-%1d ",i,j); + fprintf(ficresprobcor," p%1d-%1d ",i,j); + } + fprintf(ficresprob,"\n"); + fprintf(ficresprobcov,"\n"); + fprintf(ficresprobcor,"\n"); xp=vector(1,npar); dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath)); @@ -1859,15 +1986,20 @@ void varprob(char fileres[], double **ma if (cptcovn>0) { fprintf(ficresprob, "\n#********** Variable "); + fprintf(ficresprobcov, "\n#********** Variable "); + fprintf(ficresprobcor, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); fprintf(ficresprob, "**********\n#"); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprobcov, "**********\n#"); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprobcor, "**********\n#"); } for (age=bage; age<=fage; age ++){ cov[2]=age; for (k=1; k<=cptcovn;k++) { cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]]; - } for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; for (k=1; k<=cptcovprod;k++) @@ -1925,16 +2057,33 @@ void varprob(char fileres[], double **ma } } - /*printf("\n%d ",(int)age); + /*printf("\n%d ",(int)age); for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){ printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); }*/ fprintf(ficresprob,"\n%d ",(int)age); + fprintf(ficresprobcov,"\n%d ",(int)age); + fprintf(ficresprobcor,"\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,"%12.3e (%12.3e) ",gm[i],sqrt(doldm[i][j])); + for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ + fprintf(ficresprobcov,"%12.3e ",gm[i]); + fprintf(ficresprobcor,"%12.3e ",gm[i]); + } + i=0; + for (k=1; k<=(nlstate);k++){ + for (l=1; l<=(nlstate+ndeath);l++){ + i=i++; + fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l); + fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l); + for (j=1; j<=i;j++){ + fprintf(ficresprobcov," %12.3e",doldm[i][j]); + fprintf(ficresprobcor," %12.3e",doldm[i][j]/sqrt(doldm[i][i])/sqrt(doldm[j][j])); + } + } + } } } free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); @@ -1944,15 +2093,19 @@ void varprob(char fileres[], double **ma } free_vector(xp,1,npar); fclose(ficresprob); - + fclose(ficresprobcov); + fclose(ficresprobcor); } + /******************* Printing html file ***********/ void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \ - int lastpass, int stepm, int weightopt, char model[],\ - int imx,int jmin, int jmax, double jmeanint,char optionfile[], \ - char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\ - char version[], int popforecast, int estepm ){ + int lastpass, int stepm, int weightopt, char model[],\ + int imx,int jmin, int jmax, double jmeanint,char optionfile[], \ + char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\ + char version[], int popforecast, int estepm ,\ + double jprev1, double mprev1,double anprev1, \ + double jprev2, double mprev2,double anprev2){ int jj1, k1, i1, cpt; FILE *fichtm; /*char optionfilehtm[FILENAMELENGTH];*/ @@ -1963,25 +2116,32 @@ void printinghtml(char fileres[], char t printf("Problem with %s \n",optionfilehtm), exit(0); } - fprintf(fichtm," %s
\n + fprintf(fichtm," %s
\n Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n \n Total number of observations=%d
\n Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n -
-