--- imach/src/imach.c 2024/04/30 10:59:22 1.360
+++ imach/src/imach.c 2024/05/12 20:29:32 1.361
@@ -1,6 +1,18 @@
-/* $Id: imach.c,v 1.360 2024/04/30 10:59:22 brouard Exp $
+/* $Id: imach.c,v 1.361 2024/05/12 20:29:32 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ 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..
@@ -1393,12 +1405,12 @@ double gnuplotversion=GNUPLOTVERSION;
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.360 2024/04/30 10:59:22 brouard Exp $ */
+/* $Id: imach.c,v 1.361 2024/05/12 20:29:32 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
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.360 $ $Date: 2024/04/30 10:59:22 $";
+char fullversion[]="$Revision: 1.361 $ $Date: 2024/05/12 20:29:32 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -4414,14 +4426,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);
@@ -4455,6 +4467,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++) {
@@ -4479,8 +4506,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);
@@ -8556,7 +8583,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)
@@ -8573,7 +8602,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 e.. nlstate+1 to nlstate+ndeath */
double ***p3mat;
double age,agelim, hf;
/* double ***mobaverage; */
@@ -8641,7 +8670,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)
@@ -8710,7 +8739,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
@@ -8719,14 +8748,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];
}
@@ -8748,9 +8777,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];
@@ -8758,37 +8787,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
@@ -8800,8 +8831,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++){
@@ -8810,20 +8844,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 19 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 */
@@ -8846,15 +8881,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");
@@ -14605,6 +14640,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;*/
@@ -15801,9 +15837,19 @@ Interval (in months) between two waves:
gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1) */
#endif
#ifdef POWELL
+#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);
#endif
fclose(ficrespow);
+#ifdef LINMINORIGINAL
+#else
+ free_ivector(flatdir,1,npar);
+#endif /* LINMINORIGINAL*/
hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz);
@@ -16778,6 +16824,8 @@ Please run with mle=-1 to get a correct
cptcod= 0; /* To be deleted */
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 state\n\
# (these are weighted average of eij where weights are ");
@@ -16814,26 +16862,35 @@ 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[j][i] is the variance of epj */
+ /* 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}]$$ */
- /* \sigma^2_x/\mu_y^2 +\sigma^2_y \mu^2x/\mu_y^4 */
- /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp */
- /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstata+1]^4 */
+ /* \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[nlstate+1]/epj[nlstate+1] + 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( 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");
}