version 1.361, 2024/05/12 20:29:32
|
version 1.364, 2024/06/28 12:27:05
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.364 2024/06/28 12:27:05 brouard |
|
* imach.c (Module): fixing some bugs in gnuplot and quantitative variables, but not completely solved |
|
|
|
Revision 1.363 2024/06/28 09:31:55 brouard |
|
Summary: Adding log lines too |
|
|
|
Revision 1.362 2024/06/28 08:00:31 brouard |
|
Summary: 0.99s6 |
|
|
|
* imach.c (Module): s6 errors with age*age (harmless). |
|
|
Revision 1.361 2024/05/12 20:29:32 brouard |
Revision 1.361 2024/05/12 20:29:32 brouard |
Summary: Version 0.99s5 |
Summary: Version 0.99s5 |
|
|
Line 1240 Important routines
|
Line 1251 Important routines
|
- Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities) |
- Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities) |
and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually. |
and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually. |
- printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables |
- printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables |
o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, eĢliminating 1 1 if |
o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, eliminating 1 1 if |
race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless. |
race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless. |
|
|
|
|
Line 1531 char *endptr;
|
Line 1542 char *endptr;
|
long lval; |
long lval; |
double dval; |
double dval; |
|
|
|
/* This for praxis gegen */ |
|
/* int prin=1; */ |
|
double h0=0.25; |
|
double macheps; |
|
double ffmin; |
|
|
#define NR_END 1 |
#define NR_END 1 |
#define FREE_ARG char* |
#define FREE_ARG char* |
#define FTOL 1.0e-10 |
#define FTOL 1.0e-10 |
Line 3007 static void print2() /* print a line of
|
Line 3024 static void print2() /* print a line of
|
/* printf("... after %u function calls ...\n", nf); */ |
/* printf("... after %u function calls ...\n", nf); */ |
/* printf("... including %u linear searches ...\n", nl); */ |
/* printf("... including %u linear searches ...\n", nl); */ |
/* printf("%10d %10d%14.7g",nl, nf, fx); */ |
/* printf("%10d %10d%14.7g",nl, nf, fx); */ |
printf ( "\n" ); |
/* printf ( "\n" ); */ |
printf ( " Linear searches %d", nl ); |
printf ( " Linear searches %d", nl ); |
|
fprintf (ficlog, " Linear searches %d", nl ); |
/* printf ( " Linear searches %d\n", nl ); */ |
/* printf ( " Linear searches %d\n", nl ); */ |
/* printf ( " Function evaluations %d\n", nf ); */ |
/* printf ( " Function evaluations %d\n", nf ); */ |
/* printf ( " Function value FX = %g\n", fx ); */ |
/* printf ( " Function value FX = %g\n", fx ); */ |
printf ( " Function evaluations %d", nf ); |
printf ( " Function evaluations %d", nf ); |
printf ( " Function value FX = %.12lf\n", fx ); |
printf ( " Function value FX = %.12lf\n", fx ); |
|
fprintf (ficlog, " Function evaluations %d", nf ); |
|
fprintf (ficlog, " Function value FX = %.12lf\n", fx ); |
#ifdef DEBUGPRAX |
#ifdef DEBUGPRAX |
printf("n=%d prin=%d\n",n,prin); |
printf("n=%d prin=%d\n",n,prin); |
#endif |
#endif |
if(fx <= fmin) printf(" UNDEFINED "); else printf("%14.7g",log(fx-fmin)); |
/* if(fx <= fmin) printf(" UNDEFINED "); else printf("%14.7g",log(fx-fmin)); */ |
if ( n <= 4 || 2 < prin ) |
if ( n <= 4 || 2 < prin ) |
{ |
{ |
/* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */ |
/* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */ |
for(i=1;i<=n;i++)printf("%14.7g",x[i]); |
for(i=1;i<=n;i++){ |
|
printf(" %14.7g",x[i]); |
|
fprintf(ficlog," %14.7g",x[i]); |
|
} |
/* r8vec_print ( n, x, " X:" ); */ |
/* r8vec_print ( n, x, " X:" ); */ |
} |
} |
printf("\n"); |
printf("\n"); |
|
fprintf(ficlog,"\n"); |
} |
} |
|
|
|
|
Line 3231 L1: /* L1 or try loop */
|
Line 3255 L1: /* L1 or try loop */
|
if (k > 0) *d2 = 0; |
if (k > 0) *d2 = 0; |
} |
} |
#ifdef DEBUGPRAX |
#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 |
#endif |
if (*d2 <= small_windows) *d2 = small_windows; |
if (*d2 <= small_windows) *d2 = small_windows; |
*x1 = x2; fx = fm; |
*x1 = x2; fx = fm; |
Line 3707 mloop:
|
Line 3731 mloop:
|
printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t); |
printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t); |
#endif |
#endif |
/* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */ |
/* 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 |
#ifdef DEBUGPRAX |
printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx); |
printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx); |
#endif |
#endif |
Line 4218 void powell(double p[], double **xi, int
|
Line 4242 void powell(double p[], double **xi, int
|
printf(" + age*age "); |
printf(" + age*age "); |
fprintf(ficlog," + 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) { |
if(Typevar[j]==0) { |
printf(" + V%d ",Tvar[j]); |
printf(" + V%d ",Tvar[j]); |
fprintf(ficlog," + V%d ",Tvar[j]); |
fprintf(ficlog," + V%d ",Tvar[j]); |
Line 6593 void mlikeli(FILE *ficres,double p[], in
|
Line 6617 void mlikeli(FILE *ficres,double p[], in
|
#else /* FLATSUP */ |
#else /* FLATSUP */ |
/* powell(p,xi,npar,ftol,&iter,&fret,func);*/ |
/* powell(p,xi,npar,ftol,&iter,&fret,func);*/ |
/* praxis ( t0, h0, n, prin, x, beale_f ); */ |
/* praxis ( t0, h0, n, prin, x, beale_f ); */ |
int prin=1; |
int prin=4; |
double h0=0.25; |
/* double h0=0.25; */ |
double macheps; |
/* double macheps; */ |
double fmin; |
/* double fmin; */ |
macheps=pow(16.0,-13.0); |
macheps=pow(16.0,-13.0); |
/* #include "praxis.h" */ |
/* #include "praxis.h" */ |
/* Be careful that praxis start at x[0] and powell start at p[1] */ |
/* Be careful that praxis start at x[0] and powell start at p[1] */ |
Line 6606 printf("Praxis Gegenfurtner \n");
|
Line 6630 printf("Praxis Gegenfurtner \n");
|
fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog); |
fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog); |
/* praxis ( ftol, h0, npar, prin, p1, func ); */ |
/* praxis ( ftol, h0, npar, prin, p1, func ); */ |
/* fmin = praxis(1.e-5,macheps, h, n, prin, x, 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"); |
printf("End Praxis\n"); |
#endif /* FLATSUP */ |
#endif /* FLATSUP */ |
|
|
Line 8602 void concatwav(int wav[], int **dh, int
|
Line 8626 void concatwav(int wav[], int **dh, int
|
double ***gradg, ***trgradg; /**< for var eij */ |
double ***gradg, ***trgradg; /**< for var eij */ |
double **gradgp, **trgradgp; /**< for var p point j */ |
double **gradgp, **trgradgp; /**< for var p point j */ |
double *gpp, *gmp; /**< 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 ***p3mat; |
double age,agelim, hf; |
double age,agelim, hf; |
/* double ***mobaverage; */ |
/* double ***mobaverage; */ |
Line 10767 set ter svg size 640, 480\nunset log y\n
|
Line 10791 set ter svg size 640, 480\nunset log y\n
|
/* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */ |
/* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */ |
/* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */ |
/* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */ |
kl++; |
kl++; |
|
/* Problem with quantitative variables TinvDoQresult[nres] */ |
/* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */ |
/* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */ |
sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,lv, kl+1, vlv ); |
sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%lg " ,kl,lv, kl+1, vlv );/* Solved but quantitative must be shifted */ |
kl++; |
kl++; |
if(k <cptcovs && cptcovs>1) |
if(k <cptcovs && cptcovs>1) |
sprintf(gplotcondition+strlen(gplotcondition)," && "); |
sprintf(gplotcondition+strlen(gplotcondition)," && "); |
Line 10782 set ter svg size 640, 480\nunset log y\n
|
Line 10807 set ter svg size 640, 480\nunset log y\n
|
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \ |
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \ |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt ); |
fprintf(ficgp,",\\\n '' "); |
fprintf(ficgp,",\\\n '' "); |
fprintf(ficgp," u %d:(",iagec); |
fprintf(ficgp," u %d:(",iagec); /* Below iyearc should be increades if quantitative variable in the reult line */ |
|
/* $7==6 && $8==2.47 ) && (($9-$10) == 1953 ) ? $12/(1.-$24) : 1/0):7 with labels center not */ |
|
/* but was && $7==6 && $8==2 ) && (($7-$8) == 1953 ) ? $12/(1.-$24) : 1/0):7 with labels center not */ |
fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \ |
fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \ |
iyearc, iagec, offyear, \ |
iyearc, iagec, offyear, \ |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc ); |
Line 12283 double gompertz(double x[])
|
Line 12310 double gompertz(double x[])
|
A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); |
A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); |
} else if (cens[i] == 0){ |
} else if (cens[i] == 0){ |
A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) |
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 |
} else |
printf("Gompertz cens[%d] neither 1 nor 0\n",i); |
printf("Gompertz cens[%d] neither 1 nor 0\n",i); |
/*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ |
/*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ |
Line 15843 Interval (in months) between two waves:
|
Line 15871 Interval (in months) between two waves:
|
flatdir=ivector(1,npar); |
flatdir=ivector(1,npar); |
for (j=1;j<=npar;j++) flatdir[j]=0; |
for (j=1;j<=npar;j++) flatdir[j]=0; |
#endif /*LINMINORIGINAL */ |
#endif /*LINMINORIGINAL */ |
powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); |
/* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */ |
#endif |
/* 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); |
fclose(ficrespow); |
#ifdef LINMINORIGINAL |
#ifdef LINMINORIGINAL |
#else |
#else |
free_ivector(flatdir,1,npar); |
free_ivector(flatdir,1,npar); |
#endif /* LINMINORIGINAL*/ |
#endif /* LINMINORIGINAL*/ |
|
#endif /* POWELL */ |
hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz); |
hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz); |
|
|
for(i=1; i <=NDIM; i++) |
for(i=1; i <=NDIM; i++) |
Line 15984 Please run with mle=-1 to get a correct
|
Line 16020 Please run with mle=-1 to get a correct
|
fprintf(ficlog," + age*age "); |
fprintf(ficlog," + age*age "); |
fprintf(fichtm, "<th>+ age*age</th>"); |
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) { |
if(Typevar[j]==0) { |
printf(" + V%d ",Tvar[j]); |
printf(" + V%d ",Tvar[j]); |
fprintf(ficres," + V%d ",Tvar[j]); |
fprintf(ficres," + V%d ",Tvar[j]); |
Line 16055 Please run with mle=-1 to get a correct
|
Line 16091 Please run with mle=-1 to get a correct
|
fprintf(ficlog," + age*age "); |
fprintf(ficlog," + age*age "); |
fprintf(fichtm, "<th>+ age*age</th>"); |
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) { |
if(Typevar[j]==0) { |
printf(" + V%d ",Tvar[j]); |
printf(" + V%d ",Tvar[j]); |
fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]); |
fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]); |