#define rint(a) floor(a+0.5)
/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */
#define mytinydouble 1.0e-16
+#define myhugeexpdouble 100. /* exp(100)=2.68*e43
/* #define DEQUAL(a,b) (fabs((a)-(b))<mytinydouble) */
/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/mynrutils_8h-source.html */
/* static double dsqrarg; */
*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]) */
+#ifdef DEBUGLINMIN
+ printf(" DEBUGLINMIN mnbrak *fa=%12.7f *fb=%12.7f, *ax=%12.7f *bx=%12.7f\n", *fa, *fb, *ax, *bx);
+#endif
/* while(*fb != *fb){ /\* *ax should be ok, reducing distance to *ax *\/ */
/* printf("Warning mnbrak *fb = %lf, *bx=%lf *ax=%lf *fa==%lf iter=%d\n",*fb, *bx, *ax, *fa, iterscale++); */
/* double l; */
{
int i;
+ double f2deb;
/* #ifndef MSDOS */
/* double tflin[N]; */
/* #endif */
return (*fun)((tflin-1), n);
#else
/* return (*fun)(tflin, n);*/
+ f2deb= (*fun)(tflin);
+ if(isnan(f2deb)){
+ printf("Isnan flin, f2deb=%14.8e \n",f2deb);
+ }
+ printf("flin, f2deb=%14.8e \n",f2deb);
return (*fun)(tflin);
#endif
}
matprint(" min vectors:",v,n,n);
#endif
/* find step size */
+#ifdef DEBUGPRAX
+ printf(" Find step size\n");
+#endif
s = 0.;
/* for (i=0; i<n; i++) s += x[i]*x[i]; */
for (i=1; i<=n; i++) s += x[i]*x[i];
s = sqrt(s);
- if (dz)
+ if (dz){
t2 = m4*sqrt(fabs(fx)/dmin + s*ldt) + m2*ldt;
- else
+#ifdef DEBUGPRAXSS
+ printf(" additional flin dz s=%14.7f t2=%14.7f fx=%14.7f dmin=%14.7f ldt=%d fm=%14.7f f1=%14.7f\n",s,t2,fx,dmin,ldt,fm,f1);
+#endif
+ }else{
t2 = m4*sqrt(fabs(fx)/(*d2) + s*ldt) + m2*ldt;
+#ifdef DEBUGPRAXSS
+ printf(" additional flin !dz s=%14.7f t2=%14.7f fx=%14.7f *d2=%14.7f ldt=%d fm=%14.7f f1=%14.7f\n",s,t2,fx,*d2,ldt,fm,f1);
+#endif
+ }
s = s*m4 + t;
- if (dz && t2 > s) t2 = s;
- if (t2 < small_windows) t2 = small_windows;
- if (t2 > 0.01*h) t2 = 0.01 * h;
+ if (dz && t2 > s){
+ t2 = s;
+#ifdef DEBUGPRAXSS
+ printf(" additional flin dz&&t2>s s=%14.7f t2=%14.7f \n",s,t2);
+#endif
+ }
+ if (t2 < small_windows){
+ t2 = small_windows;
+#ifdef DEBUGPRAXSS
+ printf(" additional flin t2 < small_windows s=%14.7f t2=%14.7f \n",s,t2);
+#endif
+ }
+ if (t2 > 0.01*h){
+ t2 = 0.01 * h;
+ }
if (fk && f1 <= fm) {
+#ifdef DEBUGPRAXSS
+ printf(" additional flin fk && f1 <= fm X1=%14.7f t2=%14.7f f1=%14.7f fm=%14.7f fk=%d\n",*x1,t2,f1,fm,fk);
+#endif
xm = *x1;
fm = f1;
}
#endif
f2 = flin(x2, j);
#ifdef DEBUGPRAX
+ if(isnan(f2)){
+ printf("Isnan flin, f2=%14.8e \n",f2);
+ printf(" second additional second flin x2=%14.8e x1=%14.8e f0=%14.8e f1=%18.12e dirj=%d\n",x2,*x1,f0,f1,j);
+ f2 = flin(x2, j);
+ }
printf(" additional second flin x2=%16.10e x1=%16.10e f1=%18.12e f0=%18.10e f2=%18.10e fm=%18.10e\n",x2, *x1, f1,f0,f2,fm);
#endif
if (f2 <= fm) {
#ifdef DEBUGPRAX
double d11,d12;
d11=(f1-f0)/(*x1);d12=(f2-f0)/x2;
- printf(" d11=%18.12e d12=%18.12e d11-d12=%18.12e x1-x2=%18.12e (d11-d12)/(x2-(*x1))=%18.12e\n", d11 ,d12, d11-d12, x2-(*x1), (d11-d12)/(x2-(*x1)));
+ printf(" d11=%18.12e d12=%18.12e d11-d12=%18.12e x1-x2=%18.12e (d11-d12)/((*x1)-x2)=%18.12e\n", d11 ,d12, d11-d12, x2-(*x1), (d11-d12)/((*x1)-x2));
printf(" original computing f1=%18.12e *d2=%16.10e f0=%18.12e f1-f0=%16.10e f2-f0=%16.10e\n",f1,*d2,f0,f1-f0, f2-f0);
- double ff1=7.783920622852e+04;
- double f1mf0=9.0344736236e-05;
- *d2 = (f1mf0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2);
+ /* double ff1=7.783920622852e+04; */
+ /* double f1mf0=9.0344736236e-05; */
+ /* *d2 = (f1mf0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); */
/* *d2 = (ff1-f0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); */
- printf(" simpliff computing *d2=%16.10e f1mf0=%18.12e,f1=f0+f1mf0=%18.12e\n",*d2,f1mf0,f0+f1mf0);
+ /* printf(" simpliff computing *d2=%16.10e f1mf0=%18.12e,f1=f0+f1mf0=%18.12e\n",*d2,f1mf0,f0+f1mf0); */
*d2 = ((f1-f0)/ (*x1) - (f2-f0)/x2)/((*x1)-x2);
printf(" overlifi computing *d2=%16.10e\n",*d2);
#endif
#endif
f2 = flin(x2, j); /* x[i]+x2*v[i][j] */
#ifdef DEBUGPRAX
- printf(" after flin f0=%14.8e f1=%14.8e f2=%14.8e fm=%14.8e\n",f0,f1,f2, fm);
+ printf(" after flin f0=%14.8e f1=%14.8e f2=%14.8e fm=%14.8e x2=%14.8e *x1=%14.8e\n",f0,f1,f2, fm, x2, *x1);
#endif
if ((k < nits) && (f2 > f0)) {
#ifdef DEBUGPRAX
/* Convergence test will use last linmin estimation (fret) and compare to former iteration (fp) */
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */
/* New value of last point Pn is not computed, P(n-1) */
+#ifdef FLATSUP
for(j=1;j<=n;j++) {
if(flatdir[j] >0){
printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
/* printf("\n"); */
/* fprintf(ficlog,"\n"); */
}
+#endif /* FLATSUP */
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
s1=0;
for(j=1; j<i; j++){
/* printf("debug1 %d %d ps=%lf exp(ps)=%lf \n",i,j,ps[i][j],exp(ps[i][j])); */
+#ifdef INFINITYORIGINAL
s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+#else
+ if(ps[i][j] >= myhugeexpdouble) /* should return close to HUGE_VAL */
+ s1+=exp(myhugeexpdouble);
+ else
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+#endif
}
for(j=i+1; j<=nlstate+ndeath; j++){
/* printf("debug2 %d %d ps=%lf exp(ps)=%lf \n",i,j,ps[i][j],exp(ps[i][j])); */
+#ifdef INFINITYORIGINAL
s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+#else
+ if(ps[i][j] >= myhugeexpdouble) /* should return close to HUGE_VAL */
+ s1+=exp(myhugeexpdouble);
+ else
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+#endif
}
/* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
- ps[i][i]=1./(s1+1.);
+ ps[i][i]=1./(s1+1.); /* in fact pii */
/* Computing other pijs */
- for(j=1; j<i; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];/* Bug valgrind */
+ for(j=1; j<i; j++){
+#ifdef INFINITYORIGINAL
+ ps[i][j]= exp(ps[i][j])*ps[i][i];/* in fact pij, Bug valgrind */
+#else
+ if(ps[i][j] >= myhugeexpdouble) /* should return close to HUGE_VAL */
+ ps[i][j]= exp(myhugeexpdouble)*ps[i][i];/* in fact pij, Bug valgrind */
+ else
+ ps[i][j]= exp(ps[i][j])*ps[i][i];/* in fact pij, Bug valgrind */
+#endif
+ }
for(j=i+1; j<=nlstate+ndeath; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
+ {
+#ifdef INFINITYORIGINAL
+ ps[i][j]= exp(ps[i][j])*ps[i][i];/* in fact pij, Bug valgrind */
+#else
+ if(ps[i][j] >= myhugeexpdouble) /* should return close to HUGE_VAL */
+ ps[i][j]= exp(myhugeexpdouble)*ps[i][i];/* in fact pij, Bug valgrind */
+ else
+ ps[i][j]= exp(ps[i][j])*ps[i][i];/* in fact pij, Bug valgrind */
+#endif
+ }
/* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
} /* end i */
it should be very low but not zero otherwise the log go to
infinity.
*/
-/* #ifdef INFINITYORIGINAL */
-/* lli=log(out[s1][s2] - savm[s1][s2]); */
-/* #else */
-/* if ((out[s1][s2] - savm[s1][s2]) < mytinydouble) */
-/* lli=log(mytinydouble); */
-/* else */
-/* lli=log(out[s1][s2] - savm[s1][s2]); */
-/* #endif */
+#ifdef INFINITYORIGINAL
lli=log(out[s1][s2] - savm[s1][s2]);
+#else
+ if ((out[s1][s2] - savm[s1][s2]) < mytinydouble)
+ lli=log(mytinydouble);
+ else
+ lli=log(out[s1][s2] - savm[s1][s2]);
+#endif
+ /* lli=log(out[s1][s2] - savm[s1][s2]); */
} else if ( s2==-1 ) { /* alive */
for (j=1,survp=0. ; j<=nlstate; j++)
/* lli= log(survp); */
/* } */
else{
+#ifdef INFINITYORIGINAL
lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+#else
+ if((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2] <mytinydouble){
+ lli=log(mytinydouble);
+ }else{
+ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+ }
+ if (isinf(lli)){
+ printf("isinf (lli)\n");
+ printf("isinfa lli=%14.8e i=%d mw[mi][i]=%d s1=%d s2=%d out[s1][s2]=%14.8e savm[s1][s2]=%14.8e\n",lli,i,mw[mi][i], s1, s2,out[s1][s2],savm[s1][s2] );
+ }
+#endif
/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
}
/*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
/* printf("num[i], i=%d, 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];
+#ifdef DEBUGPRAX
+ if (isinf(lli)){
+ printf("isnan (lli)\n");
+ printf("isnana lli=%14.8e i=%d mw[mi][i]=%d s1=%d s2=%d out[s1][s2]=%14.8e savm[s1][s2]=%14.8e\n",lli,i,mw[mi][i], s1, s2,out[s1][s2],savm[s1][s2] );
+ }
+#endif
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
/* if (lli < log(mytinydouble)){ */
/* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */
s1=s[mw[mi][i]][i];
s2=s[mw[mi+1][i]][i];
if( s2 > nlstate){
+#ifdef INFINITYORIGINAL
lli=log(out[s1][s2] - savm[s1][s2]);
+#else
+ if ((out[s1][s2] - savm[s1][s2]) < mytinydouble)
+ lli=log(mytinydouble);
+ else
+ lli=log(out[s1][s2] - savm[s1][s2]);
+#endif
+ /* lli=log(out[s1][s2] - savm[s1][s2]); */
} else if ( s2==-1 ) { /* alive */
for (j=1,survp=0. ; j<=nlstate; j++)
survp += out[s1][j];
/* bias is positive if real duration
* is higher than the multiple of stepm and negative otherwise.
*/
- if( s2 > nlstate && (mle <5) ){ /* Jackson */
+ if( s2 > nlstate && (mle <5) ){ /* Jackson *//* if p13 is very close to 1, out - savm is close to 0 or negative */
+#ifdef INFINITYORIGINAL
lli=log(out[s1][s2] - savm[s1][s2]);
+#else
+ if ((out[s1][s2] - savm[s1][s2]) < mytinydouble)
+ lli=log(mytinydouble);
+ else
+ lli=log(out[s1][s2] - savm[s1][s2]);
+#endif
} else if ( s2==-1 ) { /* alive */
for (j=1,survp=0. ; j<=nlstate; j++)
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \
int firstpass, int lastpass, int stepm, int weightopt, char model[])
-{ /* Some frequencies as well as proposing some starting values */
+{
+ /* Some frequencies as well as proposing some starting values */
/* Frequencies of any combination of dummy covariate used in the model equation */
- int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0, s1=1, s2=1;
+ int i, m, jk, j1, ibool, z1,j, nj, nl, k, iv, jj=0, s1=1, s2=1;
int iind=0, iage=0;
int mi; /* Effective wave */
int first;
/* For that combination of covariates j1 (V4=1 V3=0 for example), we count and print the frequencies in one pass */
for (iind=1; iind<=imx; iind++) { /* For each individual iind */
- bool=1;
+ ibool=1;
if(j !=0){
if(anyvaryingduminmodel==0){ /* If All fixed covariates */
if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=1+age+%s\n",Tvaraff[z1],NCOVMAX, cptcoveff, model);
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]){ /* for combination j1 of covariates */
/* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */
- bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
- /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", */
- /* bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),*/
+ ibool=0; /* ibool should be equal to 1 to be selected, one covariate value failed */
+ /* printf("ibool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", */
+ /* ibool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),*/
/* j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
} /* Onlyf fixed */
} /* cptcoveff > 0 */
} /* end any */
}/* end j==0 */
- if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */
+ if (ibool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */
/* for(m=firstpass; m<=lastpass; m++){ */
for(mi=firstpass; mi<wav[iind];mi++){ /* mi=1; mi<wav[iind]For each wave */
m=mw[mi][iind];
if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's
value is -1, we don't select. It differs from the
constant and age model which counts them. */
- bool=0; /* not selected */
+ ibool=0; /* not selected */
}else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
/* i1=Tvaraff[z1]; */
/* i2=TnsdVar[i1]; */
/* i4=covar[i1][iind]; */
/* if(i4 != i3){ */
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) { /* Bug valgrind */
- bool=0;
+ ibool=0;
}
}
}
}/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */
} /* end j==0 */
- /* bool =0 we keep that guy which corresponds to the combination of dummy values */
- if(bool==1){ /*Selected */
+ /* ibool =0 we keep that guy which corresponds to the combination of dummy values */
+ if(ibool==1){ /*Selected */
/* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
and mw[mi+1][iind]. dh depends on stepm. */
agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
/* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
}
}else{
- bool=1;
- }/* end bool 2 */
+ ibool=1;
+ }/* end ibool 2 */
} /* end m */
/* for (z1=1; z1<= nqfveff; z1++) { /\* Quantitative variables, calculating mean *\/ */
/* idq[z1]=idq[z1]+weight[iind]; */
/* meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /\* Computes mean of quantitative with selected filter *\/ */
/* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /\* *weight[iind];*\/ /\* Computes mean of quantitative with selected filter *\/ */
/* } */
- } /* end bool */
+ } /* end ibool */
} /* end iind = 1 to imx */
/* prop[s][age] is fed for any initial and valid live state as well as
freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
We still use firstpass and lastpass as another selection.
*/
- int i, m, jk, j1, bool, z1,j, iv;
+ int i, m, jk, j1, ibool, z1,j, iv;
int mi; /* Effective wave */
int iage;
double agebegin; /*, ageend;*/
fprintf(ficlog,"Prevalence combination of varying and fixed dummies %d\n",j1);
int ierroragemin=0;int ierroragemax=0;
for (i=1; i<=imx; i++) { /* Each individual */
- bool=1;
+ ibool=1;
/* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
for(mi=firstpass; mi<wav[i];mi++){ /* mi=1 For this wave too look where individual can be counted V4=0 V3=0 */
m=mw[mi][i];
if( Fixed[Tmodelind[z1]]==1){
iv= Tvar[Tmodelind[z1]];/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality */
- bool=0;
+ ibool=0;
}else if( Fixed[Tmodelind[z1]]== 0) /* fixed */
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) {
- bool=0;
+ ibool=0;
}
}
- if(bool==1){ /* Otherwise we skip that wave/person */
+ if(ibool==1){ /* Otherwise we skip that wave/person */
agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
/* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
if(m >=firstpass && m <=lastpass){
} /* end valid statuses */
} /* end selection of dates */
} /* end selection of waves */
- } /* end bool */
+ } /* end ibool */
} /* end wave */
} /* end individual */
if(ierroragemin>=1 || ierroragemax>=1){
}
L=0.0;
/* agegomp=AGEGOMP; */
- /* for (i=0; i<=imx; i++)
- if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/
-
+ /* for (i=1; i<=imx; i++) */
+ /* printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d wav=%d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]); */
for (i=1;i<=imx ; i++) {
- /* mu(a)=mu(agecomp)*exp(teta*(age-agegomp))
+ /* mu(a)=mu(agegomp)*exp(teta*(age-agegomp))
mu(a)=x[1]*exp(x[2]*(age-agegomp)); x[1] and x[2] are per year.
* L= Product mu(agedeces)exp(-\int_ageexam^agedc mu(u) du ) for a death between agedc (in month)
* and agedc +1 month, cens[i]=0: log(x[1]/YEARM)
if (cens[i] == 1){
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(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+ A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
+ +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) {*/ /* ??? */
L=L+A*weight[i];
- /* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/
+#ifdef DEBUGGOMP
+ printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);
+#endif
}
}
-
- /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/
-
+#ifdef DEBUGGOMP
+ printf("x1=%2.9lf x2=%2.9lf x3=%2.9lf L=%lf\n",x[1],x[2],x[3],L);
+#endif
return -2*L*num/sump;
}
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");
+ 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
git://henry.ined.fr / .git / commitdiff