--- imach/src/imach.c 2015/09/22 19:45:16 1.202 +++ imach/src/imach.c 2015/11/17 22:12:03 1.209 @@ -1,6 +1,38 @@ -/* $Id: imach.c,v 1.202 2015/09/22 19:45:16 brouard Exp $ +/* $Id: imach.c,v 1.209 2015/11/17 22:12:03 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.209 2015/11/17 22:12:03 brouard + Summary: Adding ftolpl parameter + Author: N Brouard + + We had difficulties to get smoothed confidence intervals. It was due + to the period prevalence which wasn't computed accurately. The inner + parameter ftolpl is now an outer parameter of the .imach parameter + file after estepm. If ftolpl is small 1.e-4 and estepm too, + computation are long. + + Revision 1.208 2015/11/17 14:31:57 brouard + Summary: temporary + + 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 *** + + Revision 1.205 2015/10/23 15:50:53 brouard + Summary: 0.98r3 some clarification for graphs on likelihood contributions + + Revision 1.204 2015/10/01 16:20:26 brouard + Summary: Some new graphs of contribution to likelihood + + Revision 1.203 2015/09/30 17:45:14 brouard + Summary: looking at better estimation of the hessian + + Also a better criteria for convergence to the period prevalence And + therefore adding the number of years needed to converge. (The + prevalence in any alive state shold sum to one + Revision 1.202 2015/09/22 19:45:16 brouard Summary: Adding some overall graph on contribution to likelihood. Might change @@ -647,7 +679,10 @@ /* #define DEBUG */ /* #define DEBUGBRENT */ -#define DEBUGLINMIN +/* #define DEBUGLINMIN */ +/* #define DEBUGHESS */ +#define DEBUGHESSIJ +/* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan)*\/ */ #define POWELL /* Instead of NLOPT */ #define POWELLF1F3 /* Skip test */ /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ @@ -735,12 +770,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.202 2015/09/22 19:45:16 brouard Exp $ */ +/* $Id: imach.c,v 1.209 2015/11/17 22:12:03 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; -char fullversion[]="$Revision: 1.202 $ $Date: 2015/09/22 19:45:16 $"; +char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; +char fullversion[]="$Revision: 1.209 $ $Date: 2015/11/17 22:12:03 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -930,7 +965,7 @@ static int split( char *path, char *dirc } /* got dirc from getcwd*/ printf(" DIRC = %s \n",dirc); - } else { /* strip direcotry from path */ + } else { /* strip directory from path */ ss++; /* after this, the filename */ l2 = strlen( ss ); /* length of filename */ if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); @@ -1602,8 +1637,11 @@ void linmin(double p[], double xi[], int double xx,xmin,bx,ax; double fx,fb,fa; - double scale=10., axs, xxs, xxss; /* Scale added for infinity */ - +#ifdef LINMINORIGINAL +#else + double scale=10., axs, xxs; /* Scale added for infinity */ +#endif + ncom=n; pcom=vector(1,n); xicom=vector(1,n); @@ -1613,12 +1651,15 @@ void linmin(double p[], double xi[], int xicom[j]=xi[j]; /* Former scale xi[j] of currrent direction i */ } - /* axs=0.0; */ - /* xxss=1; /\* 1 and using scale *\/ */ - xxs=1; - /* do{ */ - ax=0.; +#ifdef LINMINORIGINAL + xx=1.; +#else + axs=0.0; + xxs=1.; + do{ xx= xxs; +#endif + ax=0.; mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */ /* brackets with inputs ax=0 and xx=1, but points, pcom=p, and directions values, xicom=xi, are sent via f1dim(x) */ /* xt[x,j]=pcom[j]+x*xicom[j] f(ax) = f(xt(a,j=1,n)) = f(p(j) + 0 * xi(j)) and f(xx) = f(xt(x, j=1,n)) = f(p(j) + 1 * xi(j)) */ @@ -1626,12 +1667,19 @@ void linmin(double p[], double xi[], int /* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */ /* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */ /* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus [0:xi[j]]*/ - /* if (fx != fx){ */ - /* xxs=xxs/scale; /\* Trying a smaller xx, closer to initial ax=0 *\/ */ - /* printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); */ - /* } */ - /* }while(fx != fx); */ - +#ifdef LINMINORIGINAL +#else + if (fx != fx){ + xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */ + printf("|"); + fprintf(ficlog,"|"); +#ifdef DEBUGLINMIN + printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); +#endif + } + }while(fx != fx); +#endif + #ifdef DEBUGLINMIN printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); @@ -1650,14 +1698,23 @@ void linmin(double p[], double xi[], int fprintf(ficlog,"linmin end "); #endif for (j=1;j<=n;j++) { - /* printf(" before xi[%d]=%12.8f", j,xi[j]); */ - xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ - /* if(xxs <1.0) */ - /* printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); */ +#ifdef LINMINORIGINAL + xi[j] *= xmin; +#else +#ifdef DEBUGLINMIN + if(xxs <1.0) + printf(" before xi[%d]=%12.8f", j,xi[j]); +#endif + xi[j] *= xmin*xxs; /* xi rescaled by xmin and number of loops: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ +#ifdef DEBUGLINMIN + if(xxs <1.0) + printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); +#endif +#endif p[j] += xi[j]; /* Parameters values are updated accordingly */ } - /* printf("\n"); */ #ifdef DEBUGLINMIN + printf("\n"); printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); fprintf(ficlog,"Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); for (j=1;j<=n;j++) { @@ -1668,6 +1725,7 @@ void linmin(double p[], double xi[], int fprintf(ficlog,"\n"); } } +#else #endif free_vector(xicom,1,n); free_vector(pcom,1,n); @@ -1745,10 +1803,10 @@ void powell(double p[], double **xi, int for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */ fptt=(*fret); #ifdef DEBUG - printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); - fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); + printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); + fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); #endif - printf("%d",i);fflush(stdout); /* print direction (parameter) i */ + printf("%d",i);fflush(stdout); /* print direction (parameter) i */ fprintf(ficlog,"%d",i);fflush(ficlog); linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/ /* Outputs are fret(new point p) p is updated and xit rescaled */ @@ -1917,19 +1975,40 @@ void powell(double p[], double **xi, int /**** Prevalence limit (stable or period prevalence) ****************/ -double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) +double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij) { /* Computes the prevalence limit in each live state at age x by left multiplying the unit - matrix by transitions matrix until convergence is reached */ - + matrix by transitions matrix until convergence is reached with precision ftolpl */ + /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ + /* Wx is row vector: population in state 1, population in state 2, population dead */ + /* or prevalence in state 1, prevalence in state 2, 0 */ + /* newm is the matrix after multiplications, its rows are identical at a factor */ + /* Initial matrix pimij */ + /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ + /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ + /* 0, 0 , 1} */ + /* + * and after some iteration: */ + /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */ + /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */ + /* 0, 0 , 1} */ + /* And prevalence by suppressing the deaths are close to identical rows in prlim: */ + /* {0.51571254859325999, 0.4842874514067399, */ + /* 0.51326036147820708, 0.48673963852179264} */ + /* If we start from prlim again, prlim tends to a constant matrix */ + int i, ii,j,k; - double min, max, maxmin, maxmax,sumnew=0.; + double *min, *max, *meandiff, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ double **out, cov[NCOVMAX+1], **pmij(); double **newm; - double agefin, delaymax=100 ; /* Max number of years to converge */ - long int ncvyear=0, ncvloop=0; + double agefin, delaymax=200. ; /* 100 Max number of years to converge */ + int ncvloop=0; + min=vector(1,nlstate); + max=vector(1,nlstate); + meandiff=vector(1,nlstate); + for (ii=1;ii<=nlstate+ndeath;ii++) for (j=1;j<=nlstate+ndeath;j++){ oldm[ii][j]=(ii==j ? 1.0 : 0.0); @@ -1967,29 +2046,45 @@ double **prevalim(double **prlim, int nl savm=oldm; oldm=newm; - maxmax=0.; - for(j=1;j<=nlstate;j++){ - min=1.; - max=0.; - 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++){ + max[j]=0.; + min[j]=1.; + } + 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); - 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); */ + max[j]=FMAX(max[j],prlim[i][j]); + min[j]=FMIN(min[j],prlim[i][j]); } - maxmin=max-min; - maxmax=FMAX(maxmax,maxmin); + } + + maxmax=0.; + for(j=1; j<=nlstate; j++){ + meandiff[j]=(max[j]-min[j])/(max[j]+min[j])*2.; /* mean difference for each column */ + maxmax=FMAX(maxmax,meandiff[j]); + /* printf(" age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, j, meandiff[j],(int)agefin, j, max[j], j, min[j],maxmax); */ } /* 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); */ if(maxmax < ftolpl){ - /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age-(int)agefin); */ + /* printf("maxmax=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */ + free_vector(min,1,nlstate); + free_vector(max,1,nlstate); + free_vector(meandiff,1,nlstate); return prlim; } } /* age loop */ - printf("Warning: the stable prevalence did not converge with the required precision ftolpl=6*10^5*ftol=%g. \n\ -Earliest age to start was %d-%d=%d, ncvloop=%ld, ncvyear=%d\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); + /* 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. Try to lower 'ftolpl'. \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); + /* 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); */ + free_vector(min,1,nlstate); + free_vector(max,1,nlstate); + free_vector(meandiff,1,nlstate); + return prlim; /* should not reach here */ } @@ -2565,9 +2660,9 @@ double funcone( double *x) ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ if(globpr){ - fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f\ + fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ %11.6f %11.6f %11.6f ", \ - num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], + num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); for(k=1,llt=0.,l=0.; k<=nlstate; k++){ llt +=ll[k]*gipmx/gsw; @@ -2605,8 +2700,8 @@ void likelione(FILE *ficres,double p[], 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_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 age i s1 s2 mi mw dh likeli weight 2wlli out sav "); + fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -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 age i s1 s2 mi mw dh likeli weight %%weight 2wlli 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," -2*gipw/gsw*weight*ll[%d]++",k); @@ -2616,11 +2711,23 @@ void likelione(FILE *ficres,double p[], *fretone=(*funcone)(p); if(*globpri !=0){ fclose(ficresilk); - fprintf(fichtm,"\n
File of contributions to the likelihood computed with initial parameters and mle >= 1. You should at least run with mle >= 1 and starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); - fprintf(fichtm,"
- The first 3 individuals are drawn with lines. The function drawn is -2Log(L) in log scale: %s.png
\ -",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + if (mle ==0) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with initial parameters and mle = %d.",mle); + else if(mle >=1) + fprintf(fichtm,"\n
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: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); + + + for (k=1; k<= nlstate ; k++) { + fprintf(fichtm,"
- Probability p%dj by origin %d and destination j %s-p%dj.png
\ +",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); + } + fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + fprintf(fichtm,"
- and by state of destination %s-dest.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); fflush(fichtm); - } + } return; } @@ -2694,32 +2801,32 @@ void mlikeli(FILE *ficres,double p[], in #endif free_matrix(xi,1,npar,1,npar); fclose(ficrespow); - printf("#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); - fprintf(ficlog,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); } /**** Computes Hessian and covariance matrix ***/ -void hesscov(double **matcov, double p[], int npar, double delti[], double ftolhess, double (*func)(double [])) +void hesscov(double **matcov, double **hess, double p[], int npar, double delti[], double ftolhess, double (*func)(double [])) { double **a,**y,*x,pd; - double **hess; + /* double **hess; */ int i, j; int *indx; double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar); - double hessij(double p[], double delti[], int i, int j,double (*func)(double []),int npar); + double hessij(double p[], double **hess, double delti[], int i, int j,double (*func)(double []),int npar); void lubksb(double **a, int npar, int *indx, double b[]) ; void ludcmp(double **a, int npar, int *indx, double *d) ; double gompertz(double p[]); - hess=matrix(1,npar,1,npar); + /* hess=matrix(1,npar,1,npar); */ printf("\nCalculation of the hessian matrix. Wait...\n"); fprintf(ficlog,"\nCalculation of the hessian matrix. Wait...\n"); for (i=1;i<=npar;i++){ - printf("%d",i);fflush(stdout); - fprintf(ficlog,"%d",i);fflush(ficlog); + printf("%d-",i);fflush(stdout); + fprintf(ficlog,"%d-",i);fflush(ficlog); hess[i][i]=hessii(p,ftolhess,i,delti,func,npar); @@ -2730,9 +2837,9 @@ void hesscov(double **matcov, double p[] for (i=1;i<=npar;i++) { for (j=1;j<=npar;j++) { if (j>i) { - printf(".%d%d",i,j);fflush(stdout); - fprintf(ficlog,".%d%d",i,j);fflush(ficlog); - hess[i][j]=hessij(p,delti,i,j,func,npar); + printf(".%d-%d",i,j);fflush(stdout); + fprintf(ficlog,".%d-%d",i,j);fflush(ficlog); + hess[i][j]=hessij(p,hess, delti,i,j,func,npar); hess[j][i]=hess[i][j]; /*printf(" %lf ",hess[i][j]);*/ @@ -2766,53 +2873,78 @@ void hesscov(double **matcov, double p[] fprintf(ficlog,"\n#Hessian matrix#\n"); for (i=1;i<=npar;i++) { for (j=1;j<=npar;j++) { - printf("%.3e ",hess[i][j]); - fprintf(ficlog,"%.3e ",hess[i][j]); + printf("%.6e ",hess[i][j]); + fprintf(ficlog,"%.6e ",hess[i][j]); } printf("\n"); fprintf(ficlog,"\n"); } + /* printf("\n#Covariance matrix#\n"); */ + /* fprintf(ficlog,"\n#Covariance matrix#\n"); */ + /* for (i=1;i<=npar;i++) { */ + /* for (j=1;j<=npar;j++) { */ + /* printf("%.6e ",matcov[i][j]); */ + /* fprintf(ficlog,"%.6e ",matcov[i][j]); */ + /* } */ + /* printf("\n"); */ + /* fprintf(ficlog,"\n"); */ + /* } */ + /* Recompute Inverse */ - for (i=1;i<=npar;i++) - for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; - ludcmp(a,npar,indx,&pd); + /* for (i=1;i<=npar;i++) */ + /* for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; */ + /* ludcmp(a,npar,indx,&pd); */ + + /* printf("\n#Hessian matrix recomputed#\n"); */ + + /* for (j=1;j<=npar;j++) { */ + /* for (i=1;i<=npar;i++) x[i]=0; */ + /* x[j]=1; */ + /* lubksb(a,npar,indx,x); */ + /* for (i=1;i<=npar;i++){ */ + /* y[i][j]=x[i]; */ + /* printf("%.3e ",y[i][j]); */ + /* fprintf(ficlog,"%.3e ",y[i][j]); */ + /* } */ + /* printf("\n"); */ + /* fprintf(ficlog,"\n"); */ + /* } */ + + /* Verifying the inverse matrix */ +#ifdef DEBUGHESS + y=matprod2(y,hess,1,npar,1,npar,1,npar,matcov); - /* printf("\n#Hessian matrix recomputed#\n"); + printf("\n#Verification: multiplying the matrix of covariance by the Hessian matrix, should be unity:#\n"); + fprintf(ficlog,"\n#Verification: multiplying the matrix of covariance by the Hessian matrix. Should be unity:#\n"); for (j=1;j<=npar;j++) { - for (i=1;i<=npar;i++) x[i]=0; - x[j]=1; - lubksb(a,npar,indx,x); for (i=1;i<=npar;i++){ - y[i][j]=x[i]; - printf("%.3e ",y[i][j]); - fprintf(ficlog,"%.3e ",y[i][j]); + printf("%.2f ",y[i][j]); + fprintf(ficlog,"%.2f ",y[i][j]); } printf("\n"); fprintf(ficlog,"\n"); } - */ +#endif free_matrix(a,1,npar,1,npar); free_matrix(y,1,npar,1,npar); free_vector(x,1,npar); free_ivector(indx,1,npar); - free_matrix(hess,1,npar,1,npar); + /* free_matrix(hess,1,npar,1,npar); */ } /*************** hessian matrix ****************/ double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar) -{ +{ /* Around values of x, computes the function func and returns the scales delti and hessian */ int i; int l=1, lmax=20; - double k1,k2; + double k1,k2, res, fx; double p2[MAXPARM+1]; /* identical to x */ - double res; double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; - double fx; int k=0,kmax=10; double l1; @@ -2828,9 +2960,9 @@ double hessii(double x[], double delta, p2[theta]=x[theta]-delt; k2=func(p2)-fx; /*res= (k1-2.0*fx+k2)/delt/delt; */ - res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ + res= (k1+k2)/delt/delt/2.; /* Divided by 2 because L and not 2*L */ -#ifdef DEBUGHESS +#ifdef DEBUGHESSII printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); #endif @@ -2844,48 +2976,125 @@ double hessii(double x[], double delta, else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ delts=delt; } - } + } /* End loop k */ } delti[theta]=delts; return res; } -double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) +double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,double (*func)(double []),int npar) { int i; int l=1, lmax=20; double k1,k2,k3,k4,res,fx; double p2[MAXPARM+1]; - int k; + int k, kmax=1; + double v1, v2, cv12, lc1, lc2; + int firstime=0; + fx=func(x); - for (k=1; k<=2; k++) { + for (k=1; k<=kmax; k=k+10) { for (i=1;i<=npar;i++) p2[i]=x[i]; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; + p2[thetai]=x[thetai]+delti[thetai]*k; + p2[thetaj]=x[thetaj]+delti[thetaj]*k; k1=func(p2)-fx; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; + p2[thetai]=x[thetai]+delti[thetai]*k; + p2[thetaj]=x[thetaj]-delti[thetaj]*k; k2=func(p2)-fx; - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; + p2[thetai]=x[thetai]-delti[thetai]*k; + p2[thetaj]=x[thetaj]+delti[thetaj]*k; k3=func(p2)-fx; - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; + p2[thetai]=x[thetai]-delti[thetai]*k; + p2[thetaj]=x[thetaj]-delti[thetaj]*k; k4=func(p2)-fx; - res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */ -#ifdef DEBUG - 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); + 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 || 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]; + v2=hess[thetaj][thetaj]; + cv12=res; + /* Computing eigen value of Hessian matrix */ + lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + if ((lc2 <0) || (lc1 <0) ){ + printf("Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj); + fprintf(ficlog, "Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj); + 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); + } #endif } return res; } + /* Not done yet: Was supposed to fix if not exactly at the maximum */ +/* double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) */ +/* { */ +/* int i; */ +/* int l=1, lmax=20; */ +/* double k1,k2,k3,k4,res,fx; */ +/* double p2[MAXPARM+1]; */ +/* double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; */ +/* int k=0,kmax=10; */ +/* double l1; */ + +/* fx=func(x); */ +/* for(l=0 ; l <=lmax; l++){ /\* Enlarging the zone around the Maximum *\/ */ +/* l1=pow(10,l); */ +/* delts=delt; */ +/* for(k=1 ; k khi/nkhif) || (k2 >khi/nkhif) || (k4 >khi/nkhif) || (k4 >khi/nkhif)){ /\* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. *\/ */ +/* k=kmax; l=lmax*10; */ +/* } */ +/* else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ */ +/* delts=delt; */ +/* } */ +/* } /\* End loop k *\/ */ +/* } */ +/* delti[theta]=delts; */ +/* return res; */ +/* } */ + + /************** Inverse of matrix **************/ void ludcmp(double **a, int n, int *indx, double *d) { @@ -3818,7 +4027,7 @@ void cvevsij(double ***eij, double x[], } /************ Variance ******************/ -void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) + void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) { /* Variance of health expectancies */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ @@ -3872,7 +4081,6 @@ void varevsij(char optionfilefiname[], d 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); @@ -3883,6 +4091,7 @@ void varevsij(char optionfilefiname[], d 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);*/ @@ -3920,9 +4129,9 @@ void varevsij(char optionfilefiname[], d /* 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 because of memory size limitations, + 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 @@ -3943,8 +4152,8 @@ void varevsij(char optionfilefiname[], d for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); if (popbased==1) { if(mobilav ==0){ @@ -3956,13 +4165,14 @@ void varevsij(char optionfilefiname[], d } } + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Returns p3mat[i][j][h] for h=1 to nhstepm */ for(j=1; j<= nlstate; j++){ for(h=0; h<=nhstepm; h++){ for(i=1, gp[h][j]=0.;i<=nlstate;i++) gp[h][j] += prlim[i][i]*p3mat[i][j][h]; } } - /* This for computing probability of death (h=1 means + /* Next for computing probability of death (h=1 means computed over hstepm matrices product = hstepm*stepm months) as a weighted average of prlim. */ @@ -3974,8 +4184,8 @@ void varevsij(char optionfilefiname[], d for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij); if (popbased==1) { if(mobilav ==0){ @@ -3987,6 +4197,8 @@ void varevsij(char optionfilefiname[], d } } + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ for(h=0; h<=nhstepm; h++){ for(i=1, gm[h][j]=0.;i<=nlstate;i++) @@ -4049,8 +4261,8 @@ void varevsij(char optionfilefiname[], d varppt[j][i]=doldmp[j][i]; /* end ppptj */ /* x centered again */ - hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); + + prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij); if (popbased==1) { if(mobilav ==0){ @@ -4066,6 +4278,7 @@ void varevsij(char optionfilefiname[], d computed over hstepm (estepm) matrices product = hstepm*stepm months) as a weighted average of prlim. */ + hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); for(j=nlstate+1;j<=nlstate+ndeath;j++){ for(i=1,gmp[j]=0.;i<= nlstate; i++) gmp[j] += prlim[i][i]*p3mat[i][j][1]; @@ -4128,9 +4341,9 @@ void varevsij(char optionfilefiname[], d } /* end varevsij */ /************ Variance of prevlim ******************/ -void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, char strstart[]) + void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[]) { - /* Variance of prevalence limit */ + /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ double **dnewm,**doldm; @@ -4138,6 +4351,7 @@ void varprevlim(char fileres[], double * double *xp; double *gp, *gm; double **gradg, **trgradg; + double **mgm, **mgp; double age,agelim; int theta; @@ -4160,6 +4374,8 @@ void varprevlim(char fileres[], double * if (stepm >= YEARM) hstepm=1; nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ 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); @@ -4167,18 +4383,27 @@ void varprevlim(char fileres[], double * 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,ij); - for(i=1;i<=nlstate;i++) + if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + else + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + 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,ij); - for(i=1;i<=nlstate;i++) + if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + else + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + 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]; + /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */ } /* End theta */ trgradg =matrix(1,nlstate,1,npar); @@ -4186,11 +4411,34 @@ void varprevlim(char fileres[], double * for(j=1; j<=nlstate;j++) for(theta=1; theta <=npar; theta++) trgradg[j][theta]=gradg[theta][j]; + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* 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==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* 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==79 ||(int)age== 80 ||(int)age== 81){ matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + }else{ + matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + } for(i=1;i<=nlstate;i++) varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ @@ -4200,6 +4448,8 @@ void varprevlim(char fileres[], double * 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 */ @@ -4611,7 +4861,7 @@ divided by h: hPij/h : \n- Survival functions from state %d in any different live states and total.\ + fprintf(fichtm,"
\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.\ %s%d_%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1); } @@ -4621,7 +4871,7 @@ divided by h: hPij/h : ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s%d%d.svg
\ + fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s_%d%d.svg
\ ",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); } /* } /\* end i1 *\/ */ @@ -4631,7 +4881,7 @@ divided by h: hPij/h :

  • Result files (second order: variances)

    \n\ - Parameter file with estimated parameters and covariance matrix: %s
    \ - - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.
    \ + - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file if optimization has been done (mle != 0).
    \ But because parameters are usually highly correlated (a higher incidence of disability \ and a higher incidence of recovery can give very close observed transition) it might \ be very useful to look not only at linear confidence intervals estimated from the \ @@ -4691,15 +4941,15 @@ See page 'Matrix of variance-covariance } for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"
    - Observed (cross-sectional) and period (incidence based) \ -prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg
    \ -",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1); +prevalence (with 95%% confidence interval) in state (%d): %s_%d-%d.svg
    \ +",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1); } fprintf(fichtm,"\n
    - Total life expectancy by age and \ health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \ true period expectancies (those weighted with period prevalences are also\ drawn in addition to the population based expectancies computed using\ - observed and cahotic prevalences: %s_%d.svg
    \ -",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); + observed and cahotic prevalences:     %s_%d.svg
    \ +",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); /* } /\* end i1 *\/ */ }/* End k1 */ fprintf(fichtm,""); @@ -4728,15 +4978,28 @@ void printinggnuplot(char fileresu[], ch fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";"); /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */ - fprintf(ficgp,"\nset ter png size 640, 480"); -/* good for mle=4 plot by number of matrix products. + fprintf(ficgp,"\nset ter pngcairo size 640, 480"); +/* nice for mle=4 plot by number of matrix products. replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */ /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */ /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */ - fprintf(ficgp,"\nset out \"%s.png\";",subdirf2(optionfilefiname,"ILK_")); - fprintf(ficgp,"\nplot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); - fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); - fprintf(ficgp,"\nset out\n"); + fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_")); + fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk)); + fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_")); + fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk)); + for (i=1; i<= nlstate ; i ++) { + fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i); + fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk)); + fprintf(ficgp," u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1); + for (j=2; j<= nlstate+ndeath ; j ++) { + fprintf(ficgp,",\\\n \"\" u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j); + } + fprintf(ficgp,";\nset out; unset ylabel;\n"); + } + /* unset log; plot "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */ + /* fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */ + /* fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */ + fprintf(ficgp,"\nset out;unset log\n"); /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */ strcpy(dirfileres,optionfilefiname); @@ -6317,11 +6580,12 @@ void syscompilerinfo(int logged) } - int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl){ + int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){ /*--------------- Prevalence limit (period or stable prevalence) --------------*/ int i, j, k, i1 ; /* double ftolpl = 1.e-10; */ double age, agebase, agelim; + double tot; strcpy(filerespl,"PL_"); strcat(filerespl,fileresu); @@ -6332,7 +6596,7 @@ void syscompilerinfo(int logged) printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl); fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl); pstamp(ficrespl); - fprintf(ficrespl,"# Period (stable) prevalence \n"); + fprintf(ficrespl,"# Period (stable) prevalence. Precision given by ftolpl=%g \n", ftolpl); fprintf(ficrespl,"#Age "); for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); fprintf(ficrespl,"\n"); @@ -6367,18 +6631,21 @@ void syscompilerinfo(int logged) for(j=1;j<=cptcoveff;j++) { fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } - for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); - fprintf(ficrespl,"\n"); + for(i=1; i<=nlstate;i++) fprintf(ficrespl," %d-%d ",i,i); + fprintf(ficrespl,"Total Years_to_converge\n"); for (age=agebase; age<=agelim; age++){ /* for (age=agebase; age<=agebase; age++){ */ - prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); + prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k); fprintf(ficrespl,"%.0f ",age ); for(j=1;j<=cptcoveff;j++) fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - for(i=1; i<=nlstate;i++) + tot=0.; + for(i=1; i<=nlstate;i++){ + tot += prlim[i][i]; fprintf(ficrespl," %.5f", prlim[i][i]); - fprintf(ficrespl,"\n"); + } + fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp); } /* Age */ /* was end of cptcod */ } /* cptcov */ @@ -6471,7 +6738,7 @@ int main(int argc, char *argv[]) #endif int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); int i,j, k, n=MAXN,iter=0,m,size=100, cptcod; - + int ncvyear=0; /* Number of years needed for the period prevalence to converge */ int jj, ll, li, lj, lk; int numlinepar=0; /* Current linenumber of parameter file */ int num_filled; @@ -6519,6 +6786,7 @@ int main(int argc, char *argv[]) double ***param; /* Matrix of parameters */ double *p; double **matcov; /* Matrix of covariance */ + double **hess; /* Hessian matrix */ double ***delti3; /* Scale */ double *delti; /* Scale */ double ***eij, ***vareij; @@ -6578,14 +6846,23 @@ int main(int argc, char *argv[]) printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); - fgets(pathr,FILENAMELENGTH,stdin); + if(!fgets(pathr,FILENAMELENGTH,stdin)){ + printf("ERROR Empty parameter file name\n"); + goto end; + } i=strlen(pathr); if(pathr[i-1]=='\n') pathr[i-1]='\0'; i=strlen(pathr); - if(pathr[i-1]==' ') /* This may happen when dragging on oS/X! */ + if(i >= 1 && pathr[i-1]==' ') {/* This may happen when dragging on oS/X! */ pathr[i-1]='\0'; - for (tok = pathr; tok != NULL; ){ + } + i=strlen(pathr); + if( i==0 ){ + printf("ERROR Empty parameter file name\n"); + goto end; + } + for (tok = pathr; tok != NULL; ){ printf("Pathr |%s|\n",pathr); while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); printf("val= |%s| pathr=%s\n",val,pathr); @@ -6721,11 +6998,13 @@ int main(int argc, char *argv[]) if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \ &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ if (num_filled != 8) { - printf("Not 8\n"); + printf("Not 8 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); + printf("but line=%s\n",line); } printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt); } - ftolpl=6*ftol*1.e5; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ + /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ + /*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ /* Third parameter line */ while(fgets(line, MAXLINE, ficpar)) { /* If line starts with a # it is a comment */ @@ -6755,14 +7034,13 @@ int main(int argc, char *argv[]) } } /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */ + printf("model=1+age+%s\n",model);fflush(stdout); } /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */ /* numlinepar=numlinepar+3; /\* In general *\/ */ /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */ - if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */ - model[strlen(model)-1]='\0'; - fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); - fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); + fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); + fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); fflush(ficlog); /* if(model[0]=='#'|| model[0]== '\0'){ */ if(model[0]=='#'){ @@ -6827,6 +7105,7 @@ int main(int argc, char *argv[]) fprintf(ficlog," You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso); param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); matcov=matrix(1,npar,1,npar); + hess=matrix(1,npar,1,npar); } else{ /* Read guessed parameters */ @@ -6935,6 +7214,7 @@ run imach with mle=-1 to get a correct t ungetc(c,ficpar); matcov=matrix(1,npar,1,npar); + hess=matrix(1,npar,1,npar); for(i=1; i <=npar; i++) for(j=1; j <=npar; j++) matcov[i][j]=0.; @@ -7177,7 +7457,7 @@ Please run with mle=-1 to get a correct printf("Problem with file %s",optionfilegnuplot); } else{ - fprintf(ficgp,"\n# %s\n", version); + fprintf(ficgp,"\n# IMaCh-%s\n", version); fprintf(ficgp,"# %s\n", optionfilegnuplot); //fprintf(ficgp,"set missing 'NaNq'\n"); fprintf(ficgp,"set datafile missing 'NaNq'\n"); @@ -7204,13 +7484,15 @@ Please run with mle=-1 to get a correct else{ fprintf(fichtmcov,"\nIMaCh Cov %s\n %s
    %s     \
    \n\ -Title=%s
    Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
    \n",\ +Title=%s
    Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
    \n",\ optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); } - fprintf(fichtm,"\nIMaCh %s\n %s
    %s     \ + fprintf(fichtm,"\n\n\nIMaCh %s\n     IMaCh for Interpolated Markov Chain
    \nSponsored by Copyright (C) 2002-2015 INED-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015
    \ +
    \n\ +IMaCh-%s
    %s     \
    \n\ -Title=%s
    Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
    \n\ +Title=%s
    Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
    \n\ \n\
    \