--- imach/src/imach.c 2002/06/26 23:25:02 1.50 +++ imach/src/imach.c 2002/07/19 12:22:25 1.51 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.50 2002/06/26 23:25:02 brouard Exp $ +/* $Id: imach.c,v 1.51 2002/07/19 12:22:25 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -2154,7 +2154,7 @@ void varprob(char optionfilefiname[], do int k2, l2, j1, z1; int k=0,l, cptcode; 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 *xp; double *gp, *gm; @@ -2378,15 +2378,15 @@ void varprob(char optionfilefiname[], do /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ first1=1; - for (k1=1; k1<=(nlstate);k1++){ - for (l1=1; l1<=(nlstate+ndeath);l1++){ - if(l1==k1) continue; - i=(k1-1)*(nlstate+ndeath)+l1; - for (k2=1; k2<=(nlstate);k2++){ - for (l2=1; l2<=(nlstate+ndeath);l2++){ - if(l2==k2) continue; - j=(k2-1)*(nlstate+ndeath)+l2; - if(j<=i) continue; + for (k2=1; k2<=(nlstate);k2++){ + for (l2=1; l2<=(nlstate+ndeath);l2++){ + if(l2==k2) continue; + j=(k2-1)*(nlstate+ndeath)+l2; + for (k1=1; k1<=(nlstate);k1++){ + for (l1=1; l1<=(nlstate+ndeath);l1++){ + if(l1==k1) continue; + i=(k1-1)*(nlstate+ndeath)+l1; + if(i<=j) continue; for (age=bage; age<=fage; age ++){ if ((int)age %5==0){ v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM; @@ -2394,55 +2394,49 @@ void varprob(char optionfilefiname[], do cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM; mu1=mu[i][(int) age]/stepm*YEARM ; mu2=mu[j][(int) age]/stepm*YEARM; + c12=cv12/sqrt(v1*v2); /* Computing eigen value of matrix of covariance */ - lc1=(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)); - 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); + 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)))/2.; /* Eigen vectors */ 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; 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*/ /* 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){ first=0; - 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 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)\"",k1,l1,k2,l2); 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%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\"",\ - mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \ - mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age); - */ - 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)); + 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, ",k1,l1,k2,l2,optionfilefiname, j1,k1,l1,k2,l2,optionfilefiname, j1,k1,l1,k2,l2); + fprintf(fichtm,"\n
",optionfilefiname, j1,k1,l1,k2,l2); + fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k1,l1,k2,l2); + fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); + fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); + 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),\ + mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); }else{ first=0; - 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, mu2,mu1); - /* - 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\"",\ - mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \ - mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age); - */ - 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)); + 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)) not",\ + mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),\ + mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); }/* if first */ } /* age mod 5 */ } /* 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; } /*l12 */ } /* k12 */