/* $Id$
$State$
$Log$
+ Revision 1.207 2015/10/27 17:36:57 brouard
+ *** empty log message ***
+
Revision 1.206 2015/10/24 07:14:11 brouard
*** empty log message ***
/* double **matprod2(); */ /* test */
double **out, cov[NCOVMAX+1], **pmij();
double **newm;
- double agefin, delaymax=100 ; /* Max number of years to converge */
+ double agefin, delaymax=100. ; /* 100 Max number of years to converge */
int ncvloop=0;
for (ii=1;ii<=nlstate+ndeath;ii++)
prlim[i][j]= newm[i][j]/(1-sumnew);
max=FMAX(max,prlim[i][j]);
min=FMIN(min,prlim[i][j]);
- printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min);
+ /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */
}
maxmin=(max-min)/(max+min)*2;
maxmax=FMAX(maxmax,maxmin);
+ /* for(i=1; i<=nlstate; i++) { */
+ /* sumnew=0.; */
+ /* sumnew+=prlim[i][j]; */
+ /* } */
+ /* prmimj = sumnew/(float)nlstate; /\* Means of various prevalence limits. */
} /* j loop */
*ncvyear= (int)age- (int)agefin;
- printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear);
+ /* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
if(maxmax < ftolpl){
/* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
return prlim;
}
} /* age loop */
- printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g. \n\
-Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
-/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
+ /* After some age loop it doesn't converge */
+ printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g within %.0f years. \n\
+Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
+ /* printf(" age= %d newm\n",(int)age); */
+ /* for(i=1; i<=nlstate+ndeath; i++) { */
+ /* for(j=1;j<=nlstate+ndeath;j++){ */
+ /* printf(" %lf", newm[i][j]); */
+ /* } */
+ /* printf("\n"); */
+ /* } */
+ /* printf("\n"); */
+ /* printf("prlim\n"); */
+ /* for(i=1; i<=nlstate; i++) { */
+ /* sumnew=0; */
+ /* for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; */
+ /* for(j=1;j<=nlstate;j++){ */
+ /* prlim[i][j]= newm[i][j]/(1-sumnew); */
+ /* printf(" %lf", prlim[i][j]); */
+ /* } */
+ /* printf("\n"); */
+ /* } */
+ /* printf("\n"); */
+
+ /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
return prlim; /* should not reach here */
}
fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);
fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
- fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
-<img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
- fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
-<img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
- fflush(fichtm);
for (k=1; k<= nlstate ; k++) {
fprintf(fichtm,"<br>- Probability p%dj by origin %d and destination j <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
<img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
}
+ fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
+<img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
+<img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ fflush(fichtm);
}
return;
}
double p2[MAXPARM+1];
int k, kmax=1;
double v1, v2, cv12, lc1, lc2;
+
+ int firstime=0;
fx=func(x);
for (k=1; k<=kmax; k=k+10) {
k4=func(p2)-fx;
res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */
if(k1*k2*k3*k4 <0.){
+ firstime=1;
kmax=kmax+10;
- if(kmax >=10){
+ }
+ if(kmax >=10 || firstime ==1){
printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
- }
}
#ifdef DEBUGHESSIJ
v1=hess[thetai][thetai];
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
}
printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
-
fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
pstamp(ficresprobmorprev);
fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);
}
fprintf(ficresprobmorprev,"\n");
+
fprintf(ficgp,"\n# Routine varevsij");
fprintf(ficgp,"\nunset title \n");
/* fprintf(fichtm, "#Local time at start: %s", strstart);*/
/* For example 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 function hpijx to understand why (it is linked to memory size questions) */
- /* We decided (b) to get a life expectancy respecting the most precise curvature of the
+ nstepm is the number of stepm from age to agelim.
+ Look at function hpijx to understand why (it is linked to memory size questions)
+ we decided (b) to get a life expectancy respecting the most precise curvature of the
survival function given by stepm (the optimization length). Unfortunately it
means that if the survival funtion is printed every two years of age and if
you sum them up and add 1 year (area under the trapezoids) you won't get the same
double *xp;
double *gp, *gm;
double **gradg, **trgradg;
+ double **mgm, **mgp;
double age,agelim;
int theta;
nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
if (stepm >= YEARM) hstepm=1;
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
gradg=matrix(1,npar,1,nlstate);
+ mgp=matrix(1,npar,1,nlstate);
+ mgm=matrix(1,npar,1,nlstate);
gp=vector(1,nlstate);
gm=vector(1,nlstate);
for(i=1; i<=npar; i++){ /* Computes gradient */
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
+ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Missing or not useful because 1 year */
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
- for(i=1;i<=nlstate;i++)
+ for(i=1;i<=nlstate;i++){
gp[i] = prlim[i][i];
-
+ mgp[theta][i] = prlim[i][i];
+ }
for(i=1; i<=npar; i++) /* Computes gradient */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
- for(i=1;i<=nlstate;i++)
+ for(i=1;i<=nlstate;i++){
gm[i] = prlim[i][i];
-
+ mgm[theta][i] = prlim[i][i];
+ }
for(i=1;i<=nlstate;i++)
gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta];
} /* End theta */
for(j=1; j<=nlstate;j++)
for(theta=1; theta <=npar; theta++)
trgradg[j][theta]=gradg[theta][j];
+ if((int)age==68 ||(int)age== 69 ){
+ printf("\nmgm mgp %d ",(int)age);
+ for(j=1; j<=nlstate;j++){
+ printf("%d ",j);
+ for(theta=1; theta <=npar; theta++)
+ printf("%d %lf %lf",theta,mgm[theta][j],mgp[theta][j]);
+ printf("\n ");
+ }
+ }
+ if((int)age==68 ||(int)age== 69 ){
+ printf("\n gradg %d ",(int)age);
+ for(j=1; j<=nlstate;j++){
+ printf("%d ",j);
+ for(theta=1; theta <=npar; theta++)
+ printf("%d %lf ",theta,gradg[theta][j]);
+ printf("\n ");
+ }
+ }
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] =0.;
- if((int)age==67 ||(int)age== 66 ){
+ if((int)age==68 ||(int)age== 69 ){
matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
}else{
fprintf(ficresvpl,"\n");
free_vector(gp,1,nlstate);
free_vector(gm,1,nlstate);
+ free_matrix(mgm,1,npar,1,nlstate);
+ free_matrix(mgp,1,npar,1,nlstate);
free_matrix(gradg,1,npar,1,nlstate);
free_matrix(trgradg,1,nlstate,1,npar);
} /* End age */
}
/* State specific survival functions (period) */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"<br>\n- Survival functions from state %d in any different live states and total.\
+ fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\
Or probability to survive in various states (1 to %d) being in state %d at different ages.\
<a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);
}
printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
}
- printf("Computing Health Expectancies: result on file '%s' \n", filerese);
- fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
+ printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
+ fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
/*}*/
}
fclose(ficreseij);
-
+ printf("done evsij\n");fflush(stdout);
+ fprintf(ficlog,"done evsij\n");fflush(ficlog);
/*---------- Health expectancies and variances ------------*/
printf("Problem with total LE resultfile: %s\n", filerest);goto end;
fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
}
- printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
- fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
+ printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);
+ fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);
strcpy(fileresstde,"STDE_");
printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
}
- printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
- fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+ printf(" Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+ fprintf(ficlog," Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
strcpy(filerescve,"CVE_");
strcat(filerescve,fileresu);
printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
}
- printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
- fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+ printf(" Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+ fprintf(ficlog," Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
strcpy(fileresv,"V_");
strcat(fileresv,fileresu);
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
}
- printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
- fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
+ printf(" Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(stdout);
+ fprintf(ficlog," Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(ficlog);
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
for (k=1; k <= (int) pow(2,cptcoveff); k++){
- fprintf(ficrest,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficrest,"******\n");
-
- fprintf(ficresstdeij,"\n#****** ");
- fprintf(ficrescveij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- fprintf(ficresstdeij,"******\n");
- fprintf(ficrescveij,"******\n");
-
- fprintf(ficresvij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresvij,"******\n");
-
- eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
- cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
- /*
- */
- /* goto endfree; */
-
- vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- pstamp(ficrest);
-
-
- for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
- oldm=oldms;savm=savms; /* ZZ Segmentation fault */
- cptcod= 0; /* To be deleted */
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
- fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
- if(vpopbased==1)
- fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
- else
- fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
- fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
- for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
- fprintf(ficrest,"\n");
- /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
- epj=vector(1,nlstate+1);
- for(age=bage; age <=fage ;age++){
- prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */
- if (vpopbased==1) {
- if(mobilav ==0){
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=probs[(int)age][i][k];
- }else{ /* mobilav */
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=mobaverage[(int)age][i][k];
- }
- }
-
- fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
- /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
- /* printf(" age %4.0f ",age); */
- for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
- for(i=1, epj[j]=0.;i <=nlstate;i++) {
- epj[j] += prlim[i][i]*eij[i][j][(int)age];
- /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
- /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
- }
- epj[nlstate+1] +=epj[j];
+ fprintf(ficrest,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficrest,"******\n");
+
+ fprintf(ficresstdeij,"\n#****** ");
+ fprintf(ficrescveij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ fprintf(ficresstdeij,"******\n");
+ fprintf(ficrescveij,"******\n");
+
+ fprintf(ficresvij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresvij,"******\n");
+
+ eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+ printf(" cvevsij %d, ",k);
+ fprintf(ficlog, " cvevsij %d, ",k);
+ cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
+ printf(" end cvevsij \n ");
+ fprintf(ficlog, " end cvevsij \n ");
+
+ /*
+ */
+ /* goto endfree; */
+
+ vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+ pstamp(ficrest);
+
+
+ for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+ oldm=oldms;savm=savms; /* ZZ Segmentation fault */
+ cptcod= 0; /* To be deleted */
+ printf("varevsij %d \n",vpopbased);
+ fprintf(ficlog, "varevsij %d \n",vpopbased);
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
+ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
+ if(vpopbased==1)
+ fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
+ else
+ fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
+ fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
+ for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+ fprintf(ficrest,"\n");
+ /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
+ epj=vector(1,nlstate+1);
+ printf("Computing age specific period (stable) prevalences in each health state \n");
+ fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
+ for(age=bage; age <=fage ;age++){
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */
+ if (vpopbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][k];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][k];
}
- /* printf(" age %4.0f \n",age); */
-
- for(i=1, vepp=0.;i <=nlstate;i++)
- for(j=1;j <=nlstate;j++)
- vepp += vareij[i][j][(int)age];
- fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
- for(j=1;j <=nlstate;j++){
- fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
+ }
+
+ fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+ /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
+ /* printf(" age %4.0f ",age); */
+ for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
+ for(i=1, epj[j]=0.;i <=nlstate;i++) {
+ epj[j] += prlim[i][i]*eij[i][j][(int)age];
+ /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+ /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
}
- fprintf(ficrest,"\n");
+ epj[nlstate+1] +=epj[j];
+ }
+ /* printf(" age %4.0f \n",age); */
+
+ for(i=1, vepp=0.;i <=nlstate;i++)
+ for(j=1;j <=nlstate;j++)
+ vepp += vareij[i][j][(int)age];
+ fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+ for(j=1;j <=nlstate;j++){
+ fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
}
+ fprintf(ficrest,"\n");
}
- free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
- free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
- free_vector(epj,1,nlstate+1);
+ } /* End vpopbased */
+ free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+ free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+ free_vector(epj,1,nlstate+1);
+ printf("done \n");fflush(stdout);
+ fprintf(ficlog,"done\n");fflush(ficlog);
+
/*}*/
- }
+ } /* End k */
free_vector(weight,1,n);
free_imatrix(Tvard,1,NCOVMAX,1,2);
free_imatrix(s,1,maxwav+1,1,n);
fclose(ficrescveij);
fclose(ficresvij);
fclose(ficrest);
+ printf("done Health expectancies\n");fflush(stdout);
+ fprintf(ficlog,"done Health expectancies\n");fflush(ficlog);
fclose(ficpar);
/*------- Variance of period (stable) prevalence------*/
printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl);
exit(0);
}
- printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl);
+ printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
+ fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
}
fclose(ficresvpl);
+ printf("done variance-covariance of period prevalence\n");fflush(stdout);
+ fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
/*---------- End : free ----------------*/
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
git://henry.ined.fr / .git / commitdiff