version 1.131, 2009/06/20 16:22:47
|
version 1.132, 2009/07/06 08:22:05
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.132 2009/07/06 08:22:05 brouard |
|
Many tings |
|
|
Revision 1.131 2009/06/20 16:22:47 brouard |
Revision 1.131 2009/06/20 16:22:47 brouard |
Some dimensions resccaled |
Some dimensions resccaled |
|
|
Line 597 void replace_back_to_slash(char *s, char
|
Line 600 void replace_back_to_slash(char *s, char
|
} |
} |
} |
} |
|
|
|
char *trimbb(char *out, char *in) |
|
{ /* Trim multiple blanks in line */ |
|
char *s; |
|
s=out; |
|
while (*in != '\0'){ |
|
while( *in == ' ' && *(in+1) == ' ' && *(in+1) != '\0'){ |
|
in++; |
|
} |
|
*out++ = *in++; |
|
} |
|
*out='\0'; |
|
return s; |
|
} |
|
|
int nbocc(char *s, char occ) |
int nbocc(char *s, char occ) |
{ |
{ |
int i,j=0; |
int i,j=0; |
Line 1692 double funcone( double *x)
|
Line 1709 double funcone( double *x)
|
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); |
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ |
if(globpr){ |
if(globpr){ |
fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ |
fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ |
%11.6f %11.6f %11.6f ", \ |
%11.6f %11.6f %11.6f ", \ |
Line 1897 double hessii(double x[], double delta,
|
Line 1914 double hessii(double x[], double delta,
|
int i; |
int i; |
int l=1, lmax=20; |
int l=1, lmax=20; |
double k1,k2; |
double k1,k2; |
double p2[NPARMAX+1]; |
double p2[MAXPARM+1]; /* identical to x */ |
double res; |
double res; |
double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; |
double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; |
double fx; |
double fx; |
Line 1918 double hessii(double x[], double delta,
|
Line 1935 double hessii(double x[], double delta,
|
/*res= (k1-2.0*fx+k2)/delt/delt; */ |
/*res= (k1-2.0*fx+k2)/delt/delt; */ |
res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ |
res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ |
|
|
#ifdef DEBUG |
#ifdef DEBUGHESS |
printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); |
printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); |
fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); |
fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); |
#endif |
#endif |
Line 1944 double hessij( double x[], double delti[
|
Line 1961 double hessij( double x[], double delti[
|
int i; |
int i; |
int l=1, l1, lmax=20; |
int l=1, l1, lmax=20; |
double k1,k2,k3,k4,res,fx; |
double k1,k2,k3,k4,res,fx; |
double p2[NPARMAX+1]; |
double p2[MAXPARM+1]; |
int k; |
int k; |
|
|
fx=func(x); |
fx=func(x); |
Line 2155 void freqsummary(char fileres[], int ia
|
Line 2172 void freqsummary(char fileres[], int ia
|
pos += freq[jk][m][i]; |
pos += freq[jk][m][i]; |
if(pp[jk]>=1.e-10){ |
if(pp[jk]>=1.e-10){ |
if(first==1){ |
if(first==1){ |
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
} |
} |
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
}else{ |
}else{ |
Line 2464 void tricode(int *Tvar, int **nbcode, in
|
Line 2481 void tricode(int *Tvar, int **nbcode, in
|
} |
} |
|
|
for (i=0; i<=cptcode; i++) { /* i=-1 ?*/ |
for (i=0; i<=cptcode; i++) { /* i=-1 ?*/ |
if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariate. In fact ncodemax[j]=2 (dichotom. variables only) but it can be more */ |
if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j |
|
th covariate. In fact |
|
ncodemax[j]=2 |
|
(dichotom. variables only) but |
|
it can be more */ |
} /* Ndum[-1] number of undefined modalities */ |
} /* Ndum[-1] number of undefined modalities */ |
|
|
ij=1; |
ij=1; |
Line 4431 int main(int argc, char *argv[])
|
Line 4452 int main(int argc, char *argv[])
|
double ***mobaverage; |
double ***mobaverage; |
int *indx; |
int *indx; |
char line[MAXLINE], linepar[MAXLINE]; |
char line[MAXLINE], linepar[MAXLINE]; |
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; |
char linetmp[MAXLINE]; |
|
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; |
char pathr[MAXLINE], pathimach[MAXLINE]; |
char pathr[MAXLINE], pathimach[MAXLINE]; |
char **bp, *tok, *val; /* pathtot */ |
char **bp, *tok, *val; /* pathtot */ |
int firstobs=1, lastobs=10; |
int firstobs=1, lastobs=10; |
Line 4635 int main(int argc, char *argv[])
|
Line 4657 int main(int argc, char *argv[])
|
covar=matrix(0,NCOVMAX,1,n); |
covar=matrix(0,NCOVMAX,1,n); |
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ |
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ |
if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; |
if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; |
|
/* where is ncovprod ?*/ |
ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ |
ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/ |
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ |
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ |
npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ |
nforces= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ |
|
npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters like aij*/ |
if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){ |
if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){ |
printf("Too complex model for current IMaCh: npar=(nlstate+ndeath-1)*nlstate*ncovmodel=%d >= %d(MAXPARM) or nlstate=%d >= %d(NLSTATEMAX) or ndeath=%d >= %d(NDEATHMAX) or ncovmodel=(k+age+#of+signs)=%d(NCOVMAX) >= %d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); |
printf("Too complex model for current IMaCh: npar=(nlstate+ndeath-1)*nlstate*ncovmodel=%d >= %d(MAXPARM) or nlstate=%d >= %d(NLSTATEMAX) or ndeath=%d >= %d(NDEATHMAX) or ncovmodel=(k+age+#of+signs)=%d(NCOVMAX) >= %d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); |
fprintf(ficlog,"Too complex model for current IMaCh: %d >=%d(MAXPARM) or %d >=%d(NLSTATEMAX) or %d >=%d(NDEATHMAX) or %d(NCOVMAX) >=%d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); |
fprintf(ficlog,"Too complex model for current IMaCh: %d >=%d(MAXPARM) or %d >=%d(NLSTATEMAX) or %d >=%d(NDEATHMAX) or %d(NCOVMAX) >=%d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); |
Line 4858 run imach with mle=-1 to get a correct t
|
Line 4881 run imach with mle=-1 to get a correct t
|
printf("Comment line\n%s\n",line); |
printf("Comment line\n%s\n",line); |
continue; |
continue; |
} |
} |
|
trimbb(linetmp,line); /* Trims multiple blanks in line */ |
|
for (j=0; line[j]!='\0';j++){ |
|
line[j]=linetmp[j]; |
|
} |
|
|
|
|
for (j=maxwav;j>=1;j--){ |
for (j=maxwav;j>=1;j--){ |
cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
Line 5040 run imach with mle=-1 to get a correct t
|
Line 5068 run imach with mle=-1 to get a correct t
|
cptcovprod--; |
cptcovprod--; |
cutv(strb,stre,strd,'V'); |
cutv(strb,stre,strd,'V'); |
Tvar[i]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3 */ |
Tvar[i]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3 */ |
cptcovage++; /* Sum the number of covariates including ages as a product */ |
cptcovage++; /* Sums the number of covariates including age as a product */ |
Tage[cptcovage]=i; /* Tage[1] =2 */ |
Tage[cptcovage]=i; /* Tage[1] =2 */ |
/*printf("stre=%s ", stre);*/ |
/*printf("stre=%s ", stre);*/ |
} |
} |
Line 5436 Interval (in months) between two waves:
|
Line 5464 Interval (in months) between two waves:
|
} /* Endof if mle==-3 */ |
} /* Endof if mle==-3 */ |
|
|
else{ /* For mle >=1 */ |
else{ /* For mle >=1 */ |
globpr=1;/* debug */ |
globpr=0;/* debug */ |
/* likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone);*/ /* Prints the contributions to the likelihood */ |
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
for (k=1; k<=npar;k++) |
for (k=1; k<=npar;k++) |
printf(" %d %8.5f",k,p[k]); |
printf(" %d %8.5f",k,p[k]); |