general_outer_simplex.c

最大似然估计算法

#include "timeseries.h" 

double general_outer_simplex(time_series ts, data_kernel dk, noise_model *nm, int n_models, options op) {
		
	int    i, j, k, l;
	int    n_vert, ilo, ihi, inhi, info;
	int    nfun_eval, got_answer, n_params;
	int    vflag;

	double **p, *mle, *psum, *params, *ptry, *start;

	double mle_try, mle_save, min_mle;
	double fac, fac1, fac2;
	double md1, md2, md3, dist;

	vflag = op.verbose;
	op.verbose = (vflag == 2) ? 1 : 0;
	
	/***************************************/
	/* This is the general case so we will */
	/* estimate the number of parameters   */
	/***************************************/

	for (j = 0, n_params = 0; j < n_models; j++) {
		for (k = 0; k < nm[j].n_pvec; k++) {
			if (nm[j].pvec_flag[k] == 0) n_params++;
		}
	}

	n_vert = n_params+1;
	if (op.verbose) fprintf(stderr, " General_outer_simplex : n_params = %d\n", n_params);

	/***********************************/
	/* Set the initial starting values */
	/***********************************/

	ptry   = (double *)  calloc((size_t) n_params, sizeof(double));
	psum   = (double *)  calloc((size_t) n_params, sizeof(double));
	start  = (double *)  calloc((size_t) n_params, sizeof(double));
	params = (double *)  calloc((size_t) n_params, sizeof(double));
	mle    = (double *)  calloc((size_t) n_vert,   sizeof(double));

	p      = (double **) calloc((size_t)(n_params*n_vert), sizeof(double *));
	for (j = 0; j < n_vert; j++) p[j] = (double *) calloc((size_t) n_params, sizeof(double));

	for (l = 0, j = 0; j < n_models; j++) {
		for (k = 0; k < nm[j].n_pvec; k++) {
			if (nm[j].pvec_flag[k] == 0) {
				switch (nm[j].model) {
					case 'g':
						if (k == 0) start[l++] = 0.5 * ( log10(M_PI / 4.0 / ts.time_span) + log10(2.0 * M_PI * ts.fs * sec_per_year) );
						else start[l++] = -1.0;
						break;
					case 'b':
						if (k == 0) start[l++] = 1.0;
						if (k == 1) start[l++] = 0.5;
						if (k == 2) start[l++] = 2.0;
						break;
					case 'p':
						start[l++] = -1.0;
						break;
					case 'f':
						start[l++] = 0.5 * ( log10(M_PI / 4.0 / ts.time_span) + log10(2.0 * M_PI * ts.fs * sec_per_year) );
						break;
					case 's':
						start[l++] = (ts.t[0] + ts.t[ts.n_data]) / 2.0;
						break;
					case 't':
						start[l++] = -0.1;
						break;
				}
			}
		}
	}
			
	for (k = 0; k < n_vert; k++) {
		for (j = 0; j < n_params; j++) {
			p[k][j] = start[j];
			if (k == j+1) p[k][j] = op.delta * start[j];
		}
	}

	/*******************************/
	/* Begin the simplex algorithm */
	/*******************************/

	for (k = 0; k < n_vert; k++) {
		for (j = 0; j < n_params; j++) params[j] = p[k][j];
		sort_params_to_models(nm, n_models, params, n_params,0);
		mle[k] = general_inner_simplex(ts, dk, nm, n_models, op, 0);
		mle[k] = -mle[k];
	}

	for (j = 0; j < n_params; j++) {
		psum[j] = 0.0;
		for (k = 0; k < n_vert; k++) psum[j] += p[k][j];
	}
		

	nfun_eval = 0;
	got_answer = 0;
	while ((double) nfun_eval <= op.tol[0]) {

		if (vflag) {
			for (k = 0; k < n_vert; k++) {
				fprintf(op.fpout, "%2d %3d %12.6f ", k, nfun_eval,  -mle[k]);
				for (j = 0; j < n_params; j++) fprintf(op.fpout, "%12.8f ", p[k][j]);
				fprintf(op.fpout, "\n");
			}
			fprintf(op.fpout, "\n");
			fflush(op.fpout);
		}
			
		ilo = 0;

		ihi = mle[0] > mle[1] ? (inhi=1,0) : (inhi=0,1);

		for (i = 0; i < n_vert; i++) {
			if (mle[i] <= mle[ilo]) ilo = i;
			if (mle[i] > mle[ihi]) {
				inhi = ihi;
				ihi = i;
			} else if (mle[i] > mle[inhi] && i != ihi) inhi = i;
		}

		md1 = md2 = md3 = 0.0;
		for (k = 1; k < n_vert; k++) {
			for (j = 0; j < n_params; j++) {
				dist = fabs((p[k][j] - p[0][j]) / p[0][j]); 
				md1 = dist > md1 ? dist : md1;
				dist = fabs(p[k][j] - p[0][j]);
				md2 = dist > md2 ? dist : md2;
			}
			dist = fabs((mle[k] - mle[0]) / mle[0]);
			md3 = dist > md3 ? dist : md3;
		}

		if ((md1 <= op.tol[1] || md2 <= op.tol[3]) && (md3 <= op.tol[2])) {
                        if (vflag) {
                                fprintf(op.fpout, " Outer_simplex : Stopping criteria met.\n");
                                fprintf(op.fpout, " Outer_simplex : Stopping criteria :(%g <= %g ||  ", md1, op.tol[1]);
                                fprintf(op.fpout, " %g <= %g) &&  ", md2, op.tol[3]);
                                fprintf(op.fpout, " (%g <= %g)\n", md3, op.tol[2]);
                        }

			got_answer = 1;
			break;
		} else {
                        if (vflag) {
                                fprintf(op.fpout, " Outer_simplex : Stopping criteria :(%g <= %g ||  ", md1, op.tol[1]);
                                fprintf(op.fpout, " %g <= %g) &&  ", md2, op.tol[3]);
                                fprintf(op.fpout, " (%g <= %g)\n", md3, op.tol[2]);
                        }
		}

		nfun_eval += 2;

		/************************************************************/
		/* first extrapolate by a factor -1 through the face of the */
		/* simplex across from the high point                       */
		/************************************************************/

		fac  = -1.0;
		fac1 = (1.0 - fac) / (double) n_params;
		fac2 = fac1 - fac;

		for (j = 0 ; j < n_params; j++) ptry[j] = psum[j]*fac1 - p[ihi][j]*fac2;

		sort_params_to_models(nm, n_models, ptry, n_params,0);
		mle_try = general_inner_simplex(ts, dk, nm, n_models, op, 0);
		mle_try = -mle_try;

		if (mle_try < mle[ihi]) {
			mle[ihi] = mle_try;
			for (j = 0; j < n_params; j++) {
				psum[j] += ptry[j]-p[ihi][j];
				p[ihi][j] = ptry[j];
			}
		}

		/************************************************************/

		if (mle_try <= mle[ilo]) {

			fac = 2.0;
			fac1 = (1.0 - fac) / (double) n_params;
			fac2 = fac1 - fac;

			for (j = 0 ; j < n_params; j++) ptry[j] = psum[j]*fac1 - p[ihi][j]*fac2;
			sort_params_to_models(nm, n_models, ptry, n_params,0);
			mle_try = general_inner_simplex(ts, dk, nm, n_models, op, 0);
			mle_try = -mle_try;

			if (mle_try < mle[ihi]) {
				mle[ihi] = mle_try;
				for (j = 0; j < n_params; j++) {
					psum[j] += ptry[j]-p[ihi][j];
					p[ihi][j] = ptry[j];
				}
			}
		} else if (mle_try >= mle[inhi]) {
			mle_save = mle[ihi];

			fac = 0.5;
			fac1 = (1.0 - fac) / (double) n_params;
			fac2 = fac1 - fac;
	
			for (j = 0 ; j < n_params; j++) ptry[j] = psum[j]*fac1 - p[ihi][j]*fac2;
			sort_params_to_models(nm, n_models, ptry, n_params,0);
			mle_try = general_inner_simplex(ts, dk, nm, n_models, op, 0);
			mle_try = -mle_try;

			if (mle_try < mle[ihi]) {
				mle[ihi] = mle_try;
				for (j = 0; j < n_params; j++) {
					psum[j] += ptry[j]-p[ihi][j];
					p[ihi][j] = ptry[j];
				}
			}
	
			if (mle_try >= mle_save) {
				for (i = 0; i < n_vert; i++) {
					if (i != ilo) {
						for (j = 0; j < n_params; j++) {
							p[i][j] = psum[j] = 0.5 * (p[i][j] + p[ilo][j]);
						}
						for (j = 0; j < n_params; j++) params[j] = p[i][j];
						sort_params_to_models(nm, n_models, params, n_params,0);
						mle[i] = general_inner_simplex(ts, dk, nm, n_models, op, 0);
						mle[i] = -mle[i];
					}
				}
			}
			nfun_eval += n_params;
			for (j = 0; j < n_params; j++) {
				psum[j] = 0.0;
				for (k = 0; k < n_vert; k++) psum[j] += p[k][j];
			}
		} else --nfun_eval;
	}

	/**************************/
	/* The end of the simplex */
	/**************************/

	if (got_answer) {

		for (j = 0; j < n_params; j++) start[j] = p[ilo][j];

		sort_params_to_models(nm, n_models, start, n_params, 1);
		min_mle = general_inner_simplex(ts, dk, nm, n_models, op, 1);

	} else {

		fprintf(stderr, " general_outer_simplex : Could not converge on a solution\n");

		min_mle = -1e10;
		for (j = 0; j < n_params; j++) start[j] = 0.0;
	}
	
	free(ptry);
	free(start);
	free(params);
	free(mle);
	free(psum);
	for (j = 0; j < n_vert; j++) free(p[j]);
	free(p);

	return(min_mle);
}