/* $Id$
$State$
$Log$
+ Revision 1.301 2019/06/04 13:51:20 brouard
+ Summary: Error in 'r'parameter file backcast yearsbproj instead of yearsfproj
+
Revision 1.300 2019/05/22 19:09:45 brouard
Summary: version 0.99r19 of May 2019
int cptcov=0; /* Working variable */
int nobs=10; /* Number of observations in the data lastobs-firstobs */
int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */
-int npar=NPARMAX;
+int npar=NPARMAX; /* Number of parameters (nlstate+ndeath-1)*nlstate*ncovmodel; */
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */
/* New value of last point Pn is not computed, P(n-1) */
for(j=1;j<=n;j++) {
- if(flatdir[j] >0){
- printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- }
- /* printf("\n"); */
- /* fprintf(ficlog,"\n"); */
- }
+ if(flatdir[j] >0){
+ printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ }
+ /* printf("\n"); */
+ /* fprintf(ficlog,"\n"); */
+ }
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
/* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */
double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij )
{
- /* Computes the backward probability at age agefin and covariate combination ij. In fact cov is already filled and x too.
+ /* Computes the backward probability at age agefin, cov[2], and covariate combination 'ij'. In fact cov is already filled and x too.
* Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij.
*/
int i, ii, j,k;
#else
if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
if(firsthree == 0){
- printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath);
+ printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath);
firsthree=1;
}
- fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath);
+ fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath);
mw[++mi][i]=m;
mli=m;
}
*/
int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
double agec; /* generic age */
- double agelim, ppij, ppi, yp,yp1,yp2,jintmean,mintmean,aintmean;
+ double agelim, ppij, ppi, yp,yp1,yp2; /* ,jintmean,mintmean,aintmean;*/
double *popeffectif,*popcount;
double ***p3mat;
/* double ***mobaverage; */
/******************* Gompertz Likelihood ******************************/
double gompertz(double x[])
{
- double A,B,L=0.0,sump=0.,num=0.;
+ double A=0.0,B=0.,L=0.0,sump=0.,num=0.;
int i,n=0; /* n is the size of the sample */
for (i=1;i<=imx ; i++) {
/* sump=sump+1;*/
num=num+1;
}
-
-
+ L=0.0;
+ /* agegomp=AGEGOMP; */
/* 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]);*/
- for (i=1;i<=imx ; i++)
- {
- if (cens[i] == 1 && wav[i]>1)
- A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
-
- if (cens[i] == 0 && wav[i]>1)
+ for (i=1;i<=imx ; i++) {
+ /* mu(a)=mu(agecomp)*exp(teta*(age-agegomp))
+ mu(a)=x[1]*exp(x[2]*(age-agegomp)); x[1] and x[2] are per year.
+ * L= Product mu(agedeces)exp(-\int_ageexam^agedc mu(u) du ) for a death between agedc (in month)
+ * and agedc +1 month, cens[i]=0: log(x[1]/YEARM)
+ * +
+ * exp(-\int_ageexam^agecens mu(u) du ) when censored, cens[i]=1
+ */
+ if (wav[i] > 1 || agedc[i] < AGESUP) {
+ if (cens[i] == 1){
+ A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
+ } else if (cens[i] == 0){
A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
- +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM);
-
+ +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+ } else
+ printf("Gompertz cens[%d] neither 1 nor 0\n",i);
/*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
- if (wav[i] > 1 ) { /* ??? */
- L=L+A*weight[i];
+ L=L+A*weight[i];
/* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/
- }
- }
+ }
+ }
- /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/
+ /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/
return -2*L*num/sump;
}
/******************* Gompertz_f Likelihood ******************************/
double gompertz_f(const gsl_vector *v, void *params)
{
- double A,B,LL=0.0,sump=0.,num=0.;
+ double A=0.,B=0.,LL=0.0,sump=0.,num=0.;
double *x= (double *) v->data;
int i,n=0; /* n is the size of the sample */
int i=0, j=0, n=0, iv=0, v;
int lstra;
int linei, month, year,iout;
+ int noffset=0; /* This is the offset if BOM data file */
char line[MAXLINE], linetmp[MAXLINE];
char stra[MAXLINE], strb[MAXLINE];
char *stratrunc;
fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1;
}
- i=1;
+ /* Is it a BOM UTF-8 Windows file? */
+ /* First data line */
linei=0;
+ while(fgets(line, MAXLINE, fic)) {
+ noffset=0;
+ if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */
+ {
+ noffset=noffset+3;
+ printf("# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile);fflush(stdout);
+ fprintf(ficlog,"# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile);
+ fflush(ficlog); return 1;
+ }
+ /* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/
+ else if( line[0] == (char)0xFF && line[1] == (char)0xFE)
+ {
+ noffset=noffset+2;
+ printf("# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout);
+ fprintf(ficlog,"# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);
+ fflush(ficlog); return 1;
+ }
+ else if( line[0] == 0 && line[1] == 0)
+ {
+ if( line[2] == (char)0xFE && line[3] == (char)0xFF){
+ noffset=noffset+4;
+ printf("# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout);
+ fprintf(ficlog,"# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);
+ fflush(ficlog); return 1;
+ }
+ } else{
+ ;/*printf(" Not a BOM file\n");*/
+ }
+ /* If line starts with a # it is a comment */
+ if (line[noffset] == '#') {
+ linei=linei+1;
+ break;
+ }else{
+ break;
+ }
+ }
+ fclose(fic);
+ if((fic=fopen(datafile,"r"))==NULL) {
+ printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
+ fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1;
+ }
+ /* Not a Bom file */
+
+ i=1;
while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
linei=linei+1;
for(j=strlen(line); j>=0;j--){ /* Untabifies line */
return 1;
}
anint[j][i]= (double) year;
- mint[j][i]= (double)month;
+ mint[j][i]= (double)month;
+ /* if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){ */
+ /* printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */
+ /* fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */
+ /* } */
strcpy(line,stra);
} /* End loop on waves */
}
annais[i]=(double)(year);
- moisnais[i]=(double)(month);
+ moisnais[i]=(double)(month);
+ for (j=1;j<=maxwav;j++){
+ if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){
+ printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j,(int)moisnais[i],(int)annais[i]);
+ fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j, (int)moisnais[i],(int)annais[i]);
+ }
+ }
+
strcpy(line,stra);
/* Sample weight */
noffset=noffset+3;
printf("# File is an UTF8 Bom.\n"); // 0xBF
}
- else if( line[0] == (char)0xFE && line[1] == (char)0xFF)
+/* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/
+ else if( line[0] == (char)0xFF && line[1] == (char)0xFE)
{
noffset=noffset+2;
printf("# File is an UTF16BE BOM file\n");
}
if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
if (num_filled != 1){
- printf("ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line);
- fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line);
+ printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
+ fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
model[0]='\0';
goto end;
}
ximort[i][j]=(i == j ? 1.0 : 0.0);
}
- /*p[1]=0.0268; p[NDIM]=0.083;*/
- /*printf("%lf %lf", p[1], p[2]);*/
+ p[1]=0.0268; p[NDIM]=0.083;
+ /* printf("%lf %lf", p[1], p[2]); */
#ifdef GSL
printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
}
- lsurv=vector(1,AGESUP);
- lpop=vector(1,AGESUP);
- tpop=vector(1,AGESUP);
+ lsurv=vector(agegomp,AGESUP);
+ lpop=vector(agegomp,AGESUP);
+ tpop=vector(agegomp,AGESUP);
lsurv[agegomp]=100000;
for (k=agegomp;k<=AGESUP;k++) {
stepm, weightopt,\
model,imx,p,matcov,agemortsup);
- free_vector(lsurv,1,AGESUP);
- free_vector(lpop,1,AGESUP);
- free_vector(tpop,1,AGESUP);
+ free_vector(lsurv,agegomp,AGESUP);
+ free_vector(lpop,agegomp,AGESUP);
+ free_vector(tpop,agegomp,AGESUP);
free_matrix(ximort,1,NDIM,1,NDIM);
free_ivector(dcwave,firstobs,lastobs);
free_vector(agecens,firstobs,lastobs);
prvforecast = 1;
}
else if((num_filled=sscanf(line,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",&prevfcast,&yrfproj,&mobilavproj)) !=EOF){/* && (num_filled == 3))*/
- printf("prevforecast=%d yearsfproj=%lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
- fprintf(ficlog,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
- fprintf(ficres,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
+ printf("prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
+ fprintf(ficlog,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
+ fprintf(ficres,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
prvforecast = 2;
}
else {
prvbackcast = 1;
}
else if((num_filled=sscanf(line,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",&prevbcast,&yrbproj,&mobilavproj)) ==3){/* && (num_filled == 3))*/
- printf("prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
- fprintf(ficlog,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
- fprintf(ficres,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
+ printf("prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
+ fprintf(ficlog,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
+ fprintf(ficres,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
prvbackcast = 2;
}
else {
git://henry.ined.fr / .git / commitdiff