/* Seek the M and m in fitting fertility data to the */
/* Coale and Trussell's model */
/* Datafile must be given as text with lines of "age [TAB] ASMFR" */
#define _f_nums 2
/* constants above must be modified by the likelihood functions to fit */
#define MAXDATA 200
/* constants above must be modified by the number of observations */
#define REALZERO 0.0E+0
#define REALONE 1.0E+0
#define _ds_lambda_init 10.0
#define _ds_compress 0.5
#define _ds_reverse -1.0
#define _ds_expand 2.0
#define _ds_th 1.0E-15
#define _ds_itermax 100000
#define MAXRAGECLS 50
#include
#include
#include
typedef struct vec {double p[_f_nums+1];} pvec; /* pvec.p[_f_nums] is return value of function */
static pvec _f_parms={0.8, 0.2, 1.0}; /* plausible initial parameters, _M and _m */
static long age[MAXDATA], n_data;
static double fr[MAXDATA], ef[MAXDATA], __M, __m;
static pvec ivec[_f_nums+1]; /* unit vectors used in downhill simplex */
static char filename[80];
static double nf[MAXRAGECLS]={0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.175,0.225,0.275,0.325,0.375,0.421,0.460,0.475,0.477,0.475,0.470,0.465,0.460,0.455,0.449,0.442,0.435,0.428,0.420,0.410,0.400,0.389,0.375,0.360,0.343,0.325,0.305,0.280,0.247,0.207,0.167,0.126,0.087,0.055,0.035,0.021,0.011,0.003};
/* n(a): Age-specific natural fertility by Louis Henry */
static double rn[MAXRAGECLS]={0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.004,0.030,0.060,0.100,0.150,0.200,0.250,0.310,0.370,0.440,0.520,0.600,0.680,0.760,0.830,0.900,0.970,1.040,1.110,1.180,1.250,1.320,1.390,1.460,1.530,1.590,1.640,1.670,1.690,1.700};
/* v(a): Age-specific downward discordance factor from natural fertility */
/*
nf[a]=n(a) and rn[a]=v(a) were presented by Coale and Trussell, 1974.
ASMFR = M*n(a)*exp(m*v(a))
*/
double f(pvec PP);
void sort(pvec *X);
long downhills(void);
long read_data(FILE *fp);
int main(int argc, char *argv[]) {
FILE *ff;
long i, k;
int fnp;
char *fx;
if (argc<=1) {
fputs("Datafile must be given as text with lines of age[TAB]ASMFR\n",stderr);
fputs("Usage: ctmm datafile\n\n The result is obtained in stdout\n",stderr);
return 0;
}
else {
strcpy(filename,argv[1]);
ff=fopen(filename,"rt");
n_data=read_data(ff);
fclose(ff);
if (argc>=3)
_f_parms.p[0]=atof(argv[2]);
if (argc>=4)
_f_parms.p[1]=atof(argv[3]);
}
k=downhills();
if (k==_ds_itermax)
fprintf(stdout,"not convergent within %ld iterations\n",k);
else {
fprintf(stdout,"M= %lf, m= %lf, RMSE=%e\n", ivec[0].p[0], ivec[0].p[1], ivec[0].p[2]);
fprintf(stdout,"age\tObs.ASMFR\tEst.ASMFR\n");
for (i=0; iX[_s_k].p[_f_nums]) {
_s_swap=X[_s_i];
X[_s_i]=X[_s_k];
X[_s_k]=_s_swap;
}
}
}
return;
}
long downhills(void) {
long _ds_i, _ds_j, _ds_iter=0;
double _ds_rratio=1.0;
pvec _ds_centroid, _ds_reflected, _ds_compressed;
for (_ds_i=1; _ds_i<=_f_nums; _ds_i++) {
for (_ds_j=0; _ds_j<_f_nums; _ds_j++) {
ivec[0].p[_ds_j] = _f_parms.p[_ds_j];
if (_ds_i==_ds_j+1)
ivec[_ds_i].p[_ds_j] = ivec[0].p[_ds_j] + _ds_lambda_init;
else
ivec[_ds_i].p[_ds_j] = ivec[0].p[_ds_j];
}
} /* setting unit vectors */
for(_ds_i=0; _ds_i<=_f_nums; _ds_i++) {
ivec[_ds_i].p[_f_nums] = f(ivec[_ds_i]);
#ifdef DDEBUG
fprintf(stderr,"Initial L [%ld] th vector=%lf\n", _ds_i, ivec[_ds_i].p[_f_nums]);
#endif
}
sort(ivec);
while ((_ds_iter < _ds_itermax) && (_ds_rratio > _ds_th)) {
if (ivec[0].p[_f_nums] != 0.0)
_ds_rratio = 1.0 - ivec[0].p[_f_nums]/ivec[_f_nums].p[_f_nums];
/* _ds_rratio = 2.0*(ivec[_f_nums].p[_f_nums]-ivec[0].p[_f_nums])/(fabs(ivec[_f_nums].p[_f_nums])+fabs(ivec[0].p[_f_nums])); */
else
_ds_rratio = 0.0;
/* for the condition of exit loop */
#ifdef DEBUG
fprintf(stderr,"ratio = %lf in %ld th iteration\n", _ds_rratio, _ds_iter);
#endif
if (_ds_rratio <= _ds_th)
break;
/* calculating centroid of the vectors except the worst point */
for (_ds_i=0; _ds_i<_f_nums; _ds_i++) {
_ds_centroid.p[_ds_i] = REALZERO;
for (_ds_j=0; _ds_j<_f_nums; _ds_j++) {
_ds_centroid.p[_ds_i] += ivec[_ds_j].p[_ds_i]/(double)_f_nums;
}
#ifdef XXDEBUG
fprintf(stderr,"centroid [%ld] = %lf in iter. %ld\n", _ds_i, _ds_centroid.p[_ds_i], _ds_iter);
#endif
}
_ds_centroid.p[_f_nums] = f(_ds_centroid);
#ifdef XDEBUG
fprintf(stderr, "this centroid's f() = %lf\n", _ds_centroid.p[_f_nums]);
#endif
for (_ds_i=0; _ds_i<_f_nums; _ds_i++) {
_ds_reflected.p[_ds_i] = 2.0*_ds_centroid.p[_ds_i] - ivec[_f_nums].p[_ds_i];
#ifdef XXDEBUG
fprintf(stderr,"reflected [%ld] = %lf in iter. %ld\n", _ds_i, _ds_reflected.p[_ds_i], _ds_iter);
#endif
}
_ds_reflected.p[_f_nums] = f(_ds_reflected);
#ifdef XDEBUG
fprintf(stderr,"reflected f() = %lf in iter. %ld\n", _ds_reflected.p[_f_nums], _ds_iter);
#endif
if (fabs(_ds_reflected.p[_f_nums]-_ds_centroid.p[_f_nums])<_ds_th) _ds_rratio = 0.0;
if (_ds_reflected.p[_f_nums] < ivec[_f_nums].p[_f_nums]) {
for (_ds_i=0; _ds_i<=_f_nums; _ds_i++)
ivec[_f_nums].p[_ds_i] = _ds_reflected.p[_ds_i];
}
if (ivec[_f_nums].p[_f_nums] <= ivec[0].p[_f_nums]) {
/* when the reflected point becomes the best, the simplex must be expanded */
for (_ds_i=0; _ds_i<_f_nums; _ds_i++)
ivec[_f_nums].p[_ds_i] = 2.0*ivec[_f_nums].p[_ds_i] - _ds_centroid.p[_ds_i];
ivec[_f_nums].p[_f_nums] = f(ivec[_f_nums]);
#ifdef DDEBUG
fprintf(stderr,"In iter. %ld, L [%ld] =%lf\n", _ds_iter, _f_nums, ivec[_f_nums].p[_f_nums]);
#endif
}
else {
if (ivec[_f_nums].p[_f_nums] > ivec[_f_nums-1].p[_f_nums]) {
/* when the reflected point is still the worst, the simplex must be compressed */
for (_ds_i=0; _ds_i<_f_nums; _ds_i++)
_ds_compressed.p[_ds_i] = (ivec[_f_nums].p[_ds_i]+_ds_centroid.p[_ds_i])/2.0;
/* _ds_compressed.p[_ds_i] = (ivec[_f_nums].p[_ds_i]+ivec[0].p[_ds_i])/2.0; */
_ds_compressed.p[_f_nums] = f(_ds_compressed);
if (_ds_compressed.p[_f_nums] < ivec[_f_nums].p[_f_nums]) {
for (_ds_i=0; _ds_i<=_f_nums; _ds_i++)
ivec[_f_nums].p[_ds_i] = _ds_compressed.p[_ds_i];
}
if (ivec[_f_nums].p[_f_nums] > ivec[_f_nums-1].p[_f_nums]) {
/* when still worst, whole simplex must be compressed to the best point */
for (_ds_i=1; _ds_i<=_f_nums; _ds_i++)
for (_ds_j=0; _ds_j<_f_nums; _ds_j++)
ivec[_ds_i].p[_ds_j] = (ivec[0].p[_ds_j] + ivec[_ds_i].p[_ds_j])/2.0;
for(_ds_i=0; _ds_i<=_f_nums; _ds_i++)
ivec[_ds_i].p[_f_nums] = f(ivec[_ds_i]);
}
else {
/* for (_ds_i=1; _ds_i<=_f_nums; _ds_i++)
for (_ds_j=0; _ds_j<_f_nums; _ds_j++)
ivec[_ds_i].p[_ds_j] = (ivec[0].p[_ds_j] + ivec[_ds_i].p[_ds_j])/2.0;
for(_ds_i=0; _ds_i<=_f_nums; _ds_i++)
ivec[_ds_i].p[_f_nums] = f(ivec[_ds_i]); */
}
#ifdef XXDEBUG
fprintf(stderr,"In iter. %ld, L =%lf\n", _ds_iter, ivec[0].p[_f_nums]);
#endif
}
else {
/* for (_ds_i=1; _ds_i<=_f_nums; _ds_i++)
for (_ds_j=0; _ds_j<_f_nums; _ds_j++)
ivec[_ds_i].p[_ds_j] = (ivec[0].p[_ds_j] + ivec[_ds_i].p[_ds_j])/2.0;
for(_ds_i=0; _ds_i<=_f_nums; _ds_i++)
ivec[_ds_i].p[_f_nums] = f(ivec[_ds_i]); */
}
#ifdef XXDEBUG
fprintf(stderr,"In iter. %ld, L =%lf\n", _ds_iter, ivec[0].p[_f_nums]);
#endif
}
sort(ivec);
_ds_iter++;
}
return _ds_iter;
}