[BACK]Return to imach.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.359 and 1.361

version 1.359, 2024/04/24 21:21:17 version 1.361, 2024/05/12 20:29:32
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     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    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    Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes
   
Line 1394  double gnuplotversion=GNUPLOTVERSION; Line 1409  double gnuplotversion=GNUPLOTVERSION;
 /* $State$ */  /* $State$ */
 #include "version.h"  #include "version.h"
 char version[]=__IMACH_VERSION__;  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 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$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 4411  void powell(double p[], double **xi, int Line 4426  void powell(double p[], double **xi, int
 #endif  #endif
 #ifdef POWELLORIGINAL  #ifdef POWELLORIGINAL
       if (t < 0.0) { /* Then we use it for new direction */        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 */        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("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);          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,"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);          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  #endif
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
         printf("Before linmin in direction P%d-P0\n",n);          printf("Before linmin in direction P%d-P0\n",n);
Line 4452  void powell(double p[], double **xi, int Line 4467  void powell(double p[], double **xi, int
           xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */            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 */            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  #ifdef LINMINORIGINAL
 #else  #else
         for (j=1, flatd=0;j<=n;j++) {          for (j=1, flatd=0;j<=n;j++) {
Line 4476  void powell(double p[], double **xi, int Line 4506  void powell(double p[], double **xi, int
           free_vector(pt,1,n);             free_vector(pt,1,n); 
           return;            return;
 #endif  #endif
         }          }  /* endif(flatd >0) */
 #endif  #endif /* LINMINORIGINAL */
         printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          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);          fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
                   
Line 8553  void  concatwav(int wav[], int **dh, int Line 8583  void  concatwav(int wav[], int **dh, int
 /************ Variance ******************/  /************ 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)   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 **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);
     * double **newm;      * double **newm;
     * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)       * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) 
Line 8570  void  concatwav(int wav[], int **dh, int Line 8602  void  concatwav(int wav[], int **dh, int
    double ***gradg, ***trgradg; /**< for var eij */     double ***gradg, ***trgradg; /**< for var eij */
    double **gradgp, **trgradgp; /**< for var p point j */     double **gradgp, **trgradgp; /**< for var p point j */
    double *gpp, *gmp; /**< 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 ***p3mat;
    double age,agelim, hf;     double age,agelim, hf;
    /* double ***mobaverage; */     /* double ***mobaverage; */
Line 8638  void  concatwav(int wav[], int **dh, int Line 8670  void  concatwav(int wav[], int **dh, int
    fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");     fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
    fprintf(fichtm,"\n<br>%s  <br>\n",digitp);     fprintf(fichtm,"\n<br>%s  <br>\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);     pstamp(ficresvij);
    fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");     fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
    if(popbased==1)     if(popbased==1)
Line 8707  void  concatwav(int wav[], int **dh, int Line 8739  void  concatwav(int wav[], int **dh, int
              prlim[i][i]=mobaverage[(int)age][i][ij];               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 */         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         /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability
Line 8716  void  concatwav(int wav[], int **dh, int Line 8748  void  concatwav(int wav[], int **dh, int
        for(j=1; j<= nlstate; j++){         for(j=1; j<= nlstate; j++){
          for(h=0; h<=nhstepm; h++){           for(h=0; h<=nhstepm; h++){
            for(i=1, gp[h][j]=0.;i<=nlstate;i++)             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         /* Next for computing shifted+ probability of death (h=1 means
           computed over hstepm matrices product = hstepm*stepm months)             computed over hstepm matrices product = hstepm*stepm months) 
           as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 .            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++)           for(i=1,gpp[j]=0.; i<= nlstate; i++)
            gpp[j] += prlim[i][i]*p3mat[i][j][1];             gpp[j] += prlim[i][i]*p3mat[i][j][1];
        }         }
Line 8745  void  concatwav(int wav[], int **dh, int Line 8777  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(h=0; h<=nhstepm; h++){
            for(i=1, gm[h][j]=0.;i<=nlstate;i++)             for(i=1, gm[h][j]=0.;i<=nlstate;i++)
              gm[h][j] += prlim[i][i]*p3mat[i][j][h];               gm[h][j] += prlim[i][i]*p3mat[i][j][h];
Line 8755  void  concatwav(int wav[], int **dh, int Line 8787  void  concatwav(int wav[], int **dh, int
        }         }
        /* This for computing probability of death (h=1 means         /* This for computing probability of death (h=1 means
           computed over hstepm matrices product = hstepm*stepm months)             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++)           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];
        }             }    
        /* end shifting computations */         /* 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++){           for(h=0; h<=nhstepm; h++){
            gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];             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];           gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
        }         }
                                                   
      } /* End theta */       } /* End theta */
             
      /* We got the gradient matrix for each theta and state j */                       /* 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); /* veij */       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(j=1; j<=nlstate;j++)
          for(theta=1; theta <=npar; theta++)           for(theta=1; theta <=npar; theta++)
            trgradg[h][j][theta]=gradg[h][theta][j];             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++)         for(theta=1; theta <=npar; theta++)
          trgradgp[j][theta]=gradgp[theta][j];           trgradgp[j][theta]=gradgp[theta][j];
      /**< as well as its transposed matrix        /**< as well as its transposed matrix 
Line 8797  void  concatwav(int wav[], int **dh, int Line 8831  void  concatwav(int wav[], int **dh, int
          vareij[i][j][(int)age] =0.;           vareij[i][j][(int)age] =0.;
   
      /* Computing trgradg by matcov by gradg at age and summing over h       /* Computing trgradg by matcov by gradg at age and summing over h
       * and k (nhstepm) formula 15 of article        * and k (nhstepm) formula 32 of article
       * Lievre-Brouard-Heathcote        * 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++){       for(h=0;h<=nhstepm;h++){
Line 8807  void  concatwav(int wav[], int **dh, int Line 8844  void  concatwav(int wav[], int **dh, int
          matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);           matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
          for(i=1;i<=nlstate;i++)           for(i=1;i<=nlstate;i++)
            for(j=1;j<=nlstate;j++)             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       /* Mortality: pptj is p.3 or p.death = trgradgp by cov by gradgp, variance of
       * p.j overall mortality formula 49 but computed directly because        * 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        * 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(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);       matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
      for(j=nlstate+1;j<=nlstate+ndeath;j++)       for(j=nlstate+1;j<=nlstate+ndeath;j++)
        for(i=nlstate+1;i<=nlstate+ndeath;i++)         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 */       /* end ppptj */
      /*  x centered again */       /*  x centered again */
                                   
Line 8843  void  concatwav(int wav[], int **dh, int Line 8881  void  concatwav(int wav[], int **dh, int
      hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres);         hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres);  
      for(j=nlstate+1;j<=nlstate+ndeath;j++){       for(j=nlstate+1;j<=nlstate+ndeath;j++){
        for(i=1,gmp[j]=0.;i<= nlstate; i++)          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 */       /* end probability of death */
                                   
      fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);       fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
      for(j=nlstate+1; j<=(nlstate+ndeath);j++){       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++){         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");       fprintf(ficresprobmorprev,"\n");
Line 9910  prevalence (with 95%% confidence interva Line 9948  prevalence (with 95%% confidence interva
        fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \       fprintf(fichtm,"\n<br>- 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\  true period expectancies (those weighted with period prevalences are also\
  drawn in addition to the population based expectancies computed using\   drawn in addition to the population based expectancies computed using\
  observed and cahotic prevalences:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);   observed and cahotic prevalences:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
Line 10264  void printinggnuplot(char fileresu[], ch Line 10305  void printinggnuplot(char fileresu[], ch
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/        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);          fprintf(ficgp,"\nset label \"popbased %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",vpopbased,gplotlabel);
         if(vpopbased==0){          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          }else
           fprintf(ficgp,"\nreplot ");            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;            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);            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 ++) {            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)");              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)");              else fprintf(ficgp," %%*lf (%%*lf)");  /* skipping that field with a star */
           }               }   
           if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);            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);            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 ++) {            for (j=1; j<= nlstate+1 ; j ++) {
             if (j==i) fprintf(ficgp," %%lf (%%lf)");              if (j==i) fprintf(ficgp," %%lf (%%lf)");
Line 10287  void printinggnuplot(char fileresu[], ch Line 10328  void printinggnuplot(char fileresu[], ch
             if (j==i) fprintf(ficgp," %%lf (%%lf)");              if (j==i) fprintf(ficgp," %%lf (%%lf)");
             else 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");            else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
         } /* state */          } /* 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 */        } /* vpopbased */
       fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; unset label;\n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */        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 */      } /* end nres */
Line 14569  int main(int argc, char *argv[]) Line 14640  int main(int argc, char *argv[])
   double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;    double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
   double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */    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 fret;
   double dum=0.; /* Dummy variable */    double dum=0.; /* Dummy variable */
   /* double*** p3mat;*/    /* double*** p3mat;*/
Line 15765  Interval (in months) between two waves: Line 15837  Interval (in months) between two waves:
     gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1)  */      gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1)  */
 #endif  #endif
 #ifdef POWELL  #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);       powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
 #endif    #endif  
     fclose(ficrespow);      fclose(ficrespow);
   #ifdef LINMINORIGINAL
   #else
         free_ivector(flatdir,1,npar); 
   #endif  /* LINMINORIGINAL*/
           
     hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz);       hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz); 
   
Line 16740  Please run with mle=-1 to get a correct Line 16822  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*/        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
         oldm=oldms;savm=savms; /* ZZ Segmentation fault */          oldm=oldms;savm=savms; /* ZZ Segmentation fault */
         cptcod= 0; /* To be deleted */          cptcod= 0; /* To be deleted */
         printf("varevsij vpopbased=%d \n",vpopbased);          printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
         fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased);          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 */          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)          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          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? */          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 (std) ",i);
           for (i=1;i<=nlstate;i++) fprintf(ficrest," %% e.%d/e.. (std) ",i);
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
         /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\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");          printf("Computing age specific forward period (stable) prevalences in each health state \n");
Line 16775  Please run with mle=-1 to get a correct Line 16862  Please run with mle=-1 to get a correct
               /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/                /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
               /* 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); */            /* 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++)              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));            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++){            for(j=1;j <=nlstate;j++){
             fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));              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");            fprintf(ficrest,"\n");
         }          }
       } /* End vpopbased */        } /* End vpopbased */
Removed from v.1.359  
changed lines
  Added in v.1.361

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>