version 1.100, 2004/07/12 18:29:06
|
version 1.103, 2005/09/30 15:54:49
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.103 2005/09/30 15:54:49 lievre |
|
(Module): sump fixed, loop imx fixed, and simplifications. |
|
|
|
Revision 1.102 2004/09/15 17:31:30 brouard |
|
Add the possibility to read data file including tab characters. |
|
|
|
Revision 1.101 2004/09/15 10:38:38 brouard |
|
Fix on curr_time |
|
|
Revision 1.100 2004/07/12 18:29:06 brouard |
Revision 1.100 2004/07/12 18:29:06 brouard |
Add version for Mac OS X. Just define UNIX in Makefile |
Add version for Mac OS X. Just define UNIX in Makefile |
|
|
Line 241
|
Line 250
|
/* $Id$ */ |
/* $Id$ */ |
/* $State$ */ |
/* $State$ */ |
|
|
char version[]="Imach version 0.97b, May 2004, INED-EUROREVES "; |
char version[]="Imach version 0.97c, September 2004, INED-EUROREVES "; |
char fullversion[]="$Revision$ $Date$"; |
char fullversion[]="$Revision$ $Date$"; |
int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ |
int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ |
int nvar; |
int nvar; |
Line 436 int nbocc(char *s, char occ)
|
Line 445 int nbocc(char *s, char occ)
|
|
|
void cutv(char *u,char *v, char*t, char occ) |
void cutv(char *u,char *v, char*t, char occ) |
{ |
{ |
/* cuts string t into u and v where u is ended by char occ excluding it |
/* cuts string t into u and v where u ends before first occurence of char 'occ' |
and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2) |
and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2') |
gives u="abcedf" and v="ghi2j" */ |
gives u="abcedf" and v="ghi2j" */ |
int i,lg,j,p=0; |
int i,lg,j,p=0; |
i=0; |
i=0; |
Line 879 void powell(double p[], double **xi, int
|
Line 888 void powell(double p[], double **xi, int
|
fprintf(ficrespow,"\n");fflush(ficrespow); |
fprintf(ficrespow,"\n");fflush(ficrespow); |
if(*iter <=3){ |
if(*iter <=3){ |
tm = *localtime(&curr_time.tv_sec); |
tm = *localtime(&curr_time.tv_sec); |
strcpy(strcurr,asctime(&tmf)); |
strcpy(strcurr,asctime(&tm)); |
/* asctime_r(&tm,strcurr); */ |
/* asctime_r(&tm,strcurr); */ |
forecast_time=curr_time; |
forecast_time=curr_time; |
itmp = strlen(strcurr); |
itmp = strlen(strcurr); |
if(strcurr[itmp-1]=='\n') |
if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ |
strcurr[itmp-1]='\0'; |
strcurr[itmp-1]='\0'; |
printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); |
printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); |
fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); |
fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); |
Line 895 void powell(double p[], double **xi, int
|
Line 904 void powell(double p[], double **xi, int
|
itmp = strlen(strfor); |
itmp = strlen(strfor); |
if(strfor[itmp-1]=='\n') |
if(strfor[itmp-1]=='\n') |
strfor[itmp-1]='\0'; |
strfor[itmp-1]='\0'; |
printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); |
printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); |
fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); |
fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); |
} |
} |
} |
} |
for (i=1;i<=n;i++) { |
for (i=1;i<=n;i++) { |
Line 1235 double func( double *x)
|
Line 1244 double func( double *x)
|
} /* end mult */ |
} /* end mult */ |
|
|
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ |
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ |
/* But now since version 0.9 we anticipate for bias and large stepm. |
/* But now since version 0.9 we anticipate for bias at large stepm. |
* If stepm is larger than one month (smallest stepm) and if the exact delay |
* If stepm is larger than one month (smallest stepm) and if the exact delay |
* (in months) between two waves is not a multiple of stepm, we rounded to |
* (in months) between two waves is not a multiple of stepm, we rounded to |
* the nearest (and in case of equal distance, to the lowest) interval but now |
* the nearest (and in case of equal distance, to the lowest) interval but now |
* we keep into memory the bias bh[mi][i] and also the previous matrix product |
* we keep into memory the bias bh[mi][i] and also the previous matrix product |
* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the |
* (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the |
* probability in order to take into account the bias as a fraction of the way |
* probability in order to take into account the bias as a fraction of the way |
* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies |
* from savm to out if bh is negative or even beyond if bh is positive. bh varies |
* -stepm/2 to stepm/2 . |
* -stepm/2 to stepm/2 . |
* For stepm=1 the results are the same as for previous versions of Imach. |
* For stepm=1 the results are the same as for previous versions of Imach. |
* For stepm > 1 the results are less biased than in previous versions. |
* For stepm > 1 the results are less biased than in previous versions. |
Line 1250 double func( double *x)
|
Line 1259 double func( double *x)
|
s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
bbh=(double)bh[mi][i]/(double)stepm; |
bbh=(double)bh[mi][i]/(double)stepm; |
/* bias is positive if real duration |
/* bias bh is positive if real duration |
* is higher than the multiple of stepm and negative otherwise. |
* is higher than the multiple of stepm and negative otherwise. |
*/ |
*/ |
/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ |
/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ |
if( s2 > nlstate){ |
if( s2 > nlstate){ |
/* i.e. if s2 is a death state and if the date of death is known then the contribution |
/* i.e. if s2 is a death state and if the date of death is known then the contribution |
to the likelihood is the probability to die between last step unit time and current |
to the likelihood is the probability to die between last step unit time and current |
step unit time, which is also the differences between probability to die before dh |
step unit time, which is also equal to probability to die before dh |
and probability to die before dh-stepm . |
minus probability to die before dh-stepm . |
In version up to 0.92 likelihood was computed |
In version up to 0.92 likelihood was computed |
as if date of death was unknown. Death was treated as any other |
as if date of death was unknown. Death was treated as any other |
health state: the date of the interview describes the actual state |
health state: the date of the interview describes the actual state |
Line 3915 double gompertz(double x[])
|
Line 3924 double gompertz(double x[])
|
int i,n=0; /* n is the size of the sample */ |
int i,n=0; /* n is the size of the sample */ |
for (i=0;i<=imx-1 ; i++) { |
for (i=0;i<=imx-1 ; i++) { |
sump=sump+weight[i]; |
sump=sump+weight[i]; |
sump=sump+1; |
/* sump=sump+1;*/ |
num=num+1; |
num=num+1; |
} |
} |
|
|
|
|
/* for (i=1; i<=imx; i++) |
/* for (i=0; i<=imx; i++) |
if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ |
if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ |
|
|
for (i=0;i<=imx-1 ; i++) |
for (i=1;i<=imx ; i++) |
{ |
{ |
if (cens[i]==1 & wav[i]>1) |
if (cens[i]==1 & wav[i]>1) |
A=-x[1]/(x[2])* |
A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); |
(exp(x[2]/YEARM*(agecens[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))); |
|
|
|
if (cens[i]==0 & wav[i]>1) |
if (cens[i]==0 & wav[i]>1) |
A=-x[1]/(x[2])* |
A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) |
(exp(x[2]/YEARM*(agedc[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))) |
+log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); |
+log(x[1]/YEARM)+x[2]/YEARM*(agedc[i]*12-agegomp*12)+log(YEARM); |
|
|
|
if (wav[i]>1 & agecens[i]>15) { |
if (wav[i]>1 & agecens[i]>15) { |
L=L+A*weight[i]; |
L=L+A*weight[i]; |
Line 4396 int main(int argc, char *argv[])
|
Line 4403 int main(int argc, char *argv[])
|
i=1; |
i=1; |
while (fgets(line, MAXLINE, fic) != NULL) { |
while (fgets(line, MAXLINE, fic) != NULL) { |
if ((i >= firstobs) && (i <=lastobs)) { |
if ((i >= firstobs) && (i <=lastobs)) { |
|
for(j=0; line[j] != '\n';j++){ /* Untabifies line */ |
|
if(line[j] == '\t') |
|
line[j] = ' '; |
|
} |
for (j=maxwav;j>=1;j--){ |
for (j=maxwav;j>=1;j--){ |
cutv(stra, strb,line,' '); s[j][i]=atoi(strb); |
cutv(stra, strb,line,' '); s[j][i]=atoi(strb); |
strcpy(line,stra); |
strcpy(line,stra); |