--- imach/src/imach.c 2002/06/10 13:12:49 1.48 +++ imach/src/imach.c 2002/06/20 14:03:39 1.49 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.48 2002/06/10 13:12:49 brouard Exp $ +/* $Id: imach.c,v 1.49 2002/06/20 14:03:39 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -84,7 +84,7 @@ char version[80]="Imach version 0.8h, May 2002, INED-EUROREVES "; int erreur; /* Error number */ int nvar; -int cptcovn, cptcovage=0, cptcoveff=0,cptcov; +int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ @@ -1939,7 +1939,7 @@ void varprevlim(char fileres[], double * /************ Variance of one-step probabilities ******************/ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax) { - int i, j, i1, k1, l1; + int i, j=0, i1, k1, l1, t, tj; int k2, l2, j1, z1; int k=0,l, cptcode; int first=1; @@ -2012,37 +2012,43 @@ void varprob(char optionfilefiname[], do exit(0); } else{ - fprintf(fichtm,"\n

Computing matrix of variance-covariance of step probabilities

\n"); - fprintf(fichtm,"\n
We have drawn ellipsoids of confidence around the pij, pkl to understand the covariance between two incidences. They are expressed in year-1 in order to be less dependent of stepm.
\n"); + fprintf(fichtm,"\n
  • Computing matrix of variance-covariance of step probabilities

  • \n"); + fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the pij, pkl to understand the covariance between two incidences. They are expressed in year-1 in order to be less dependent of stepm.
    \n"); fprintf(fichtm,"\n
    We have drawn x'cov-1x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis.
    When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.
    \n"); } + + cov[1]=1; - j=cptcoveff; - if (cptcovn<1) {j=1;ncodemax[1]=1;} + tj=cptcoveff; + if (cptcovn<1) {tj=1;ncodemax[1]=1;} j1=0; - for(k1=1; k1<=1;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ - j1++; - - if (cptcovn>0) { - fprintf(ficresprob, "\n#********** Variable "); - fprintf(ficresprobcov, "\n#********** Variable "); - fprintf(ficgp, "\n#********** Variable "); - fprintf(fichtm, "\n

    ********** Variable

    \n "); - fprintf(ficresprobcor, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficresprob, "**********\n#"); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficresprobcov, "**********\n#"); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficgp, "**********\n#"); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficgp, "**********\n#"); - for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(fichtm, "**********\n#"); - } - + for(t=1; t<=tj;t++){ + for(i1=1; i1<=ncodemax[t];i1++){ + j1++; + + if (cptcovn>0) { + fprintf(ficresprob, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprob, "**********\n#"); + fprintf(ficresprobcov, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprobcov, "**********\n#"); + + fprintf(ficgp, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficgp, "**********\n#"); + + + fprintf(fichtm, "\n
    ********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(fichtm, "**********\n
    "); + + fprintf(ficresprobcor, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficgp, "**********\n#"); + } + for (age=bage; age<=fage; age ++){ cov[2]=age; for (k=1; k<=cptcovn;k++) { @@ -2169,9 +2175,9 @@ void varprob(char optionfilefiname[], do fprintf(ficgp,"\nset parametric;set nolabel"); 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 ter png small\nset size 0.65,0.65"); - fprintf(fichtm,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1 :varpijgr%s%1d%1d-%1d%1d.png, ",k2,l2,k1,l1,optionfilefiname,k2,l2,k1,l1,optionfilefiname,k2,l2,k1,l1); - fprintf(fichtm,"\n
    , ",optionfilefiname,k2,l2,k1,l1); - fprintf(ficgp,"\nset out \"varpijgr%s%1d%1d-%1d%1d.png\"",optionfilefiname,k2,l2,k1,l1); + fprintf(fichtm,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1 :varpijgr%s%d%1d%1d-%1d%1d.png, ",k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1); + fprintf(fichtm,"\n
    ",optionfilefiname, j1,k2,l2,k1,l1); + fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k2,l2,k1,l1); 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, k2,l2,k1,l1); 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\"",\ @@ -2187,7 +2193,7 @@ void varprob(char optionfilefiname[], do }/* if first */ } /* age mod 5 */ } /* end loop age */ - fprintf(ficgp,"\nset out \"varpijgr%s%1d%1d-%1d%1d.png\";replot;",optionfilefiname,k2,l2,k1,l1); + fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k2,l2,k1,l1); first=1; } /*l12 */ } /* k12 */ @@ -2223,7 +2229,7 @@ void printinghtml(char fileres[], char t printf("Problem with %s \n",optionfilehtm), exit(0); } - fprintf(fichtm,"