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

Diff for /imach/src/imach.c between versions 1.10 and 1.12

version 1.10, 2001/05/09 14:25:42 version 1.12, 2002/02/20 16:57:00
Line 8 Line 8
   Health expectancies are computed from the transistions observed between    Health expectancies are computed from the transistions observed between
   waves and are computed for each degree of severity of disability (number    waves and are computed for each degree of severity of disability (number
   of life states). More degrees you consider, more time is necessary to    of life states). More degrees you consider, more time is necessary to
   reach the Maximum Likelihood of the parameters involved in the model.    reach the Maximum Likelihood of the parameters involved in the model.
   The simplest model is the multinomial logistic model where pij is    The simplest model is the multinomial logistic model where pij is
   the probabibility to be observed in state j at the second wave conditional    the probabibility to be observed in state j at the second wave conditional
   to be observed in state i at the first wave. Therefore the model is:    to be observed in state i at the first wave. Therefore the model is:
Line 22 Line 22
   delay between waves is not identical for each individual, or if some    delay between waves is not identical for each individual, or if some
   individual missed an interview, the information is not rounded or lost, but    individual missed an interview, the information is not rounded or lost, but
   taken into account using an interpolation or extrapolation.    taken into account using an interpolation or extrapolation.
   hPijx is the probability to be    hPijx is the probability to be
   observed in state i at age x+h conditional to the observed state i at age    observed in state i at age x+h conditional to the observed state i at age
   x. The delay 'h' can be split into an exact number (nh*stepm) of    x. The delay 'h' can be split into an exact number (nh*stepm) of
   unobserved intermediate  states. This elementary transition (by month or    unobserved intermediate  states. This elementary transition (by month or
Line 693  double **prevalim(double **prlim, int nl Line 693  double **prevalim(double **prlim, int nl
   }    }
 }  }
   
 /*************** transition probabilities **********/  /*************** transition probabilities ***************/
   
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 {  {
Line 716  double **pmij(double **ps, double *cov, Line 716  double **pmij(double **ps, double *cov,
         s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];          s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
         /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/          /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
       }        }
       ps[i][j]=s2;        ps[i][j]=(s2);
     }      }
   }    }
       /*ps[3][2]=1;*/
   
   for(i=1; i<= nlstate; i++){    for(i=1; i<= nlstate; i++){
      s1=0;       s1=0;
     for(j=1; j<i; j++)      for(j=1; j<i; j++)
Line 740  double **pmij(double **ps, double *cov, Line 742  double **pmij(double **ps, double *cov,
     }      }
   }    }
   
   
   /*   for(ii=1; ii<= nlstate+ndeath; ii++){    /*   for(ii=1; ii<= nlstate+ndeath; ii++){
     for(jj=1; jj<= nlstate+ndeath; jj++){      for(jj=1; jj<= nlstate+ndeath; jj++){
      printf("%lf ",ps[ii][jj]);       printf("%lf ",ps[ii][jj]);
Line 804  double ***hpxij(double ***po, int nhstep Line 807  double ***hpxij(double ***po, int nhstep
       cov[1]=1.;        cov[1]=1.;
       cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;        cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
       for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];        for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
 for (k=1; k<=cptcovage;k++)        for (k=1; k<=cptcovage;k++)
         cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];          cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
    for (k=1; k<=cptcovprod;k++)        for (k=1; k<=cptcovprod;k++)
         cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];          cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
   
   
Line 913  void hesscov(double **matcov, double p[] Line 916  void hesscov(double **matcov, double p[]
   void lubksb(double **a, int npar, int *indx, double b[]) ;    void lubksb(double **a, int npar, int *indx, double b[]) ;
   void ludcmp(double **a, int npar, int *indx, double *d) ;    void ludcmp(double **a, int npar, int *indx, double *d) ;
   
   
   hess=matrix(1,npar,1,npar);    hess=matrix(1,npar,1,npar);
   
   printf("\nCalculation of the hessian matrix. Wait...\n");    printf("\nCalculation of the hessian matrix. Wait...\n");
Line 921  void hesscov(double **matcov, double p[] Line 923  void hesscov(double **matcov, double p[]
     printf("%d",i);fflush(stdout);      printf("%d",i);fflush(stdout);
     hess[i][i]=hessii(p,ftolhess,i,delti);      hess[i][i]=hessii(p,ftolhess,i,delti);
     /*printf(" %f ",p[i]);*/      /*printf(" %f ",p[i]);*/
       /*printf(" %lf ",hess[i][i]);*/
   }    }
    
   for (i=1;i<=npar;i++) {    for (i=1;i<=npar;i++) {
     for (j=1;j<=npar;j++)  {      for (j=1;j<=npar;j++)  {
       if (j>i) {        if (j>i) {
         printf(".%d%d",i,j);fflush(stdout);          printf(".%d%d",i,j);fflush(stdout);
         hess[i][j]=hessij(p,delti,i,j);          hess[i][j]=hessij(p,delti,i,j);
         hess[j][i]=hess[i][j];          hess[j][i]=hess[i][j];    
           /*printf(" %lf ",hess[i][j]);*/
       }        }
     }      }
   }    }
Line 1032  double hessii( double x[], double delta, Line 1036  double hessii( double x[], double delta,
     }      }
   }    }
   delti[theta]=delts;    delti[theta]=delts;
   return res;    return res;
     
 }  }
   
Line 1270  void  concatwav(int wav[], int **dh, int Line 1274  void  concatwav(int wav[], int **dh, int
   /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;    /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
      double sum=0., jmean=0.;*/       double sum=0., jmean=0.;*/
   
 int j, k=0,jk, ju, jl;    int j, k=0,jk, ju, jl;
      double sum=0.;    double sum=0.;
 jmin=1e+5;    jmin=1e+5;
  jmax=-1;    jmax=-1;
 jmean=0.;    jmean=0.;
   for(i=1; i<=imx; i++){    for(i=1; i<=imx; i++){
     mi=0;      mi=0;
     m=firstpass;      m=firstpass;
Line 1309  jmean=0.; Line 1313  jmean=0.;
           if(j==0) j=1;  /* Survives at least one month after exam */            if(j==0) j=1;  /* Survives at least one month after exam */
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           else if (j <= jmin)jmin=j;            if (j <= jmin) jmin=j;
           sum=sum+j;            sum=sum+j;
           if (j<0) printf("j=%d num=%d ",j,i);            /* if (j<10) printf("j=%d num=%d ",j,i); */
           }            }
         }          }
         else{          else{
Line 1319  jmean=0.; Line 1323  jmean=0.;
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           else if (j <= jmin)jmin=j;            else if (j <= jmin)jmin=j;
             /*   if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
           sum=sum+j;            sum=sum+j;
         }          }
         jk= j/stepm;          jk= j/stepm;
Line 1335  jmean=0.; Line 1340  jmean=0.;
   }    }
   jmean=sum/k;    jmean=sum/k;
   printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);    printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
 }   }
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)  void tricode(int *Tvar, int **nbcode, int imx)
 {  {
Line 1373  void tricode(int *Tvar, int **nbcode, in Line 1378  void tricode(int *Tvar, int **nbcode, in
   
  for (k=0; k<19; k++) Ndum[k]=0;   for (k=0; k<19; k++) Ndum[k]=0;
   
  for (i=1; i<=ncovmodel; i++) {   for (i=1; i<=ncovmodel-2; i++) {
       ij=Tvar[i];        ij=Tvar[i];
       Ndum[ij]++;        Ndum[ij]++;
     }      }
Line 1448  void varevsij(char fileres[], double *** Line 1453  void varevsij(char fileres[], double ***
   double **dnewm,**doldm;    double **dnewm,**doldm;
   int i, j, nhstepm, hstepm, h;    int i, j, nhstepm, hstepm, h;
   int k, cptcode;    int k, cptcode;
    double *xp;    double *xp;
   double **gp, **gm;    double **gp, **gm;
   double ***gradg, ***trgradg;    double ***gradg, ***trgradg;
   double ***p3mat;    double ***p3mat;
Line 1663  int main() Line 1668  int main()
   int ju,jl, mi;    int ju,jl, mi;
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;    int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;    int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
     
   int hstepm, nhstepm;    int hstepm, nhstepm;
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
Line 1772  split(pathtot, path,optionfile); Line 1777  split(pathtot, path,optionfile);
       fprintf(ficparo,"\n");        fprintf(ficparo,"\n");
     }      }
     
   npar= (nlstate+ndeath-1)*nlstate*ncovmodel;      npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
   
   p=param[1][1];    p=param[1][1];
     
   /* Reads comments: lines beginning with '#' */    /* Reads comments: lines beginning with '#' */
Line 1862  split(pathtot, path,optionfile); Line 1868  split(pathtot, path,optionfile);
     tab=ivector(1,NCOVMAX);      tab=ivector(1,NCOVMAX);
     ncodemax=ivector(1,8);      ncodemax=ivector(1,8);
   
     i=1;      i=1;
     while (fgets(line, MAXLINE, fic) != NULL)    {      while (fgets(line, MAXLINE, fic) != NULL)    {
       if ((i >= firstobs) && (i <=lastobs)) {        if ((i >= firstobs) && (i <=lastobs)) {
                 
Line 1884  split(pathtot, path,optionfile); Line 1890  split(pathtot, path,optionfile);
           cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);            cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
         }          }
         num[i]=atol(stra);          num[i]=atol(stra);
          
         /*printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/          /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
             printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
   
         i=i+1;          i=i+1;
       }        }
     }      }
       /* printf("ii=%d", ij);
     /*scanf("%d",i);*/         scanf("%d",i);*/
   imx=i-1; /* Number of individuals */    imx=i-1; /* Number of individuals */
   
     /* for (i=1; i<=imx; i++){
       if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
       if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
       if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
     }
     for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/
   
   /* Calculation of the number of parameter from char model*/    /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,15);    Tvar=ivector(1,15);
   Tprod=ivector(1,15);    Tprod=ivector(1,15);
Line 1975  split(pathtot, path,optionfile); Line 1989  split(pathtot, path,optionfile);
     }      }
     /*-calculation of age at interview from date of interview and age at death -*/      /*-calculation of age at interview from date of interview and age at death -*/
     agev=matrix(1,maxwav,1,imx);      agev=matrix(1,maxwav,1,imx);
   
      for (i=1; i<=imx; i++)
        for(m=2; (m<= maxwav); m++)
          if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
            anint[m][i]=9999;
            s[m][i]=-1;
          }
         
     for (i=1; i<=imx; i++)  {      for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);        agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
Line 2083  printf("Total number of individuals= %d, Line 2104  printf("Total number of individuals= %d,
     newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */      oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
           
     /* For Powell, parameters are in a vector p[] starting at p[1]      /* For Powell, parameters are in a vector p[] starting at p[1]
        so we point p on param[1][1] so that p[1] maps on param[1][1][1] */         so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
     p=param[1][1]; /* *(*(*(param +1)+1)+0) */      p=param[1][1]; /* *(*(*(param +1)+1)+0) */
Line 2666  strcpy(fileresvpl,"vpl"); Line 2687  strcpy(fileresvpl,"vpl");
   /*printf("Total time was %d uSec.\n", total_usecs);*/    /*printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/    /*------ End -----------*/
   
   
  end:   end:
 #ifdef windows  #ifdef windows
  chdir(pathcd);   chdir(pathcd);
 #endif  #endif
  /*system("wgnuplot graph.plt");*/   /*system("wgnuplot graph.plt");*/
  system("../gp37mgw/wgnuplot graph.plt");   /*system("../gp37mgw/wgnuplot graph.plt");*/
    /*system("cd ../gp37mgw");*/
    system("..\\gp37mgw\\wgnuplot graph.plt");
   
 #ifdef windows  #ifdef windows
   while (z[0] != 'q') {    while (z[0] != 'q') {
Removed from v.1.10  
changed lines
  Added in v.1.12

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