--- imach/src/imach.c 2002/03/08 16:17:18 1.30 +++ imach/src/imach.c 2002/03/10 13:43:02 1.31 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.30 2002/03/08 16:17:18 lievre Exp $ +/* $Id: imach.c,v 1.31 2002/03/10 13:43:02 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -1527,7 +1527,7 @@ void tricode(int *Tvar, int **nbcode, in void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij) { /* Health expectancies */ - int i, j, nhstepm, hstepm, h; + int i, j, nhstepm, hstepm, h, nstepm, k; double age, agelim,hf; double ***p3mat; @@ -1538,33 +1538,40 @@ void evsij(char fileres[], double ***eij fprintf(ficreseij," %1d-%1d",i,j); fprintf(ficreseij,"\n"); - hstepm=1*YEARM; /* Every j years of age (in month) */ - hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */ + k=1; /* For example stepm=6 months */ + hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */ + hstepm=1; /* or (b) We decided to compute the life expectancy with the smallest unit */ + /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. + nhstepm is the number of hstepm from age to agelim + nstepm is the number of stepm from age to agelin. + Look at hpijx to understand the reason of that which relies in memory size + and note for a fixed period like k years */ + /* We decided (b) to get a life expectancy respecting the most precise curvature of the + survival function given par stepm (the optimization length). Unfortunately it + means that if the survival funtion is printed only each two years of age and if + you sum them up and add 1 year (area under the trapezoids) you won't get the same + results. So we changed our mind and took the option of the best precision. + */ + hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */ agelim=AGESUP; for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ /* nhstepm age range expressed in number of stepm */ - nhstepm=(int) rint((agelim-age)*YEARM/stepm); - /* Typically if 20 years = 20*12/6=40 stepm */ - /*if (stepm >= YEARM) hstepm=1;*/ -hstepm=1; - nhstepm = nhstepm/hstepm;/* Expressed in hstepm, typically 40/4=10 */ + nstepm=(int) rint((agelim-age)*YEARM/stepm); + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + 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); /* 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); -hf=stepm/YEARM; -/*printf("stepm=%d nhstepm=%d hstepm=%d age=%lf ",stepm, nhstepm, hstepm, age);*/ + hf=hstepm/YEARM; /* Duration of hstepm expressed in year unit. */ 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] +=hf*(p3mat[i][j][h]+p3mat[i][j][h+1])/2.0; + 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]);*/ } - - /* hf=1; - if (stepm >= YEARM) hf=stepm/YEARM;*/ -hf=stepm/YEARM; - fprintf(ficreseij,"%3.0f",age ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ @@ -1896,7 +1903,7 @@ void printinghtml(char fileres[], char t printf("Problem with %s \n",optionfilehtm), exit(0); } - fprintf(fichtm,"