Summary: 0.99s6

* imach.c (Module): s6 errors with age*age (harmless).

diff --git a/src/imach.c b/src/imach.c
index 3c02ae394fdcb24598625901093cee0e87f57625..7c335d95463813c2a5d6137f6466fb9866f1f740 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,18 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.361  2024/05/12 20:29:32  brouard
+  Summary: Version 0.99s5
+
+  * src/imach.c Version 0.99s5 In fact, the covariance of total life
+  expectancy e.. with a partial life expectancy e.j is high,
+  therefore the complete matrix of variance covariance has to be
+  included in the formula of the standard error of the proportion of
+  total life expectancy spent in a specific state:
+  var(X/Y)=mu_x^2/mu_y^2*(sigma_x^2/mu_x^2 -2
+  sigma_xy/mu_x/mu_y+sigma^2/mu_y^2).  Also an error with mle=-3
+  made the program core dump. It is fixed in this version.
+
   Revision 1.360  2024/04/30 10:59:22  brouard
   Summary: Version 0.99s4 and estimation of std of e.j/e..
 
@@ -1519,6 +1531,12 @@ char *endptr;
 long lval;
 double dval;
 
+/* This for praxis gegen */
+  /* int prin=1; */
+  double h0=0.25;
+  double macheps;
+  double ffmin;
+
 #define NR_END 1
 #define FREE_ARG char*
 #define FTOL 1.0e-10
@@ -3219,7 +3237,7 @@ L1:  /* L1 or try loop */
       if (k > 0) *d2 = 0;
    }
 #ifdef DEBUGPRAX
-   printf(" bebe end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
+   printf(" bebe end of min x1 might be very wrong x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
 #endif
    if (*d2 <= small_windows) *d2 = small_windows;
    *x1 = x2; fx = fm;
@@ -3695,7 +3713,7 @@ mloop:
    printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t); 
 #endif
    /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */
-   minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global */
+   minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global it seems that fx doesn't correspond to f(s=*x1) */
 #ifdef DEBUGPRAX
    printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx); 
 #endif
@@ -4206,7 +4224,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
        printf("  + age*age  ");
        fprintf(ficlog,"  + age*age  ");
     }
-    for(j=1;j <=ncovmodel-2;j++){
+    for(j=1;j <=ncovmodel-2-nagesqr;j++){
       if(Typevar[j]==0) {
        printf("  +      V%d  ",Tvar[j]);
        fprintf(ficlog,"  +      V%d  ",Tvar[j]);
@@ -6581,10 +6599,10 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub
 #else  /* FLATSUP */
 /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/
 /*   praxis ( t0, h0, n, prin, x, beale_f ); */
-  int prin=1;
-  double h0=0.25;
-  double macheps;
-  double fmin;
+  /* int prin=1; */
+  /* double h0=0.25; */
+  /* double macheps; */
+  /* double fmin; */
   macheps=pow(16.0,-13.0);
 /* #include "praxis.h" */
   /* Be careful that praxis start at x[0] and powell start at p[1] */
@@ -6594,7 +6612,7 @@ printf("Praxis Gegenfurtner \n");
 fprintf(ficlog, "Praxis  Gegenfurtner\n");fflush(ficlog);
 /* praxis ( ftol, h0, npar, prin, p1, func ); */
   /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */
-  fmin = praxis(ftol,macheps, h0, npar, prin, p, func);
+  ffmin = praxis(ftol,macheps, h0, npar, prin, p, func);
 printf("End Praxis\n");
 #endif  /* FLATSUP */
 
@@ -8590,7 +8608,7 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
    double ***gradg, ***trgradg; /**< for var eij */
    double **gradgp, **trgradgp; /**< for var p point j */
    double *gpp, *gmp; /**< for var p point j */
-   double **varppt; /**< for var e.. nlstate+1 to nlstate+ndeath */
+   double **varppt; /**< for var p.3 p.death nlstate+1 to nlstate+ndeath */
    double ***p3mat;
    double age,agelim, hf;
    /* double ***mobaverage; */
@@ -12271,7 +12289,8 @@ double gompertz(double x[])
         A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
        } else if (cens[i] == 0){
        A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
-         +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+         +log(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+       /* +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); */  /* To be seen */
       } else
         printf("Gompertz cens[%d] neither 1 nor 0\n",i);
       /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
@@ -15831,8 +15850,16 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
   flatdir=ivector(1,npar); 
   for (j=1;j<=npar;j++) flatdir[j]=0; 
 #endif /*LINMINORIGINAL */
-     powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
-#endif  
+    /* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */
+  /* double h0=0.25; */
+  macheps=pow(16.0,-13.0);
+  printf("Praxis Gegenfurtner mle=%d\n",mle);
+  fprintf(ficlog, "Praxis  Gegenfurtner mle=%d\n", mle);fflush(ficlog);
+   /* ffmin = praxis(ftol,macheps, h0, npar, prin, p, gompertz); */
+  /* For the Gompertz we use only two parameters */
+  int _npar=2;
+   ffmin = praxis(ftol,macheps, h0, _npar, 4, p, gompertz);
+  printf("End Praxis\n");
     fclose(ficrespow);
 #ifdef LINMINORIGINAL
 #else
@@ -15972,7 +15999,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
       fprintf(ficlog,"  + age*age  ");
       fprintf(fichtm, "<th>+ age*age</th>");
     }
-    for(j=1;j <=ncovmodel-2;j++){
+    for(j=1;j <=ncovmodel-2-nagesqr;j++){
       if(Typevar[j]==0) {
        printf("  +      V%d  ",Tvar[j]);
        fprintf(ficres,"  +      V%d  ",Tvar[j]);
@@ -16043,7 +16070,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
        fprintf(ficlog,"  + age*age  ");
        fprintf(fichtm, "<th>+ age*age</th>");
       }
-      for(j=1;j <=ncovmodel-2;j++){
+      for(j=1;j <=ncovmodel-2-nagesqr;j++){
        if(Typevar[j]==0) {
          printf("  +      V%d  ",Tvar[j]);
          fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);