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

Diff for /imach/src/imach.c between versions 1.63 and 1.64

version 1.63, 2002/11/20 17:35:59 version 1.64, 2002/11/21 08:46:21
Line 944  double func( double *x) Line 944  double func( double *x)
           for (kk=1; kk<=cptcovage;kk++) {            for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];              cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
           }            }
           
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                         1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           savm=oldm;            savm=oldm;
           oldm=newm;            oldm=newm;
           
           
         } /* end mult */          } /* end mult */
               
         /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */          /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
Line 968  double func( double *x) Line 965  double func( double *x)
          */           */
         s1=s[mw[mi][i]][i];          s1=s[mw[mi][i]][i];
         s2=s[mw[mi+1][i]][i];          s2=s[mw[mi+1][i]][i];
         bbh=(double)bh[mi][i]/(double)stepm;          bbh=(double)bh[mi][i]/(double)stepm; 
           /* bias is positive if real duration
         /*      lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/           * is higher than the multiple of stepm and negative otherwise.
            */
           /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
           lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));  /* linear interpolation */
           /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
           /*if(lli ==000.0)*/
           /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
           ipmx +=1;
           sw += weight[i];
           ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
         } /* end of wave */
       } /* end of individual */
     }  else if(mle==2){
       for (i=1,ipmx=0, sw=0.; i<=imx; i++){
         for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
         for(mi=1; mi<= wav[i]-1; mi++){
           for (ii=1;ii<=nlstate+ndeath;ii++)
             for (j=1;j<=nlstate+ndeath;j++){
               oldm[ii][j]=(ii==j ? 1.0 : 0.0);
               savm[ii][j]=(ii==j ? 1.0 : 0.0);
             }
           for(d=0; d<=dh[mi][i]; d++){
             newm=savm;
             cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
             for (kk=1; kk<=cptcovage;kk++) {
               cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
             }
             out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                          1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
             savm=oldm;
             oldm=newm;
           } /* end mult */
         
           /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
           /* But now since version 0.9 we anticipate for bias and large stepm.
            * If stepm is larger than one month (smallest stepm) and if the exact delay 
            * (in months) between two waves is not a multiple of stepm, we rounded to 
            * the nearest (and in case of equal distance, to the lowest) interval but now
            * we keep into memory the bias bh[mi][i] and also the previous matrix product
            * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
            * probability in order to take into account the bias as a fraction of the way
            * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
            * -stepm/2 to stepm/2 .
            * For stepm=1 the results are the same as for previous versions of Imach.
            * For stepm > 1 the results are less biased than in previous versions. 
            */
           s1=s[mw[mi][i]][i];
           s2=s[mw[mi+1][i]][i];
           bbh=(double)bh[mi][i]/(double)stepm; 
           /* bias is positive if real duration
            * is higher than the multiple of stepm and negative otherwise.
            */
         lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */          lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
           /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
         /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));*/          /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-+bh)*out[s1][s2])); */ /* exponential interpolation */
           /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
           /*if(lli ==000.0)*/
           /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
           ipmx +=1;
           sw += weight[i];
           ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
         } /* end of wave */
       } /* end of individual */
     }  else if(mle==3){  /* exponential inter-extrapolation */
       for (i=1,ipmx=0, sw=0.; i<=imx; i++){
         for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
         for(mi=1; mi<= wav[i]-1; mi++){
           for (ii=1;ii<=nlstate+ndeath;ii++)
             for (j=1;j<=nlstate+ndeath;j++){
               oldm[ii][j]=(ii==j ? 1.0 : 0.0);
               savm[ii][j]=(ii==j ? 1.0 : 0.0);
             }
           for(d=0; d<dh[mi][i]; d++){
             newm=savm;
             cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
             for (kk=1; kk<=cptcovage;kk++) {
               cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
             }
             out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                          1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
             savm=oldm;
             oldm=newm;
           } /* end mult */
         
           /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
           /* But now since version 0.9 we anticipate for bias and large stepm.
            * If stepm is larger than one month (smallest stepm) and if the exact delay 
            * (in months) between two waves is not a multiple of stepm, we rounded to 
            * the nearest (and in case of equal distance, to the lowest) interval but now
            * we keep into memory the bias bh[mi][i] and also the previous matrix product
            * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
            * probability in order to take into account the bias as a fraction of the way
            * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
            * -stepm/2 to stepm/2 .
            * For stepm=1 the results are the same as for previous versions of Imach.
            * For stepm > 1 the results are less biased than in previous versions. 
            */
           s1=s[mw[mi][i]][i];
           s2=s[mw[mi+1][i]][i];
           bbh=(double)bh[mi][i]/(double)stepm; 
           /* bias is positive if real duration
            * is higher than the multiple of stepm and negative otherwise.
            */
           /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); */ /* linear interpolation */
           lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
         /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/          /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
         /*if(lli ==000.0)*/          /*if(lli ==000.0)*/
         /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */          /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
Line 983  double func( double *x) Line 1080  double func( double *x)
         ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       } /* end of wave */        } /* end of wave */
     } /* end of individual */      } /* end of individual */
   }  else{     }else{  /* ml=4 no inter-extrapolation */
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){      for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];        for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){        for(mi=1; mi<= wav[i]-1; mi++){
Line 1618  void  concatwav(int wav[], int **dh, int Line 1715  void  concatwav(int wav[], int **dh, int
         jk= j/stepm;          jk= j/stepm;
         jl= j -jk*stepm;          jl= j -jk*stepm;
         ju= j -(jk+1)*stepm;          ju= j -(jk+1)*stepm;
         if(jl <= -ju){          if(mle <=1){ 
           dh[mi][i]=jk;            if(jl==0){
           bh[mi][i]=jl;              dh[mi][i]=jk;
         }              bh[mi][i]=0;
         else{            }else{ /* We want a negative bias in order to only have interpolation ie
           dh[mi][i]=jk+1;                    * at the price of an extra matrix product in likelihood */
           bh[mi][i]=ju;              dh[mi][i]=jk+1;
         }              bh[mi][i]=ju;
         if(dh[mi][i]==0){            }
           dh[mi][i]=1; /* At least one step */          }else{
           bh[mi][i]=ju; /* At least one step */            if(jl <= -ju){
           printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);              dh[mi][i]=jk;
               bh[mi][i]=jl;       /* bias is positive if real duration
                                    * is higher than the multiple of stepm and negative otherwise.
                                    */
             }
             else{
               dh[mi][i]=jk+1;
               bh[mi][i]=ju;
             }
             if(dh[mi][i]==0){
               dh[mi][i]=1; /* At least one step */
               bh[mi][i]=ju; /* At least one step */
               printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);
             }
             if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);
         }          }
         if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);        } /* end if mle */
       }      } /* end wave */
     }  
   }    }
   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);
Removed from v.1.63  
changed lines
  Added in v.1.64

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