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

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

version 1.60, 2002/11/19 00:17:38 version 1.64, 2002/11/21 08:46:21
Line 928  double func( double *x) Line 928  double func( double *x)
   cov[1]=1.;    cov[1]=1.;
   
   for(k=1; k<=nlstate; k++) ll[k]=0.;    for(k=1; k<=nlstate; k++) ll[k]=0.;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){  
     for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];    if(mle==1){
     for(mi=1; mi<= wav[i]-1; mi++){      for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (ii=1;ii<=nlstate+ndeath;ii++)        for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
         for (j=1;j<=nlstate+ndeath;j++){        for(mi=1; mi<= wav[i]-1; mi++){
           oldm[ii][j]=(ii==j ? 1.0 : 0.0);          for (ii=1;ii<=nlstate+ndeath;ii++)
           savm[ii][j]=(ii==j ? 1.0 : 0.0);            for (j=1;j<=nlstate+ndeath;j++){
         }              oldm[ii][j]=(ii==j ? 1.0 : 0.0);
       for(d=0; d<dh[mi][i]; d++){              savm[ii][j]=(ii==j ? 1.0 : 0.0);
         newm=savm;            }
         cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;          for(d=0; d<dh[mi][i]; d++){
         for (kk=1; kk<=cptcovage;kk++) {            newm=savm;
           cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];            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));            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
         savm=oldm;                         1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
         oldm=newm;            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]>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]>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]);*/
           /*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{  /* ml=4 no 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];
             }
                   
       } /* end mult */            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 */          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.          ipmx +=1;
        * If stepm is larger than one month (smallest stepm) and if the exact delay           sw += weight[i];
        * (in months) between two waves is not a multiple of stepm, we rounded to           ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
        * the nearest (and in case of equal distance, to the lowest) interval but now        } /* end of wave */
        * we keep into memory the bias bh[mi][i] and also the previous matrix product      } /* end of individual */
        * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the    } /* End of if */
        * 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;  
       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]));  
       /*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=(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 */  
   
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];    for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */    /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
   l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */    l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
Line 1584  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);
Line 3122  populforecast(char fileres[], double anp Line 3266  populforecast(char fileres[], double anp
   
 int main(int argc, char *argv[])  int main(int argc, char *argv[])
 {  {
     int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
   int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;    int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
   double agedeb, agefin,hf;    double agedeb, agefin,hf;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;    double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
Line 3160  int main(int argc, char *argv[]) Line 3304  int main(int argc, char *argv[])
   double *epj, vepp;    double *epj, vepp;
   double kk1, kk2;    double kk1, kk2;
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;    double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
   /*int *movingaverage; */  
   
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
Line 3663  int main(int argc, char *argv[]) Line 3806  int main(int argc, char *argv[])
      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) */
   
   if(mle==1){    if(mle>=1){ /* Could be 1 or 2 */
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);      mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
   }    }
           
Removed from v.1.60  
changed lines
  Added in v.1.64

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