|
version 1.360, 2024/04/30 10:59:22
|
version 1.366, 2024/07/02 09:42:10
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.366 2024/07/02 09:42:10 brouard |
| |
Summary: trying clang on Linux |
| |
|
| |
Revision 1.365 2024/06/28 13:53:38 brouard |
| |
* imach.c (Module): fixing some bugs in gnuplot and quantitative variables, but not completely solved |
| |
|
| |
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 |
| |
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 |
Revision 1.360 2024/04/30 10:59:22 brouard |
| Summary: Version 0.99s4 and estimation of std of e.j/e.. |
Summary: Version 0.99s4 and estimation of std of e.j/e.. |
| |
|
|
Line 1228 Important routines
|
Line 1257 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 1457 int countcallfunc=0; /* Count the numbe
|
Line 1486 int countcallfunc=0; /* Count the numbe
|
| int selected(int kvar); /* Is covariate kvar selected for printing results */ |
int selected(int kvar); /* Is covariate kvar selected for printing results */ |
| |
|
| double jmean=1; /* Mean space between 2 waves */ |
double jmean=1; /* Mean space between 2 waves */ |
| double **matprod2(); /* test */ |
double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b); /* test */ |
| |
/* double **matprod2(); *//* test */ |
| double **oldm, **newm, **savm; /* Working pointers to matrices */ |
double **oldm, **newm, **savm; /* Working pointers to matrices */ |
| double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ |
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ |
| double **ddnewms, **ddoldms, **ddsavms; /* for freeing later */ |
double **ddnewms, **ddoldms, **ddsavms; /* for freeing later */ |
|
Line 1504 char optionfilegnuplot[FILENAMELENGTH],
|
Line 1534 char optionfilegnuplot[FILENAMELENGTH],
|
| /* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */ |
/* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */ |
| /* struct timezone tzp; */ |
/* struct timezone tzp; */ |
| /* extern int gettimeofday(); */ |
/* extern int gettimeofday(); */ |
| struct tm tml, *gmtime(), *localtime(); |
|
| |
|
| extern time_t time(); |
/* extern time_t time(); */ /* Commented out for clang */ |
| |
/* struct tm tml, *gmtime(), *localtime(); */ |
| |
|
| |
|
| struct tm start_time, end_time, curr_time, last_time, forecast_time; |
struct tm start_time, end_time, curr_time, last_time, forecast_time; |
| time_t rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */ |
time_t rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */ |
| time_t rlast_btime; /* raw time */ |
time_t rlast_btime; /* raw time */ |
| struct tm tm; |
/* struct tm tm; */ |
| |
struct tm tm, tml; |
| |
|
| char strcurr[80], strfor[80]; |
char strcurr[80], strfor[80]; |
| |
|
|
Line 1519 char *endptr;
|
Line 1551 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 1969 int nboccstr(char *textin, char *chain)
|
Line 2007 int nboccstr(char *textin, char *chain)
|
| /* in="+V7*V4+age*V2+age*V3+age*V4" chain="age" */ |
/* in="+V7*V4+age*V2+age*V3+age*V4" chain="age" */ |
| char *strloc; |
char *strloc; |
| |
|
| int i,j=0; |
int j=0; |
| |
|
| i=0; |
|
| |
|
| strloc=textin; /* strloc points to "^+V7*V4+age+..." in textin */ |
strloc=textin; /* strloc points to "^+V7*V4+age+..." in textin */ |
| for(;;) { |
for(;;) { |
|
Line 2105 int **imatrix(long nrl, long nrh, long n
|
Line 2141 int **imatrix(long nrl, long nrh, long n
|
| } |
} |
| |
|
| /****************** free_imatrix *************************/ |
/****************** free_imatrix *************************/ |
| void free_imatrix(m,nrl,nrh,ncl,nch) |
/* void free_imatrix(m,nrl,nrh,ncl,nch); */ |
| int **m; |
/* int **m; */ |
| long nch,ncl,nrh,nrl; |
/* long nch,ncl,nrh,nrl; */ |
| /* free an int matrix allocated by imatrix() */ |
void free_imatrix(int **m,long nrl, long nrh, long ncl, long nch) |
| { |
/* free an int matrix allocated by imatrix() */ |
| free((FREE_ARG) (m[nrl]+ncl-NR_END)); |
{ |
| free((FREE_ARG) (m+nrl-NR_END)); |
free((FREE_ARG) (m[nrl]+ncl-NR_END)); |
| } |
free((FREE_ARG) (m+nrl-NR_END)); |
| |
} |
| |
|
| /******************* matrix *******************************/ |
/******************* matrix *******************************/ |
| double **matrix(long nrl, long nrh, long ncl, long nch) |
double **matrix(long nrl, long nrh, long ncl, long nch) |
|
Line 2372 values at the three points, fa, fb , and
|
Line 2409 values at the three points, fa, fb , and
|
| double ulim,u,r,q, dum; |
double ulim,u,r,q, dum; |
| double fu; |
double fu; |
| |
|
| double scale=10.; |
/* double scale=10.; */ |
| int iterscale=0; |
/* int iterscale=0; */ |
| |
|
| *fa=(*func)(*ax); /* xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/ |
*fa=(*func)(*ax); /* xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/ |
| *fb=(*func)(*bx); /* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */ |
*fb=(*func)(*bx); /* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */ |
|
Line 2763 double *e; /* used in minfit, don't konw
|
Line 2800 double *e; /* used in minfit, don't konw
|
| static int prin; /* added */ |
static int prin; /* added */ |
| static int n; |
static int n; |
| static double *x; |
static double *x; |
| static double (*fun)(); |
static double (*fun)(double *x); /* New for clang */ |
| |
/* static double (*fun)(); */ |
| /* static double (*fun)(double *x, int n); */ |
/* static double (*fun)(double *x, int n); */ |
| |
|
| /* these will be set by praxis to the global control parameters */ |
/* these will be set by praxis to the global control parameters */ |
|
Line 2988 static void print() /* print a line of
|
Line 3026 static void print() /* print a line of
|
| /* static void print2(int n, double *x, int prin, double fx, int nf, int nl) */ /* print a line of traces */ |
/* static void print2(int n, double *x, int prin, double fx, int nf, int nl) */ /* print a line of traces */ |
| static void print2() /* print a line of traces */ |
static void print2() /* print a line of traces */ |
| { |
{ |
| int i; double fmin=0.; |
int i; /* double fmin=0.; */ |
| |
|
| /* printf("\n"); */ |
/* printf("\n"); */ |
| /* printf("... chi square reduced to ... %20.10e\n", fx); */ |
/* printf("... chi square reduced to ... %20.10e\n", fx); */ |
| /* 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 3219 L1: /* L1 or try loop */
|
Line 3264 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 3695 mloop:
|
Line 3740 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 3760 next:
|
Line 3805 next:
|
| #ifdef NR_SHIFT |
#ifdef NR_SHIFT |
| fx = (*fun)((x-1), n); |
fx = (*fun)((x-1), n); |
| #else |
#else |
| fx = (*fun)(x, n); |
fx = (*fun)(x); |
| #endif |
#endif |
| /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/ |
/* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/ |
| nf++; |
nf++; |
|
Line 4206 void powell(double p[], double **xi, int
|
Line 4251 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 4414 void powell(double p[], double **xi, int
|
Line 4459 void powell(double p[], double **xi, int
|
| #endif |
#endif |
| #ifdef POWELLORIGINAL |
#ifdef POWELLORIGINAL |
| if (t < 0.0) { /* Then we use it for new direction */ |
if (t < 0.0) { /* Then we use it for new direction */ |
| #else |
#else /* Not POWELLOriginal but Brouard's */ |
| if (directest*t < 0.0) { /* Contradiction between both tests */ |
if (directest*t < 0.0) { /* Contradiction between both tests */ |
| printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); |
printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); |
| printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); |
printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); |
| fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); |
fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); |
| fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); |
fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); |
| } |
} |
| if (directest < 0.0) { /* Then we use it for new direction */ |
if (directest < 0.0) { /* Then we use (P0, Pn) for new direction Xi_n or Xi_iBig */ |
| #endif |
#endif |
| #ifdef DEBUGLINMIN |
#ifdef DEBUGLINMIN |
| printf("Before linmin in direction P%d-P0\n",n); |
printf("Before linmin in direction P%d-P0\n",n); |
|
Line 4455 void powell(double p[], double **xi, int
|
Line 4500 void powell(double p[], double **xi, int
|
| xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ |
xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ |
| xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ |
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ |
| } |
} |
| |
|
| |
/* #else */ |
| |
/* for (i=1;i<=n-1;i++) { */ |
| |
/* for (j=1;j<=n;j++) { */ |
| |
/* xi[j][i]=xi[j][i+1]; /\* Standard method of conjugate directions, not Powell who changes the nth direction by p0 pn . *\/ */ |
| |
/* } */ |
| |
/* } */ |
| |
/* for (j=1;j<=n;j++) { */ |
| |
/* xi[j][n]=xit[j]; /\* and this nth direction by the by the average p_0 p_n *\/ */ |
| |
/* } */ |
| |
/* /\* for (j=1;j<=n-1;j++) { *\/ */ |
| |
/* /\* xi[j][1]=xi[j][j+1]; /\\* Standard method of conjugate directions *\\/ *\/ */ |
| |
/* /\* xi[j][n]=xit[j]; /\\* and this nth direction by the by the average p_0 p_n *\\/ *\/ */ |
| |
/* /\* } *\/ */ |
| |
/* #endif */ |
| #ifdef LINMINORIGINAL |
#ifdef LINMINORIGINAL |
| #else |
#else |
| for (j=1, flatd=0;j<=n;j++) { |
for (j=1, flatd=0;j<=n;j++) { |
|
Line 4479 void powell(double p[], double **xi, int
|
Line 4539 void powell(double p[], double **xi, int
|
| free_vector(pt,1,n); |
free_vector(pt,1,n); |
| return; |
return; |
| #endif |
#endif |
| } |
} /* endif(flatd >0) */ |
| #endif |
#endif /* LINMINORIGINAL */ |
| printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); |
printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); |
| fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); |
fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); |
| |
|
|
Line 4539 void powell(double p[], double **xi, int
|
Line 4599 void powell(double p[], double **xi, int
|
| |
|
| int i, ii,j,k, k1; |
int i, ii,j,k, k1; |
| double *min, *max, *meandiff, maxmax,sumnew=0.; |
double *min, *max, *meandiff, maxmax,sumnew=0.; |
| /* double **matprod2(); */ /* test */ |
double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b); /* test */ /* for clang */ |
| double **out, cov[NCOVMAX+1], **pmij(); /* **pmmij is a global variable feeded with oldms etc */ |
/* double **matprod2(); */ /* test */ |
| |
/* double **out, cov[NCOVMAX+1], **pmij(); */ /* **pmmij is a global variable feeded with oldms etc */ |
| |
double **out, cov[NCOVMAX+1], **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate); /* **pmmij is a global variable feeded with oldms etc */ |
| double **newm; |
double **newm; |
| double agefin, delaymax=200. ; /* 100 Max number of years to converge */ |
double agefin, delaymax=200. ; /* 100 Max number of years to converge */ |
| int ncvloop=0; |
int ncvloop=0; |
|
Line 4746 void powell(double p[], double **xi, int
|
Line 4808 void powell(double p[], double **xi, int
|
| int first=0; |
int first=0; |
| double *min, *max, *meandiff, maxmax,sumnew=0.; |
double *min, *max, *meandiff, maxmax,sumnew=0.; |
| /* double **matprod2(); */ /* test */ |
/* double **matprod2(); */ /* test */ |
| double **out, cov[NCOVMAX+1], **bmij(); |
double **out, cov[NCOVMAX+1], **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij); |
| |
/* double **out, cov[NCOVMAX+1], **bmij(); */ /* Deprecated in clang */ |
| double **newm; |
double **newm; |
| double **dnewm, **doldm, **dsavm; /* for use */ |
double **dnewm, **doldm, **dsavm; /* for use */ |
| double **oldm, **savm; /* for use */ |
double **oldm, **savm; /* for use */ |
|
Line 5011 double **pmij(double **ps, double *cov,
|
Line 5074 double **pmij(double **ps, double *cov,
|
| */ |
*/ |
| int ii, j; |
int ii, j; |
| |
|
| double **pmij(); |
double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate); |
| |
/* double **pmij(); */ /* No more for clang */ |
| double sumnew=0.; |
double sumnew=0.; |
| double agefin; |
double agefin; |
| double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */ |
double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */ |
|
Line 5406 double ***hbxij(double ***po, int nhstep
|
Line 5470 double ***hbxij(double ***po, int nhstep
|
| */ |
*/ |
| |
|
| int i, j, d, h, k1; |
int i, j, d, h, k1; |
| double **out, cov[NCOVMAX+1], **bmij(); |
double **out, cov[NCOVMAX+1], **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij); |
| |
/* double **out, cov[NCOVMAX+1], **bmij(); */ /* No more for clang */ |
| double **newm, ***newmm; |
double **newm, ***newmm; |
| double agexact; |
double agexact; |
| /*double agebegin, ageend;*/ |
/*double agebegin, ageend;*/ |
|
Line 6566 void mlikeli(FILE *ficres,double p[], in
|
Line 6631 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 6579 printf("Praxis Gegenfurtner \n");
|
Line 6644 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 8556 void concatwav(int wav[], int **dh, int
|
Line 8621 void concatwav(int wav[], int **dh, int
|
| /************ Variance ******************/ |
/************ Variance ******************/ |
| void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres) |
void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres) |
| { |
{ |
| /** Variance of health expectancies |
/** Computes the matrix of variance covariance of health expectancies e.j= sum_i w_i e_ij where w_i depends of popbased, |
| |
* either cross-sectional or implied. |
| |
* return vareij[i][j][(int)age]=cov(e.i,e.j)=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20 |
| * double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl); |
* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl); |
| * double **newm; |
* double **newm; |
| * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) |
* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) |
|
Line 8573 void concatwav(int wav[], int **dh, int
|
Line 8640 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 p point j nlstate 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 8641 void concatwav(int wav[], int **dh, int
|
Line 8708 void concatwav(int wav[], int **dh, int
|
| fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n"); |
fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n"); |
| fprintf(fichtm,"\n<br>%s <br>\n",digitp); |
fprintf(fichtm,"\n<br>%s <br>\n",digitp); |
| |
|
| varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); |
varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); /* In fact, currently a double */ |
| pstamp(ficresvij); |
pstamp(ficresvij); |
| fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); |
fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); |
| if(popbased==1) |
if(popbased==1) |
|
Line 8710 void concatwav(int wav[], int **dh, int
|
Line 8777 void concatwav(int wav[], int **dh, int
|
| prlim[i][i]=mobaverage[(int)age][i][ij]; |
prlim[i][i]=mobaverage[(int)age][i][ij]; |
| } |
} |
| } |
} |
| /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h. |
/**< Computes the shifted plus (gp) transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h. |
| */ |
*/ |
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */ |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */ |
| /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability |
/**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability |
|
Line 8719 void concatwav(int wav[], int **dh, int
|
Line 8786 void concatwav(int wav[], int **dh, int
|
| for(j=1; j<= nlstate; j++){ |
for(j=1; j<= nlstate; j++){ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| for(i=1, gp[h][j]=0.;i<=nlstate;i++) |
for(i=1, gp[h][j]=0.;i<=nlstate;i++) |
| gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gp[h][j] += prlim[i][i]*p3mat[i][j][h]; /* gp[h][j]= w_i h_pij */ |
| } |
} |
| } |
} |
| /* Next for computing shifted+ probability of death (h=1 means |
/* Next for computing shifted+ probability of death (h=1 means |
| computed over hstepm matrices product = hstepm*stepm months) |
computed over hstepm matrices product = hstepm*stepm months) |
| as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 . |
as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 . |
| */ |
*/ |
| for(j=nlstate+1;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1;j<=nlstate+ndeath;j++){ /* Currently only once for theta plus p.3(age) Sum_i wi pi3*/ |
| for(i=1,gpp[j]=0.; i<= nlstate; i++) |
for(i=1,gpp[j]=0.; i<= nlstate; i++) |
| gpp[j] += prlim[i][i]*p3mat[i][j][1]; |
gpp[j] += prlim[i][i]*p3mat[i][j][1]; |
| } |
} |
|
Line 8748 void concatwav(int wav[], int **dh, int
|
Line 8815 void concatwav(int wav[], int **dh, int
|
| } |
} |
| } |
} |
| |
|
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Still minus */ |
| |
|
| for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ |
for(j=1; j<= nlstate; j++){ /* gm[h][j]= Sum_i of wi * pij = h_p.j */ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| for(i=1, gm[h][j]=0.;i<=nlstate;i++) |
for(i=1, gm[h][j]=0.;i<=nlstate;i++) |
| gm[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gm[h][j] += prlim[i][i]*p3mat[i][j][h]; |
|
Line 8758 void concatwav(int wav[], int **dh, int
|
Line 8825 void concatwav(int wav[], int **dh, int
|
| } |
} |
| /* This for computing probability of death (h=1 means |
/* This for computing probability of death (h=1 means |
| computed over hstepm matrices product = hstepm*stepm months) |
computed over hstepm matrices product = hstepm*stepm months) |
| as a weighted average of prlim. |
as a weighted average of prlim. j is death. gmp[3]=sum_i w_i*p_i3=p.3 minus theta |
| */ |
*/ |
| for(j=nlstate+1;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1;j<=nlstate+ndeath;j++){ /* Currently only once theta_minus p.3=Sum_i wi pi3*/ |
| for(i=1,gmp[j]=0.; i<= nlstate; i++) |
for(i=1,gmp[j]=0.; i<= nlstate; i++) |
| gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
| } |
} |
| /* end shifting computations */ |
/* end shifting computations */ |
| |
|
| /**< Computing gradient matrix at horizon h |
/**< Computing gradient of p.j matrix at horizon h and still for one parameter of vector theta |
| |
* equation 31 and 32 |
| */ |
*/ |
| for(j=1; j<= nlstate; j++) /* vareij */ |
for(j=1; j<= nlstate; j++) /* computes grad p.j(x, over each h) where p.j is Sum_i w_i*pij(x over h) |
| |
* equation 24 */ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; |
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; |
| } |
} |
| /**< Gradient of overall mortality p.3 (or p.j) |
/**< Gradient of overall mortality p.3 (or p.death) |
| */ |
*/ |
| for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */ |
for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* computes grad of p.3 from wi+pi3 grad p.3 (theta) */ |
| gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; |
gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; |
| } |
} |
| |
|
| } /* End theta */ |
} /* End theta */ |
| |
|
| /* We got the gradient matrix for each theta and state j */ |
/* We got the gradient matrix for each theta and each state j of gradg(h]theta][j)=grad(_hp.j(theta) */ |
| trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */ |
trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); |
| |
|
| for(h=0; h<=nhstepm; h++) /* veij */ |
for(h=0; h<=nhstepm; h++) /* veij */ /* computes the transposed of grad (_hp.j(theta)*/ |
| for(j=1; j<=nlstate;j++) |
for(j=1; j<=nlstate;j++) |
| for(theta=1; theta <=npar; theta++) |
for(theta=1; theta <=npar; theta++) |
| trgradg[h][j][theta]=gradg[h][theta][j]; |
trgradg[h][j][theta]=gradg[h][theta][j]; |
| |
|
| for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ |
for(j=nlstate+1; j<=nlstate+ndeath;j++) /* computes transposed of grad p.3 (theta)*/ |
| for(theta=1; theta <=npar; theta++) |
for(theta=1; theta <=npar; theta++) |
| trgradgp[j][theta]=gradgp[theta][j]; |
trgradgp[j][theta]=gradgp[theta][j]; |
| /**< as well as its transposed matrix |
/**< as well as its transposed matrix |
|
Line 8800 void concatwav(int wav[], int **dh, int
|
Line 8869 void concatwav(int wav[], int **dh, int
|
| vareij[i][j][(int)age] =0.; |
vareij[i][j][(int)age] =0.; |
| |
|
| /* Computing trgradg by matcov by gradg at age and summing over h |
/* Computing trgradg by matcov by gradg at age and summing over h |
| * and k (nhstepm) formula 15 of article |
* and k (nhstepm) formula 32 of article |
| * Lievre-Brouard-Heathcote |
* Lievre-Brouard-Heathcote so that for each j, computes the cov(e.j,e.k) (formula 31). |
| |
* for given h and k computes trgradg[h](i,j) matcov (theta) gradg(k)(i,j) into vareij[i][j] which is |
| |
cov(e.i,e.j) and sums on h and k |
| |
* including the covariances. |
| */ |
*/ |
| |
|
| for(h=0;h<=nhstepm;h++){ |
for(h=0;h<=nhstepm;h++){ |
|
Line 8810 void concatwav(int wav[], int **dh, int
|
Line 8882 void concatwav(int wav[], int **dh, int
|
| matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]); |
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]); |
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| for(j=1;j<=nlstate;j++) |
for(j=1;j<=nlstate;j++) |
| vareij[i][j][(int)age] += doldm[i][j]*hf*hf; |
vareij[i][j][(int)age] += doldm[i][j]*hf*hf; /* This is vareij=sum_h sum_k trgrad(h_pij) V(theta) grad(k_pij) |
| |
including the covariances of e.j */ |
| } |
} |
| } |
} |
| |
|
| /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of |
/* Mortality: pptj is p.3 or p.death = trgradgp by cov by gradgp, variance of |
| * p.j overall mortality formula 19 but computed directly because |
* p.3=1-p..=1-sum i p.i overall mortality computed directly because |
| * we compute the grad (wix pijx) instead of grad (pijx),even if |
* we compute the grad (wix pijx) instead of grad (pijx),even if |
| * wix is independent of theta. |
* wix is independent of theta. |
| */ |
*/ |
| matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); |
matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); |
| matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); |
matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); |
| for(j=nlstate+1;j<=nlstate+ndeath;j++) |
for(j=nlstate+1;j<=nlstate+ndeath;j++) |
| for(i=nlstate+1;i<=nlstate+ndeath;i++) |
for(i=nlstate+1;i<=nlstate+ndeath;i++) |
| varppt[j][i]=doldmp[j][i]; |
varppt[j][i]=doldmp[j][i]; /* This is the variance of p.3 */ |
| /* end ppptj */ |
/* end ppptj */ |
| /* x centered again */ |
/* x centered again */ |
| |
|
|
Line 8846 void concatwav(int wav[], int **dh, int
|
Line 8919 void concatwav(int wav[], int **dh, int
|
| hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres); |
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres); |
| for(j=nlstate+1;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1;j<=nlstate+ndeath;j++){ |
| for(i=1,gmp[j]=0.;i<= nlstate; i++) |
for(i=1,gmp[j]=0.;i<= nlstate; i++) |
| gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; /* gmp[j] is p.3 */ |
| } |
} |
| /* end probability of death */ |
/* end probability of death */ |
| |
|
| fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); |
fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); |
| for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
| fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j])); |
fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));/* p.3 (STD p.3) */ |
| for(i=1; i<=nlstate;i++){ |
for(i=1; i<=nlstate;i++){ |
| fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); |
fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); /* wi, pi3 */ |
| } |
} |
| } |
} |
| fprintf(ficresprobmorprev,"\n"); |
fprintf(ficresprobmorprev,"\n"); |
|
Line 9934 void printinggnuplot(char fileresu[], ch
|
Line 10007 void printinggnuplot(char fileresu[], ch
|
| char dirfileres[256],optfileres[256]; |
char dirfileres[256],optfileres[256]; |
| char gplotcondition[256], gplotlabel[256]; |
char gplotcondition[256], gplotlabel[256]; |
| int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0; |
int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0; |
| int lv=0, vlv=0, kl=0; |
/* int lv=0, vlv=0, kl=0; */ |
| |
int lv=0, kl=0; |
| |
double vlv=0; |
| int ng=0; |
int ng=0; |
| int vpopbased; |
int vpopbased; |
| int ioffset; /* variable offset for columns */ |
int ioffset; /* variable offset for columns */ |
|
Line 10732 set ter svg size 640, 480\nunset log y\n
|
Line 10807 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 10747 set ter svg size 640, 480\nunset log y\n
|
Line 10823 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 12248 double gompertz(double x[])
|
Line 12326 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 14605 int main(int argc, char *argv[])
|
Line 14684 int main(int argc, char *argv[])
|
| double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW; |
double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW; |
| double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */ |
double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */ |
| |
|
| |
double stdpercent; /* for computing the std error of percent e.i: e.i/e.. */ |
| double fret; |
double fret; |
| double dum=0.; /* Dummy variable */ |
double dum=0.; /* Dummy variable */ |
| /* double*** p3mat;*/ |
/* double*** p3mat;*/ |
|
Line 15801 Interval (in months) between two waves:
|
Line 15881 Interval (in months) between two waves:
|
| gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1) */ |
gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1) */ |
| #endif |
#endif |
| #ifdef POWELL |
#ifdef POWELL |
| powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); |
#ifdef LINMINORIGINAL |
| #endif |
#else /* LINMINORIGINAL */ |
| |
|
| |
flatdir=ivector(1,npar); |
| |
for (j=1;j<=npar;j++) flatdir[j]=0; |
| |
#endif /*LINMINORIGINAL */ |
| |
/* 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); |
fclose(ficrespow); |
| |
#ifdef LINMINORIGINAL |
| |
#else |
| |
free_ivector(flatdir,1,npar); |
| |
#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 15938 Please run with mle=-1 to get a correct
|
Line 16036 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 16009 Please run with mle=-1 to get a correct
|
Line 16107 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]); |
|
Line 16778 Please run with mle=-1 to get a correct
|
Line 16876 Please run with mle=-1 to get a correct
|
| cptcod= 0; /* To be deleted */ |
cptcod= 0; /* To be deleted */ |
| printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased); |
printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased); |
| fprintf(ficlog, "varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased); |
fprintf(ficlog, "varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased); |
| |
/* Call to varevsij to get cov(e.i, e.j)= vareij[i][j][(int)age]=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20 */ |
| |
/* Depending of popbased which changes the prevalences, either cross-sectional or period */ |
| varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* cptcod not initialized Intel */ |
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* cptcod not initialized Intel */ |
| fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each state\n\ |
fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each state\n\ |
| # (these are weighted average of eij where weights are "); |
# (these are weighted average of eij where weights are "); |
|
Line 16814 Please run with mle=-1 to get a correct
|
Line 16914 Please run with mle=-1 to get a correct
|
| /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ |
/*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ |
| /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ |
/* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ |
| } |
} |
| epj[nlstate+1] +=epj[j]; |
epj[nlstate+1] +=epj[j]; /* epp=sum_j epj = sum_j sum_i w_i e_ij */ |
| } |
} |
| /* printf(" age %4.0f \n",age); */ |
/* printf(" age %4.0f \n",age); */ |
| |
|
| for(i=1, vepp=0.;i <=nlstate;i++) |
for(i=1, vepp=0.;i <=nlstate;i++) /* Variance of total life expectancy e.. */ |
| for(j=1;j <=nlstate;j++) |
for(j=1;j <=nlstate;j++) |
| vepp += vareij[i][j][(int)age]; |
vepp += vareij[i][j][(int)age]; /* sum_i sum_j cov(e.i, e.j) = var(e..) */ |
| fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
| /* vareij[j][i] is the variance of epj */ |
/* vareij[i][j] is the covariance cov(e.i, e.j) and vareij[j][j] is the variance of e.j */ |
| for(j=1;j <=nlstate;j++){ |
for(j=1;j <=nlstate;j++){ |
| fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); |
fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); |
| } |
} |
| /* And proportion of time spent in state j */ |
/* And proportion of time spent in state j */ |
| /* $$ E[r(X,Y)-E(r(X,Y))]^2=[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]' Var(X,Y)[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]$$ */ |
/* $$ E[r(X,Y)-E(r(X,Y))]^2=[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]' Var(X,Y)[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]$$ */ |
| /* \sigma^2_x/\mu_y^2 +\sigma^2_y \mu^2x/\mu_y^4 */ |
/* \frac{\mu_x^2}{\mu_y^2} ( \frac{\sigma^2_x}{\mu_x^2}-2\frac{\sigma_{xy}}{\mu_x\mu_y} +\frac{\sigma^2_y}{\mu_y^2}) */ |
| /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp */ |
/* \frac{e_{.i}^2}{e_{..}^2} ( \frac{\Var e_{.i}}{e_{.i}^2}-2\frac{\Var e_{.i} + \sum_{j\ne i} \Cov e_{.j},e_{.i}}{e_{.i}e_{..}} +\frac{\Var e_{..}}{e_{..}^2})*/ |
| /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstata+1]^4 */ |
/*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp \sigmaxy= */ |
| |
/* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstate+1]^4 */ |
| for(j=1;j <=nlstate;j++){ |
for(j=1;j <=nlstate;j++){ |
| /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */ |
/* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */ |
| fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[nlstate+1]/epj[nlstate+1] + vepp/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1] )); |
/* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */ |
| |
|
| |
for(i=1,stdpercent=0.;i<=nlstate;i++){ /* Computing cov(e..,e.j)=cov(sum_i e.i,e.j)=sum_i cov(e.i, e.j) */ |
| |
stdpercent += vareij[i][j][(int)age]; |
| |
} |
| |
stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]* (vareij[j][j][(int)age]/epj[j]/epj[j]-2.*stdpercent/epj[j]/epj[nlstate+1]+ vepp/epj[nlstate+1]/epj[nlstate+1]); |
| |
/* stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]*(vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[nlstate+1]/epj[nlstate+1]); */ /* Without covariance */ |
| |
/* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[nlstate+1]/epj[nlstate+1] + epj[j]*epj[j]*vepp/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1] )); */ |
| |
fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt(stdpercent)); |
| } |
} |
| fprintf(ficrest,"\n"); |
fprintf(ficrest,"\n"); |
| } |
} |
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src