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

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

version 1.360, 2024/04/30 10:59:22 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    Revision 1.360  2024/04/30 10:59:22  brouard
   Summary: Version 0.99s4 and estimation of std of e.j/e..    Summary: Version 0.99s4 and estimation of std of e.j/e..
   
Line 4414  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 4455  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 4479  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 8556  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 8573  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 8641  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 8710  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 8719  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 8748  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 8758  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 8800  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 8810  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 19 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 8846  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 14605  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 15801  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 16778  Please run with mle=-1 to get a correct Line 16824  Please run with mle=-1 to get a correct
         cptcod= 0; /* To be deleted */          cptcod= 0; /* To be deleted */
         printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);          printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
         fprintf(ficlog, "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 */          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\          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 ");  #  (these are weighted average of eij where weights are ");
Line 16814  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[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++){            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 */            /* 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}]$$ */            /* $$ 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 */            /* \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}) */
           /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp */            /* \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})*/
           /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstata+1]^4 */            /*\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++){            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] )); */
             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");            fprintf(ficrest,"\n");
         }          }
Removed from v.1.360  
changed lines
  Added in v.1.361

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