mle=1 corrects bias

mle=2 no correction

diff --git a/src/imach.c b/src/imach.c
index 0855c36798b5701b0d2c772302f715019052f4a6..a18c3c772c6f3d3eca2801716532351384d5503a 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -928,56 +928,86 @@ double func( double *x)
   cov[1]=1.;
 
   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];
-    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];
-       }
+
+  if(mle==1){
+    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;
+         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 */
+       } /* 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;
-      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 */
-
+       /*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;
+       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 */
+  }  else{ 
+    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 */
+       ipmx +=1;
+       sw += weight[i];
+       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+      } /* end of wave */
+    } /* end of individual */
+  } /* End of if */
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* 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 */
@@ -3122,7 +3152,7 @@ populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double
 
 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;
   double agedeb, agefin,hf;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
@@ -3160,7 +3190,6 @@ int main(int argc, char *argv[])
   double *epj, vepp;
   double kk1, kk2;
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
-  /*int *movingaverage; */
 
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
@@ -3663,7 +3692,7 @@ int main(int argc, char *argv[])
      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) */
 
-  if(mle==1){
+  if(mle>=1){ /* Could be 1 or 2 */
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
   }