version 1.50, 2002/06/26 23:25:02
|
version 1.51, 2002/07/19 12:22:25
|
Line 2154 void varprob(char optionfilefiname[], do
|
Line 2154 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 2378 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 2394 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(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, mu2,mu1);
|
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, k2,l2,k1,l1);
|
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)) t \"%d\"",\
|
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),\
|
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
|
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(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, mu2,mu1);
|
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",\
|
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\"",\
|
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),\
|
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
|
mu2,std,v21,sqrt(lc1),v22,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));
|
|
}/* 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 */
|