--- imach/src/imach.c 2024/04/24 21:21:17 1.359 +++ imach/src/imach.c 2024/06/28 09:31:55 1.363 @@ -1,6 +1,29 @@ -/* $Id: imach.c,v 1.359 2024/04/24 21:21:17 brouard Exp $ +/* $Id: imach.c,v 1.363 2024/06/28 09:31:55 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.363 2024/06/28 09:31:55 brouard + Summary: Adding log lines too + + Revision 1.362 2024/06/28 08:00:31 brouard + Summary: 0.99s6 + + * imach.c (Module): s6 errors with age*age (harmless). + + Revision 1.361 2024/05/12 20:29:32 brouard + Summary: Version 0.99s5 + + * src/imach.c Version 0.99s5 In fact, the covariance of total life + expectancy e.. with a partial life expectancy e.j is high, + therefore the complete matrix of variance covariance has to be + included in the formula of the standard error of the proportion of + total life expectancy spent in a specific state: + var(X/Y)=mu_x^2/mu_y^2*(sigma_x^2/mu_x^2 -2 + sigma_xy/mu_x/mu_y+sigma^2/mu_y^2). Also an error with mle=-3 + made the program core dump. It is fixed in this version. + + Revision 1.360 2024/04/30 10:59:22 brouard + Summary: Version 0.99s4 and estimation of std of e.j/e.. + Revision 1.359 2024/04/24 21:21:17 brouard Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes @@ -1390,12 +1413,12 @@ double gnuplotversion=GNUPLOTVERSION; #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.359 2024/04/24 21:21:17 brouard Exp $ */ +/* $Id: imach.c,v 1.363 2024/06/28 09:31:55 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; -char fullversion[]="$Revision: 1.359 $ $Date: 2024/04/24 21:21:17 $"; +char copyright[]="April 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024"; +char fullversion[]="$Revision: 1.363 $ $Date: 2024/06/28 09:31:55 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1516,6 +1539,12 @@ char *endptr; long lval; double dval; +/* This for praxis gegen */ + /* int prin=1; */ + double h0=0.25; + double macheps; + double ffmin; + #define NR_END 1 #define FREE_ARG char* #define FTOL 1.0e-10 @@ -2992,24 +3021,30 @@ static void print2() /* print a line of /* printf("... after %u function calls ...\n", nf); */ /* printf("... including %u linear searches ...\n", nl); */ /* printf("%10d %10d%14.7g",nl, nf, fx); */ - printf ( "\n" ); + /* printf ( "\n" ); */ printf ( " Linear searches %d", nl ); /* printf ( " Linear searches %d\n", nl ); */ /* printf ( " Function evaluations %d\n", nf ); */ /* printf ( " Function value FX = %g\n", fx ); */ printf ( " Function evaluations %d", nf ); printf ( " Function value FX = %.12lf\n", fx ); + fprintf (ficlog, " Function evaluations %d", nf ); + fprintf (ficlog, " Function value FX = %.12lf\n", fx ); #ifdef DEBUGPRAX printf("n=%d prin=%d\n",n,prin); #endif - if(fx <= fmin) printf(" UNDEFINED "); else printf("%14.7g",log(fx-fmin)); + /* if(fx <= fmin) printf(" UNDEFINED "); else printf("%14.7g",log(fx-fmin)); */ if ( n <= 4 || 2 < prin ) { /* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */ - for(i=1;i<=n;i++)printf("%14.7g",x[i]); + for(i=1;i<=n;i++){ + printf("%14.7g",x[i]); + fprintf(ficlog,"%14.7g",x[i]); + } /* r8vec_print ( n, x, " X:" ); */ } printf("\n"); + fprintf(ficlog,"\n"); } @@ -3216,7 +3251,7 @@ L1: /* L1 or try loop */ if (k > 0) *d2 = 0; } #ifdef DEBUGPRAX - printf(" bebe end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2); + printf(" bebe end of min x1 might be very wrong x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2); #endif if (*d2 <= small_windows) *d2 = small_windows; *x1 = x2; fx = fm; @@ -3692,7 +3727,7 @@ mloop: printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t); #endif /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */ - minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global */ + minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global it seems that fx doesn't correspond to f(s=*x1) */ #ifdef DEBUGPRAX printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx); #endif @@ -4203,7 +4238,7 @@ void powell(double p[], double **xi, int printf(" + age*age "); fprintf(ficlog," + age*age "); } - for(j=1;j <=ncovmodel-2;j++){ + for(j=1;j <=ncovmodel-2-nagesqr;j++){ if(Typevar[j]==0) { printf(" + V%d ",Tvar[j]); fprintf(ficlog," + V%d ",Tvar[j]); @@ -4411,14 +4446,14 @@ void powell(double p[], double **xi, int #endif #ifdef POWELLORIGINAL if (t < 0.0) { /* Then we use it for new direction */ -#else +#else /* Not POWELLOriginal but Brouard's */ if (directest*t < 0.0) { /* Contradiction between both tests */ printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); } - if (directest < 0.0) { /* Then we use it for new direction */ + if (directest < 0.0) { /* Then we use (P0, Pn) for new direction Xi_n or Xi_iBig */ #endif #ifdef DEBUGLINMIN printf("Before linmin in direction P%d-P0\n",n); @@ -4452,6 +4487,21 @@ void powell(double p[], double **xi, int xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ } + +/* #else */ +/* for (i=1;i<=n-1;i++) { */ +/* for (j=1;j<=n;j++) { */ +/* xi[j][i]=xi[j][i+1]; /\* Standard method of conjugate directions, not Powell who changes the nth direction by p0 pn . *\/ */ +/* } */ +/* } */ +/* for (j=1;j<=n;j++) { */ +/* xi[j][n]=xit[j]; /\* and this nth direction by the by the average p_0 p_n *\/ */ +/* } */ +/* /\* for (j=1;j<=n-1;j++) { *\/ */ +/* /\* xi[j][1]=xi[j][j+1]; /\\* Standard method of conjugate directions *\\/ *\/ */ +/* /\* xi[j][n]=xit[j]; /\\* and this nth direction by the by the average p_0 p_n *\\/ *\/ */ +/* /\* } *\/ */ +/* #endif */ #ifdef LINMINORIGINAL #else for (j=1, flatd=0;j<=n;j++) { @@ -4476,8 +4526,8 @@ void powell(double p[], double **xi, int free_vector(pt,1,n); return; #endif - } -#endif + } /* endif(flatd >0) */ +#endif /* LINMINORIGINAL */ printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); @@ -6563,10 +6613,10 @@ void mlikeli(FILE *ficres,double p[], in #else /* FLATSUP */ /* powell(p,xi,npar,ftol,&iter,&fret,func);*/ /* praxis ( t0, h0, n, prin, x, beale_f ); */ - int prin=1; - double h0=0.25; - double macheps; - double fmin; + /* int prin=1; */ + /* double h0=0.25; */ + /* double macheps; */ + /* double fmin; */ macheps=pow(16.0,-13.0); /* #include "praxis.h" */ /* Be careful that praxis start at x[0] and powell start at p[1] */ @@ -6576,7 +6626,7 @@ printf("Praxis Gegenfurtner \n"); fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog); /* praxis ( ftol, h0, npar, prin, p1, func ); */ /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */ - fmin = praxis(ftol,macheps, h0, npar, prin, p, func); + ffmin = praxis(ftol,macheps, h0, npar, prin, p, func); printf("End Praxis\n"); #endif /* FLATSUP */ @@ -8553,7 +8603,9 @@ void concatwav(int wav[], int **dh, int /************ 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 *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres) { - /** Variance of health expectancies + /** Computes the matrix of variance covariance of health expectancies e.j= sum_i w_i e_ij where w_i depends of popbased, + * either cross-sectional or implied. + * return vareij[i][j][(int)age]=cov(e.i,e.j)=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20 * double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl); * double **newm; * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) @@ -8570,7 +8622,7 @@ void concatwav(int wav[], int **dh, int double ***gradg, ***trgradg; /**< for var eij */ double **gradgp, **trgradgp; /**< for var p point j */ double *gpp, *gmp; /**< for var p point j */ - double **varppt; /**< for var p point j nlstate to nlstate+ndeath */ + double **varppt; /**< for var p.3 p.death nlstate+1 to nlstate+ndeath */ double ***p3mat; double age,agelim, hf; /* double ***mobaverage; */ @@ -8638,7 +8690,7 @@ void concatwav(int wav[], int **dh, int fprintf(fichtm,"\n

  • Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); fprintf(fichtm,"\n
    %s
    \n",digitp); - varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); + varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); /* In fact, currently a double */ pstamp(ficresvij); fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); if(popbased==1) @@ -8707,7 +8759,7 @@ void concatwav(int wav[], int **dh, int prlim[i][i]=mobaverage[(int)age][i][ij]; } } - /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h. + /**< Computes the shifted plus (gp) transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h. */ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */ /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability @@ -8716,14 +8768,14 @@ void concatwav(int wav[], int **dh, int 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]; + gp[h][j] += prlim[i][i]*p3mat[i][j][h]; /* gp[h][j]= w_i h_pij */ } } /* Next for computing shifted+ probability of death (h=1 means computed over hstepm matrices product = hstepm*stepm months) as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 . */ - for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(j=nlstate+1;j<=nlstate+ndeath;j++){ /* Currently only once for theta plus p.3(age) Sum_i wi pi3*/ for(i=1,gpp[j]=0.; i<= nlstate; i++) gpp[j] += prlim[i][i]*p3mat[i][j][1]; } @@ -8745,9 +8797,9 @@ void concatwav(int wav[], int **dh, int } } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Still minus */ - for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ + for(j=1; j<= nlstate; j++){ /* gm[h][j]= Sum_i of wi * pij = h_p.j */ for(h=0; h<=nhstepm; h++){ for(i=1, gm[h][j]=0.;i<=nlstate;i++) gm[h][j] += prlim[i][i]*p3mat[i][j][h]; @@ -8755,37 +8807,39 @@ void concatwav(int wav[], int **dh, int } /* This for computing probability of death (h=1 means computed over hstepm matrices product = hstepm*stepm months) - as a weighted average of prlim. + as a weighted average of prlim. j is death. gmp[3]=sum_i w_i*p_i3=p.3 minus theta */ - for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(j=nlstate+1;j<=nlstate+ndeath;j++){ /* Currently only once theta_minus p.3=Sum_i wi pi3*/ for(i=1,gmp[j]=0.; i<= nlstate; i++) gmp[j] += prlim[i][i]*p3mat[i][j][1]; } /* end shifting computations */ - /**< Computing gradient matrix at horizon h + /**< Computing gradient of p.j matrix at horizon h and still for one parameter of vector theta + * equation 31 and 32 */ - for(j=1; j<= nlstate; j++) /* vareij */ + for(j=1; j<= nlstate; j++) /* computes grad p.j(x, over each h) where p.j is Sum_i w_i*pij(x over h) + * equation 24 */ for(h=0; h<=nhstepm; h++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; } - /**< Gradient of overall mortality p.3 (or p.j) + /**< Gradient of overall mortality p.3 (or p.death) */ - for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */ + for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* computes grad of p.3 from wi+pi3 grad p.3 (theta) */ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; } } /* End theta */ - /* We got the gradient matrix for each theta and state j */ - trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */ + /* We got the gradient matrix for each theta and each state j of gradg(h]theta][j)=grad(_hp.j(theta) */ + trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); - for(h=0; h<=nhstepm; h++) /* veij */ + for(h=0; h<=nhstepm; h++) /* veij */ /* computes the transposed of grad (_hp.j(theta)*/ for(j=1; j<=nlstate;j++) for(theta=1; theta <=npar; theta++) trgradg[h][j][theta]=gradg[h][theta][j]; - for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ + for(j=nlstate+1; j<=nlstate+ndeath;j++) /* computes transposed of grad p.3 (theta)*/ for(theta=1; theta <=npar; theta++) trgradgp[j][theta]=gradgp[theta][j]; /**< as well as its transposed matrix @@ -8797,8 +8851,11 @@ void concatwav(int wav[], int **dh, int vareij[i][j][(int)age] =0.; /* Computing trgradg by matcov by gradg at age and summing over h - * and k (nhstepm) formula 15 of article - * Lievre-Brouard-Heathcote + * and k (nhstepm) formula 32 of article + * Lievre-Brouard-Heathcote so that for each j, computes the cov(e.j,e.k) (formula 31). + * for given h and k computes trgradg[h](i,j) matcov (theta) gradg(k)(i,j) into vareij[i][j] which is + cov(e.i,e.j) and sums on h and k + * including the covariances. */ for(h=0;h<=nhstepm;h++){ @@ -8807,20 +8864,21 @@ void concatwav(int wav[], int **dh, int matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]); for(i=1;i<=nlstate;i++) for(j=1;j<=nlstate;j++) - vareij[i][j][(int)age] += doldm[i][j]*hf*hf; + vareij[i][j][(int)age] += doldm[i][j]*hf*hf; /* This is vareij=sum_h sum_k trgrad(h_pij) V(theta) grad(k_pij) + including the covariances of e.j */ } } - /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of - * p.j overall mortality formula 49 but computed directly because + /* Mortality: pptj is p.3 or p.death = trgradgp by cov by gradgp, variance of + * p.3=1-p..=1-sum i p.i overall mortality computed directly because * we compute the grad (wix pijx) instead of grad (pijx),even if - * wix is independent of theta. + * wix is independent of theta. */ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); for(j=nlstate+1;j<=nlstate+ndeath;j++) for(i=nlstate+1;i<=nlstate+ndeath;i++) - varppt[j][i]=doldmp[j][i]; + varppt[j][i]=doldmp[j][i]; /* This is the variance of p.3 */ /* end ppptj */ /* x centered again */ @@ -8843,15 +8901,15 @@ void concatwav(int wav[], int **dh, int hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres); 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]; + gmp[j] += prlim[i][i]*p3mat[i][j][1]; /* gmp[j] is p.3 */ } /* end probability of death */ fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ - fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j])); + fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));/* p.3 (STD p.3) */ for(i=1; i<=nlstate;i++){ - fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); + fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); /* wi, pi3 */ } } fprintf(ficresprobmorprev,"\n"); @@ -9910,7 +9968,10 @@ prevalence (with 95%% confidence interva fprintf(fichtm,"",subdirf2(optionfilefiname,"V_"), cpt,k1,nres); } fprintf(fichtm,"\n
    - Total life expectancy by age and \ -health expectancies in each live states (1 to %d). If popbased=1 the smooth (due to the model) \ +health expectancies in each live state (1 to %d) with confidence intervals \ +on left y-scale as well as proportions of time spent in each live state \ +(with confidence intervals) on right y-scale 0 to 100%%.\ + 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-%d.svg",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres); @@ -10264,18 +10325,18 @@ void printinggnuplot(char fileresu[], ch for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ fprintf(ficgp,"\nset label \"popbased %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",vpopbased,gplotlabel); if(vpopbased==0){ - fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); + fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nunset ytics; unset y2tics; set ytics nomirror; set y2tics 0,10,100;set y2range [0:100];\nplot [%.f:%.f] ",ageminpar,fage); }else fprintf(ficgp,"\nreplot "); - for (i=1; i<= nlstate+1 ; i ++) { + for (i=1; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */ k=2*i; - fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); - for (j=1; j<= nlstate+1 ; j ++) { - if (j==i) fprintf(ficgp," %%lf (%%lf)"); - else fprintf(ficgp," %%*lf (%%*lf)"); + fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */ + for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/ + if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */ + else fprintf(ficgp," %%*lf (%%*lf)"); /* skipping that field with a star */ } if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i); - else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); + else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); /* state=i-1=1 to nlstate */ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased); for (j=1; j<= nlstate+1 ; j ++) { if (j==i) fprintf(ficgp," %%lf (%%lf)"); @@ -10287,9 +10348,39 @@ void printinggnuplot(char fileresu[], ch if (j==i) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); } - if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0"); + if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); /* ,\\\n added for th percentage graphs */ else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); } /* state */ + /* again for the percentag spent in state i-1=1 to i-1=nlstate */ + for (i=2; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */ + k=2*i; + fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4)<=1 && ($4)>=0 ?($4)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */ + for (j=1; j<= nlstate ; j ++) + fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */ + for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/ + if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */ + else fprintf(ficgp," %%*lf (%%*lf)"); /* skipping that field with a star */ + } + if (i== 1) fprintf(ficgp,"\" t\"%%TLE\" w l lt %d axis x1y2, \\\n",i); /* Not used */ + else fprintf(ficgp,"\" t\"%%LE in state (%d)\" w l lw 2 lt %d axis x1y2, \\\n",i-1,i+1); /* state=i-1=1 to nlstate */ + fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4-$5*2)<=1 && ($4-$5*2)>=0? ($4-$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased); + for (j=1; j<= nlstate ; j ++) + fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */ + for (j=1; j<= nlstate+1 ; j ++) { + if (j==i) fprintf(ficgp," %%lf (%%lf)"); + else fprintf(ficgp," %%*lf (%%*lf)"); + } + fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,"); + fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4+$5*2)<=1 && ($4+$5*2)>=0 ? ($4+$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased); + for (j=1; j<= nlstate ; j ++) + fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */ + for (j=1; j<= nlstate+1 ; j ++) { + if (j==i) fprintf(ficgp," %%lf (%%lf)"); + else fprintf(ficgp," %%*lf (%%*lf)"); + } + if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2"); + else fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,\\\n"); + } /* state for percent */ } /* vpopbased */ fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; unset label;\n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */ } /* end nres */ @@ -12212,7 +12303,8 @@ double gompertz(double x[]) 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(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); + /* +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); */ /* To be seen */ } else printf("Gompertz cens[%d] neither 1 nor 0\n",i); /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ @@ -14569,6 +14661,7 @@ int main(int argc, char *argv[]) double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW; double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */ + double stdpercent; /* for computing the std error of percent e.i: e.i/e.. */ double fret; double dum=0.; /* Dummy variable */ /* double*** p3mat;*/ @@ -15765,9 +15858,27 @@ Interval (in months) between two waves: gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1) */ #endif #ifdef POWELL - powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); -#endif +#ifdef LINMINORIGINAL +#else /* LINMINORIGINAL */ + + flatdir=ivector(1,npar); + for (j=1;j<=npar;j++) flatdir[j]=0; +#endif /*LINMINORIGINAL */ + /* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */ + /* double h0=0.25; */ + macheps=pow(16.0,-13.0); + printf("Praxis Gegenfurtner mle=%d\n",mle); + fprintf(ficlog, "Praxis Gegenfurtner mle=%d\n", mle);fflush(ficlog); + /* ffmin = praxis(ftol,macheps, h0, npar, prin, p, gompertz); */ + /* For the Gompertz we use only two parameters */ + int _npar=2; + ffmin = praxis(ftol,macheps, h0, _npar, 4, p, gompertz); + printf("End Praxis\n"); fclose(ficrespow); +#ifdef LINMINORIGINAL +#else + free_ivector(flatdir,1,npar); +#endif /* LINMINORIGINAL*/ hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz); @@ -15902,7 +16013,7 @@ Please run with mle=-1 to get a correct fprintf(ficlog," + age*age "); fprintf(fichtm, "+ age*age"); } - for(j=1;j <=ncovmodel-2;j++){ + for(j=1;j <=ncovmodel-2-nagesqr;j++){ if(Typevar[j]==0) { printf(" + V%d ",Tvar[j]); fprintf(ficres," + V%d ",Tvar[j]); @@ -15973,7 +16084,7 @@ Please run with mle=-1 to get a correct fprintf(ficlog," + age*age "); fprintf(fichtm, "+ age*age"); } - for(j=1;j <=ncovmodel-2;j++){ + for(j=1;j <=ncovmodel-2-nagesqr;j++){ if(Typevar[j]==0) { printf(" + V%d ",Tvar[j]); fprintf(fichtm, "+ V%d",Tvar[j]); @@ -16740,16 +16851,21 @@ Please run with mle=-1 to get a correct 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 vpopbased=%d \n",vpopbased); - fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased); + printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased); + fprintf(ficlog, "varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased); + /* Call to varevsij to get cov(e.i, e.j)= vareij[i][j][(int)age]=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20 */ + /* Depending of popbased which changes the prevalences, either cross-sectional or period */ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* 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 "); + fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each state\n\ +# (these are 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); + 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 forward period (stable) prevalences in each health state \n"); + fprintf(ficrest,"the age specific forward period (stable) prevalences in each state) \n"); + fprintf(ficrest,"# with proportions of time spent in each state with standard error (on the right of the table.\n "); fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */ for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); + for (i=1;i<=nlstate;i++) fprintf(ficrest," %% e.%d/e.. (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"); */ printf("Computing age specific forward period (stable) prevalences in each health state \n"); @@ -16775,17 +16891,36 @@ Please run with mle=-1 to get a correct /*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]; + epj[nlstate+1] +=epj[j]; /* epp=sum_j epj = sum_j sum_i w_i e_ij */ } /* printf(" age %4.0f \n",age); */ - for(i=1, vepp=0.;i <=nlstate;i++) + for(i=1, vepp=0.;i <=nlstate;i++) /* Variance of total life expectancy e.. */ for(j=1;j <=nlstate;j++) - vepp += vareij[i][j][(int)age]; + vepp += vareij[i][j][(int)age]; /* sum_i sum_j cov(e.i, e.j) = var(e..) */ fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); + /* vareij[i][j] is the covariance cov(e.i, e.j) and vareij[j][j] is the variance of e.j */ for(j=1;j <=nlstate;j++){ fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); } + /* And proportion of time spent in state j */ + /* $$ E[r(X,Y)-E(r(X,Y))]^2=[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]' Var(X,Y)[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]$$ */ + /* \frac{\mu_x^2}{\mu_y^2} ( \frac{\sigma^2_x}{\mu_x^2}-2\frac{\sigma_{xy}}{\mu_x\mu_y} +\frac{\sigma^2_y}{\mu_y^2}) */ + /* \frac{e_{.i}^2}{e_{..}^2} ( \frac{\Var e_{.i}}{e_{.i}^2}-2\frac{\Var e_{.i} + \sum_{j\ne i} \Cov e_{.j},e_{.i}}{e_{.i}e_{..}} +\frac{\Var e_{..}}{e_{..}^2})*/ + /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp \sigmaxy= */ + /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstate+1]^4 */ + for(j=1;j <=nlstate;j++){ + /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */ + /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */ + + for(i=1,stdpercent=0.;i<=nlstate;i++){ /* Computing cov(e..,e.j)=cov(sum_i e.i,e.j)=sum_i cov(e.i, e.j) */ + stdpercent += vareij[i][j][(int)age]; + } + stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]* (vareij[j][j][(int)age]/epj[j]/epj[j]-2.*stdpercent/epj[j]/epj[nlstate+1]+ vepp/epj[nlstate+1]/epj[nlstate+1]); + /* stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]*(vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[nlstate+1]/epj[nlstate+1]); */ /* Without covariance */ + /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[nlstate+1]/epj[nlstate+1] + epj[j]*epj[j]*vepp/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1] )); */ + fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt(stdpercent)); + } fprintf(ficrest,"\n"); } } /* End vpopbased */