From 4b9a2caf6d995706e6a5a90b04af1b1ac3c76109 Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Fri, 28 Jun 2024 08:00:31 +0000 Subject: [PATCH] Summary: 0.99s6 * imach.c (Module): s6 errors with age*age (harmless). --- src/imach.c | 55 +++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 41 insertions(+), 14 deletions(-) diff --git a/src/imach.c b/src/imach.c index 3c02ae3..7c335d9 100644 --- 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
\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, "+ 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(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, "+ 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(fichtm, "+ V%d",Tvar[j]); -- 2.43.0