simplex_color_only.c

最大似然估计算法

#include "timeseries.h" 

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

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

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

	/*****************************************************************/
	/* This is for a coloured noise only with multiple parameters to */
        /* estimate.                                                     */
	/*****************************************************************/

	n_vert = n_params+1;

	/***********************************/
	/* 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));

	switch (nm.model) {
		/***************************************************************/
		/* Now there can be only two models here so far that have more */
		/* than one parameter                                          */
		/***************************************************************/
		case 'g':
			start[0] = 0.5 * ( log10(M_PI / 4.0 / ts.time_span) + log10(2.0 * M_PI * ts.fs * sec_per_year) );
			start[1] = -1;
			break;
		case 'b':
			l = 0;
			if (nm.pvec_flag[0] == 0) start[l++] = 1.0;
			if (nm.pvec_flag[1] == 0) start[l++] = 0.5;
			if (nm.pvec_flag[2] == 0) start[l++] = 2.0;
			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];
		switch (nm.model) {
			case 'g':
				nm.pvec[0] = pow(10.0, params[0]); 
				nm.pvec[1] = params[1];
				break;
			case 'b':
				l = 0;
				if (nm.pvec_flag[0] == 0) nm.pvec[0] = params[l++];
				if (nm.pvec_flag[1] == 0) nm.pvec[1] = params[l++];
				if (nm.pvec_flag[2] == 0) nm.pvec[2] = params[l++];
				break;
		}
		mle[k] = color_only(ts, dk, nm, 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 (op.verbose) {
			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])) {
			got_answer = 1;
			break;
		}

		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;

		switch (nm.model) {
			case 'g':
				nm.pvec[0] = pow(10.0, ptry[0]); 
				nm.pvec[1] = ptry[1];
				break;
			case 'b':
				l = 0;
				if (nm.pvec_flag[0] == 0) nm.pvec[0] = ptry[l++];
				if (nm.pvec_flag[1] == 0) nm.pvec[1] = ptry[l++];
				if (nm.pvec_flag[2] == 0) nm.pvec[2] = ptry[l++];
				break;
		}
		mle_try = color_only(ts, dk, nm, 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;
			switch (nm.model) {
				case 'g':
					nm.pvec[0] = pow(10.0, ptry[0]); 
					nm.pvec[1] = ptry[1];
					break;
				case 'b':
					l = 0;
					if (nm.pvec_flag[0] == 0) nm.pvec[0] = ptry[l++];
					if (nm.pvec_flag[1] == 0) nm.pvec[1] = ptry[l++];
					if (nm.pvec_flag[2] == 0) nm.pvec[2] = ptry[l++];
					break;
			}
			mle_try = color_only(ts, dk, nm, 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;
			switch (nm.model) {
				case 'g':
					nm.pvec[0] = pow(10.0, ptry[0]); 
					nm.pvec[1] = ptry[1];
					break;
				case 'b':
					l = 0;
					if (nm.pvec_flag[0] == 0) nm.pvec[0] = ptry[l++];
					if (nm.pvec_flag[1] == 0) nm.pvec[1] = ptry[l++];
					if (nm.pvec_flag[2] == 0) nm.pvec[2] = ptry[l++];
					break;
			}
			mle_try = color_only(ts, dk, nm, 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];
						switch (nm.model) {
							case 'g':
								nm.pvec[0] = pow(10.0, params[0]); 
								nm.pvec[1] = params[1];
								break;
							case 'b':
								l = 0;
								if (nm.pvec_flag[0] == 0) nm.pvec[0] = params[l++];
								if (nm.pvec_flag[1] == 0) nm.pvec[1] = params[l++];
								if (nm.pvec_flag[2] == 0) nm.pvec[2] = params[l++];
								break;
						}
						mle[i] = color_only(ts, dk, nm, 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];
		}
		switch (nm.model) {
			case 'g':
				nm.pvec[0] = pow(10.0, start[0]); 
				nm.pvec_flag[0] = 2;
				nm.pvec[1] = start[1];
				nm.pvec_flag[1] = 2;
				break;
			case 'b':
				l = 0;
				if (nm.pvec_flag[0] == 0) {
					nm.pvec[0] = start[l++];
					nm.pvec_flag[0] = 2;
				}
				if (nm.pvec_flag[1] == 0) {
					nm.pvec[1] = start[l++];
					nm.pvec_flag[1] = 2;
				}
				if (nm.pvec_flag[2] == 0) {
					nm.pvec[2] = start[l++];
					nm.pvec_flag[2] = 2;
				}
			break;
		}

		min_mle = color_only(ts, dk, nm, op, 1);

	} else {

		fprintf(stderr, " simplex_color_only : 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);
}