brent_color_only.c

最大似然估计算法

#include "timeseries.h"

#define SHFT(a,b,c,d) (a)=(b); (b) = (c); (c) = (d);

double brent_color_only(time_series ts, data_kernel dk, noise_model nm, options op) {

	int j, jj, k, i, iter, key;

	int imax       = 100;
	int count      = 0;
	int got_answer = 0;

	double a, b, d, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, x, xm;
	double MLE, xmin;
	double m, c;

	double e         = 0.0;
	double zeps      = 1.0e-10;
	double cgold     = 0.3819660;
	double tolerance = 0.01;

	double start_value[3];

	double *value, *cn, *mle;

	key = 0;
	i = 0;
	switch (nm.model) {
		case 'f':
                	start_value[0] = log10(M_PI / 4.0 / ts.time_span);
                	start_value[2] = log10(2.0 * M_PI * ts.fs * sec_per_year);
                	start_value[1] = 0.5 * (start_value[0] + start_value[2]);
			key = 1;
			i = 0;
			break;
		case 'p':
                	start_value[0] = -3;
                	start_value[1] = -1;
                	start_value[2] =  1;
			i = 0;
			break;
		case 's':
			start_value[0] = ts.t[0];
			start_value[2] = ts.t[ts.n_data-1];
                	start_value[1] = 0.75 * start_value[0] + 0.25 * start_value[2];
			tolerance = 1e-8;	
			i = 0;
			break;
		case 't':
			start_value[0] = -2.0;
			start_value[1] = 0.0;
			start_value[2] = 2.0;
			i = 0;
			break;
		case 'b':
			if (nm.pvec_flag[0] == 0) {
				start_value[0] = 0.25;
				start_value[1] = 0.5;  
				start_value[2] = 1.0;
				i = 0;
			} else if (nm.pvec_flag[1] == 0) {
				start_value[0] = 1.0;
				start_value[1] = 2.0;  
				start_value[2] = 4.0;
				i = 1;
			} else {
				start_value[0] = 1.0;
				start_value[1] = 2.0;
				start_value[2] = 4.0;
				i = 2;
			}
			break;
		case 'g':
			if (nm.pvec_flag[0] == 0) {
                		start_value[0] = log10(M_PI / 4.0 / ts.time_span);
                		start_value[2] = log10(2.0 * M_PI * ts.fs * sec_per_year);
                		start_value[1] = 0.5 * (start_value[0] + start_value[2]);
				i = 0;
				key = 1;
			} else {
                		start_value[0] = -3;
                		start_value[1] = -1;
                		start_value[2] =  1;
				i = 1;
			}
			break;
        }


	if (op.verbose) {
		fprintf(op.fpout, " Starting the one-dimensional minimisation : initial ");
		fprintf(op.fpout, " %s = %5.2f\n", par_names[nm.names[i*2+1]], key == 1 ? pow(10.0, start_value[1]) : start_value[1]);
	}
	

	a = (start_value[0] < start_value[2] ? start_value[0] : start_value[2]);
	b = (start_value[0] > start_value[2] ? start_value[0] : start_value[2]);

	x = w = v = start_value[1];
	
	mle        = (double *) calloc((size_t)(200),     sizeof(double));
	value      = (double *) calloc((size_t)(200),     sizeof(double));
	cn         = (double *) calloc((size_t)(200),     sizeof(double));

	j = 0;

	nm.pvec[i] = (key == 1 ? pow(10.0,x) : x); 

	MLE = color_only(ts, dk, nm, op, 0);

	value[j] = x; mle[j] = MLE; cn[j] = nm.sigma[0]; j++;
	
	if (op.verbose) {
		fprintf(op.fpout, " %s  %11.6f ", par_names[nm.names[i*2+1]], key == 1 ? pow(10.0,x) : x); 
		fprintf(op.fpout, "mle = %13.8f ", MLE);
		fprintf(op.fpout, "cn = %12.6f\n", cn[j-1]); 
	}

	fw=fv=fx=-MLE;

	for (iter=0; iter < imax; iter++) {
		xm   = 0.5*(a+b);
		tol2 = 2.0*(tol1=tolerance*fabs(x)+zeps);
		if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
			got_answer = 1;
			xmin = x;
			break;
		}
		if (fabs(e) > tol1) {
			r = (x-w)*(fx-fv);
			q = (x-v)*(fx-fw);
			p = (x-v)*q-(x-w)*r;
			q = 2.0 * (q-r);
			if (q > 0.0) p = -p;
			q = fabs(q);
			etemp = e;
			e = d;
			if (fabs(p) >= fabs(0.5 * q* etemp) || p <= q*(a-x) || p >= q * (b-x)) {
				d = cgold*(e=(x >= xm ? a-x : b-x));
			} else {
				d = p / q;
				u = x+d;
				if (u-a < tol2 || b-u < tol2) {
					d = (xm-x < 0.0 ? -fabs(tol1) : fabs(tol1) );
				}
			}
		} else {
			d = cgold*(e=(x >= xm ? a-x : b-x));
		}

		u = (fabs(d) >= tol1 ? x+d : x + (d < 0.0 ? -fabs(tol1) : fabs(tol1) ));

		if (op.verbose) fprintf(op.fpout, " Next choice of %s = %9.6f\n", par_names[nm.names[i*2+1]], key == 1 ? pow(10.0,u) : u);

		nm.pvec[i] = (key == 1 ? pow(10.0,u) : u); 

		MLE = color_only(ts, dk, nm, op, 0);

		value[j] = u; mle[j] = MLE; cn[j] = nm.sigma[0]; j++;

		if (op.verbose) {
			fprintf(op.fpout, " %s = %11.6f ", par_names[nm.names[i*2+1]], key == 1 ? pow(10.0,u) : u); 
			fprintf(op.fpout, "mle = %13.8f ", MLE);
			fprintf(op.fpout, "cn = %12.6f\n", cn[j-1]); 
		}

		fu = -MLE;

		if (fu <= fx) {
			if (u >= x) a = x; else b = x;
			SHFT(v,w,x,u)
			SHFT(fv,fw,fx,fu)
		} else {
			if (u < x) a = u; else b = u;
			if (fu <= fw || w == x) {
				v = w;
				w = u;
				fv = fw;
				fw = fu;
			} else if (fu <= fv || v == x || v == w) {	
				v = u;
				fv = fu;
			}
		}
	}

	if (got_answer) {

		if (op.verbose) fprintf(op.fpout, " Finished : \n\n");
		for (k = 0; k < j; k++) {
			if (op.verbose) {
				fprintf(op.fpout, " %s = %11.6f ", par_names[nm.names[i*2+1]], key == 1 ? pow(10.0,value[k]) : value[k]); 
				fprintf(op.fpout, "mle = %13.8f ", mle[k]);
				fprintf(op.fpout, "cn = %12.6f\n", cn[k]); 
			}
		}

		nm.pvec[i] = (key == 1 ? pow(10.0,x) : x); 

		nm.pvec_flag[i] = 2;
		nm.pvec_sigma[i] = 0.0;

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

		/***********************************************/
		/* Add full covariance problem here eventually */
		/* This may be a general problem since you     */
		/* would have to do the same thing regardless  */
		/* of the number of parameters                 */
		/***********************************************/
	} else {
		MLE = 0.0;
	}

        free(value);
	free(cn);
	free(mle);

	return(MLE);
}