/* $Id$
$State$
$Log$
+ Revision 1.86 2003/06/17 20:04:08 brouard
+ (Module): Change position of html and gnuplot routines and added
+ routine fileappend.
+
Revision 1.85 2003/06/17 13:12:43 brouard
* imach.c (Repository): Check when date of death was earlier that
current date of interview. It may happen when the death was just
/* $Id$ */
/* $State$ */
-char version[]="Imach version 0.95a2, June 2003, INED-EUROREVES ";
+char version[]="Imach version 0.96, June 2003, INED-EUROREVES ";
char fullversion[]="$Revision$ $Date$";
int erreur; /* Error number */
int nvar;
int *wav; /* Number of waves for this individuual 0 is possible */
int maxwav; /* Maxim number of waves */
int jmin, jmax; /* min, max spacing between 2 waves */
+int gipmx, gsw; /* Global variables on the number of contributions
+ to the likelihood and the sum of weights (done by funcone)*/
int mle, weightopt;
int **mw; /* mw[mi][i] is number of the mi wave for this individual */
int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
/*************** log-likelihood *************/
double funcone( double *x)
{
+ /* Same as likeli but slower because of a lot of printf and if */
int i, ii, j, k, mi, d, kk;
double l, ll[NLSTATEMAX], cov[NCOVMAX];
double **out;
double lli; /* Individual log likelihood */
+ double llt;
int s1, s2;
double bbh, survp;
/*extern weight */
%10.6f %10.6f %10.6f ", \
num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
- for(k=1,l=0.; k<=nlstate; k++)
- fprintf(ficresilk," %10.6f",ll[k]);
- fprintf(ficresilk,"\n");
+ for(k=1,llt=0.,l=0.; k<=nlstate; k++){
+ llt +=ll[k]*gipmx/gsw;
+ fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
+ }
+ fprintf(ficresilk," %10.6f\n", -llt);
}
} /* end of wave */
} /* end of individual */
for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */
l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
+ if(globpr==0){ /* First time we count the contributions and weights */
+ gipmx=ipmx;
+ gsw=sw;
+ }
return -l;
}
-void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpr, long *ipmx, double *sw, double *fretone, double (*funcone)(double []))
+void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double []))
{
- /* This routine should help understanding what is done with the selection of individuals/waves and
+ /* This routine should help understanding what is done with
+ the selection of individuals/waves and
to check the exact contribution to the likelihood.
Plotting could be done.
*/
int k;
- if(globpr !=0){ /* Just counts and sums no printings */
+
+ if(*globpri !=0){ /* Just counts and sums no printings */
strcpy(fileresilk,"ilk");
strcat(fileresilk,fileres);
if((ficresilk=fopen(fileresilk,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresilk);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
}
- fprintf(ficresilk, "#individual(line's record) s1 s2 wave# effective_wave# number_of_product_matrix pij weight 2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state\n");
+ fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight out sav ");
/* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
for(k=1; k<=nlstate; k++)
- fprintf(ficresilk," ll[%d]",k);
- fprintf(ficresilk,"\n");
+ fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
+ fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n");
}
*fretone=(*funcone)(p);
- if(globpr !=0){
+ if(*globpri !=0){
fclose(ficresilk);
- if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
- printf("Problem with html file: %s\n", optionfilehtm);
- fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);
- exit(0);
- }
- else{
- fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",fileresilk);
- fclose(fichtm);
- }
- }
+ fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",fileresilk,fileresilk);
+ fflush(fichtm);
+ }
return;
}
else{
fprintf(ficgp,"\n# Routine varevsij");
}
- if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
- printf("Problem with html file: %s\n", optionfilehtm);
- fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);
- exit(0);
- }
- else{
- fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
- fprintf(fichtm,"\n<br>%s <br>\n",digitp);
- }
+/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */
+/* printf("Problem with html file: %s\n", optionfilehtm); */
+/* fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); */
+/* exit(0); */
+/* } */
+/* else{ */
+ fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
+ fprintf(fichtm,"\n<br>%s <br>\n",digitp);
+/* } */
varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
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");
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
fclose(ficresprobmorprev);
fclose(ficgp);
- fclose(fichtm);
+/* fclose(fichtm); */
} /* end varevsij */
/************ Variance of prevlim ******************/
else{
fprintf(ficgp,"\n# Routine varprob");
}
- if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
- printf("Problem with html file: %s\n", optionfilehtm);
- fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);
- exit(0);
- }
- else{
+/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */
+/* printf("Problem with html file: %s\n", optionfilehtm); */
+/* fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); */
+/* exit(0); */
+/* } */
+/* else{ */
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
fprintf(fichtm,"\n");
fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the p<inf>ij</inf>, p<inf>kl</inf> to understand the covariance between two incidences. They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
fprintf(fichtm,"\n<br> We have drawn x'cov<sup>-1</sup>x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis. <br> When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.<br> \n");
- }
+/* } */
cov[1]=1;
tj=cptcoveff;
fclose(ficresprobcov);
fclose(ficresprobcor);
fclose(ficgp);
- fclose(fichtm);
}
double jprev2, double mprev2,double anprev2){
int jj1, k1, i1, cpt;
/*char optionfilehtm[FILENAMELENGTH];*/
- if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
- printf("Problem with %s \n",optionfilehtm), exit(0);
- fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0);
- }
+/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */
+/* printf("Problem with %s \n",optionfilehtm), exit(0); */
+/* fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0); */
+/* } */
fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \
- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n \
} /* end i1 */
}/* End k1 */
fprintf(fichtm,"</ul>");
-fclose(fichtm);
+ fflush(fichtm);
}
/******************* Gnuplot file **************/
fclose(ficrespop);
} /* End of popforecast */
-int fileappend(FILE *fichier, char *optionfile)
+int fileappend(FILE *fichier, char *optionfich)
{
- if((fichier=fopen(optionfile,"a"))==NULL) {
- printf("Problem with file: %s\n", optionfile);
- fprintf(ficlog,"Problem with file: %s\n", optionfile);
- return (1);
+ if((fichier=fopen(optionfich,"a"))==NULL) {
+ printf("Problem with file: %s\n", optionfich);
+ fprintf(ficlog,"Problem with file: %s\n", optionfich);
+ return (0);
}
-
+ fflush(fichier);
+ return (1);
}
/***********************************************/
/**************** Main Program *****************/
int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
int jj;
int numlinepar=0; /* Current linenumber of parameter file */
+ /* FILE *fichtm; *//* Html File */
+ /* FILE *ficgp;*/ /*Gnuplot File */
double agedeb, agefin,hf;
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
fclose(ficgp);
/*--------- index.htm --------*/
- strcpy(optionfilehtm,optionfile);
+ strcpy(optionfilehtm,optionfilefiname);
strcat(optionfilehtm,".htm");
if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
printf("Problem with %s \n",optionfilehtm), exit(0);
}
- fprintf(fichtm,"<body> <font size=\"2\">%s <br> %s</font> \
+ fprintf(fichtm,"<body>\n<title>IMaCh %s</title>\n <font size=\"2\">%s <br> %s</font> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\
\n\
<ul><li><h4>Parameter files</h4>\n\
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\
- Log file of the run: <a href=\"%s\">%s</a><br>\n\
- - Gnuplot file name: <a href=\"%s\">%s</a>\n\
+ - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\
- Date and time at start: %s</ul>\n",\
- version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt,\
+ fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt,\
model,fileres,fileres,\
filelog,filelog,optionfilegnuplot,optionfilegnuplot,strt);
- fclose(fichtm);
+ /*fclose(fichtm);*/
+ fflush(fichtm);
/* Calculates basic frequencies. Computes observed prevalence at single age
and prints on file fileres'p'. */
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);
- if(fileappend(fichtm, optionfilehtm)){
+/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */
+/* printf("Problem with file: %s\n", optionfilehtm); */
+/* fprintf(ficlog,"Problem with file: %s\n", optionfilehtm); */
+/* } */
+
+
+/* if(fileappend(fichtm, optionfilehtm)){ */
+ fprintf(fichtm,"\n");
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
imx,agemin,agemax,jmin,jmax,jmean);
- fclose(fichtm);
- }
+/* fclose(fichtm); */
+/* } */
+
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);
fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);
/* printf("Total time was %d uSec.\n", total_usecs);*/
- if(fileappend(fichtm,optionfilehtm)){
- fprintf(fichtm,"<br>Localtime at start %s and at end=%s<br>",strt, strtend);
- fclose(fichtm);
- }
+/* if(fileappend(fichtm,optionfilehtm)){ */
+ fprintf(fichtm,"<br>Localtime at start %s and at end=%s<br>",strt, strtend);
+ fclose(fichtm);
/*------ End -----------*/
end:
git://henry.ined.fr / .git / commitdiff