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

Diff for /imach/src/imach.c between versions 1.50 and 1.52

version 1.50, 2002/06/26 23:25:02 version 1.52, 2002/07/19 18:49:30
Line 1862  void varevsij(char optionfilefiname[], d Line 1862  void varevsij(char optionfilefiname[], d
     exit(0);      exit(0);
   }    }
   else{    else{
     fprintf(fichtm,"\n<li><h4> Computing step probabilities of dying and weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");      fprintf(fichtm,"\n<li><h4> Computing probabilities of dying as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
   }    }
   varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
   
Line 2013  void varevsij(char optionfilefiname[], d Line 2013  void varevsij(char optionfilefiname[], d
     /* This for computing force of mortality (h=1)as a weighted average */      /* This for computing force of mortality (h=1)as a weighted average */
     for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){      for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){
       for(i=1; i<= nlstate; i++)        for(i=1; i<= nlstate; i++)
         gmp[j] += prlim[i][i]*p3mat[i][j][1];          gmp[j] += prlim[i][i]*p3mat[i][j][1];
     }          }    
     /* end force of mortality */      /* end force of mortality */
   
Line 2049  void varevsij(char optionfilefiname[], d Line 2049  void varevsij(char optionfilefiname[], d
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm);    fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm);
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm);    fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm);
   fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev);    fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev);
   fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);    fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,digitp,digit);
     /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);
   */
   fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit);    fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit);
   
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
Line 2154  void varprob(char optionfilefiname[], do Line 2156  void varprob(char optionfilefiname[], do
   int k2, l2, j1,  z1;    int k2, l2, j1,  z1;
   int k=0,l, cptcode;    int k=0,l, cptcode;
   int first=1, first1;    int first=1, first1;
   double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2;    double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
   double **dnewm,**doldm;    double **dnewm,**doldm;
   double *xp;    double *xp;
   double *gp, *gm;    double *gp, *gm;
Line 2378  void varprob(char optionfilefiname[], do Line 2380  void varprob(char optionfilefiname[], do
   
       /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/        /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
       first1=1;        first1=1;
       for (k1=1; k1<=(nlstate);k1++){        for (k2=1; k2<=(nlstate);k2++){
         for (l1=1; l1<=(nlstate+ndeath);l1++){          for (l2=1; l2<=(nlstate+ndeath);l2++){
           if(l1==k1) continue;            if(l2==k2) continue;
           i=(k1-1)*(nlstate+ndeath)+l1;            j=(k2-1)*(nlstate+ndeath)+l2;
           for (k2=1; k2<=(nlstate);k2++){            for (k1=1; k1<=(nlstate);k1++){
             for (l2=1; l2<=(nlstate+ndeath);l2++){              for (l1=1; l1<=(nlstate+ndeath);l1++){
               if(l2==k2) continue;                if(l1==k1) continue;
               j=(k2-1)*(nlstate+ndeath)+l2;                i=(k1-1)*(nlstate+ndeath)+l1;
               if(j<=i) continue;                if(i<=j) continue;
               for (age=bage; age<=fage; age ++){                for (age=bage; age<=fage; age ++){
                 if ((int)age %5==0){                  if ((int)age %5==0){
                   v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM;                    v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM;
Line 2394  void varprob(char optionfilefiname[], do Line 2396  void varprob(char optionfilefiname[], do
                   cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM;                    cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM;
                   mu1=mu[i][(int) age]/stepm*YEARM ;                    mu1=mu[i][(int) age]/stepm*YEARM ;
                   mu2=mu[j][(int) age]/stepm*YEARM;                    mu2=mu[j][(int) age]/stepm*YEARM;
                     c12=cv12/sqrt(v1*v2);
                   /* Computing eigen value of matrix of covariance */                    /* Computing eigen value of matrix of covariance */
                   lc1=(v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));                    lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
                   lc2=(v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));                    lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
                   if(first1==1){  
                     first1=0;  
                     printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\nOthers in log...\n",v1,v2,cv12,lc1,lc2);  
                   }  
                   fprintf(ficlog,"Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);  
                   /* Eigen vectors */                    /* Eigen vectors */
                   v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));                    v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
                   v21=sqrt(1.-v11*v11);                    /*v21=sqrt(1.-v11*v11); *//* error */
                     v21=(lc1-v1)/cv12*v11;
                   v12=-v21;                    v12=-v21;
                   v22=v11;                    v22=v11;
                     tnalp=v21/v11;
                     if(first1==1){
                       first1=0;
                       printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);
                     }
                     fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);
                   /*printf(fignu*/                    /*printf(fignu*/
                   /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */                    /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */
                   /* mu2+ v21*lc1*cost + v21*lc2*sin(t) */                    /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */
                   if(first==1){                    if(first==1){
                     first=0;                      first=0;
                     fprintf(ficgp,"\nset parametric;set nolabel");                      fprintf(ficgp,"\nset parametric;unset label");
                     fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k2,l2,k1,l1);                      fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
                     fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");                      fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
                     fprintf(fichtm,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup> :<a href=\"varpijgr%s%d%1d%1d-%1d%1d.png\">varpijgr%s%d%1d%1d-%1d%1d.png</A>, ",k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1);                      fprintf(fichtm,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup> :<a href=\"varpijgr%s%d%1d%1d-%1d%1d.png\">varpijgr%s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,optionfilefiname, j1,k1,l1,k2,l2,optionfilefiname, j1,k1,l1,k2,l2);
                     fprintf(fichtm,"\n<br><img src=\"varpijgr%s%d%1d%1d-%1d%1d.png\"> ",optionfilefiname, j1,k2,l2,k1,l1);                      fprintf(fichtm,"\n<br><img src=\"varpijgr%s%d%1d%1d-%1d%1d.png\"> ",optionfilefiname, j1,k1,l1,k2,l2);
                     fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k2,l2,k1,l1);                      fprintf(fichtm,"\n<br> Correlation at age %d (%.3f),",(int) age, c12);
                     fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);                      fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k1,l1,k2,l2);
                     fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);                      fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
                     /*              fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\                      fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
                             mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \                      fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\
                             mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);                              mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),\
                     */                              mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
                     fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\  
                             mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \  
                             mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));  
                   }else{                    }else{
                     first=0;                      first=0;
                     fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);                      fprintf(fichtm," %d (%.3f),",(int) age, c12);
                     fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);                      fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
                     /*                      fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
                     fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\                      fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\
                             mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \                              mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),\
                             mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);                              mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
                     */  
                     fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\  
                             mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \  
                             mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));  
                   }/* if first */                    }/* if first */
                 } /* age mod 5 */                  } /* age mod 5 */
               } /* end loop age */                } /* end loop age */
               fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k2,l2,k1,l1);                fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k1,l1,k2,l2);
               first=1;                first=1;
             } /*l12 */              } /*l12 */
           } /* k12 */            } /* k12 */
Removed from v.1.50  
changed lines
  Added in v.1.52

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