brent.c

最大似然估计算法

#include "timeseries.h"

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

double brent(double *start_value, time_series ts, data_kernel dk, noise_model nm, options op, int color_only, double *params_final, double *cov_final) {

	int j, jj, k, i, iter;

	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 min_cn, wh_mle, MLE, xmin, wh_rms, dummy1, dummy2;
	double m, c, delta, sx, sx2, sy, sxy;

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

	double *start, *params_mle, *cov_mle, *params, *cov;
	double *value, *wh, *cn, *mle, *par;

	if (op.verbose) {
		fprintf(op.fpout, " Starting the one-dimensional minimisation : initial ");
		fprintf(op.fpout, " %s = %5.2f\n", par_names[nm.key_num[0]], nm.model == 'f' ? 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];
	
	start      = (double *) calloc((size_t)(3),       sizeof(double));
	mle        = (double *) calloc((size_t)(200),     sizeof(double));
	value      = (double *) calloc((size_t)(200),     sizeof(double));
	wh         = (double *) calloc((size_t)(200),     sizeof(double));
	cn         = (double *) calloc((size_t)(200),     sizeof(double));

	params_mle = (double *) calloc((size_t)(dk.n_par+2),                sizeof(double));
	cov_mle    = (double *) calloc((size_t)((dk.n_par+2)*(dk.n_par+2)), sizeof(double));

	params     = (double *) calloc((size_t)(dk.n_par+1),                sizeof(double));
	cov        = (double *) calloc((size_t)((dk.n_par+1)*(dk.n_par+1)), sizeof(double));

	j = 0;

	nm.par[0] = (nm.model == 'f' ? pow(10.0,x) : x); 
	if (color_only) {
		par = (double *) calloc((size_t)(dk.n_par+1), sizeof(double));
		MLE =  cats_CN_only_stripped2(ts, dk, nm, op, par);
		for (k = 0; k < dk.n_par; k++) params_mle[k] = par[k];
		params_mle[dk.n_par] = 0.0;
		params_mle[dk.n_par+1] = par[dk.n_par];
		free(par);
	} else {
		/* MLE = cats_simplex(start, ts, dk, nm, op, params_mle, cov_mle,0); */
		MLE = brent_angle(ts, dk, nm, op, params_mle, cov_mle,0);
	}

	value[j] = x; mle[j] = MLE; wh[j] = params_mle[dk.n_par]; cn[j] = params_mle[dk.n_par+1]; j++;
	
	if (op.verbose) {
		fprintf(op.fpout, " %s = %11.6f ", par_names[nm.key_num[0]], nm.model == 'f' ? pow(10.0,x) : x); 
		fprintf(op.fpout, "mle = %13.8f ", MLE);
		fprintf(op.fpout, "wh = %12.6f ",  params_mle[dk.n_par]*1000.0);
		fprintf(op.fpout, "cn = %12.6f\n", params_mle[dk.n_par+1]*1000.0); 
	}

	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.key_num[0]], nm.model == 'f' ? pow(10.0,u) : u);

		nm.par[0] = (nm.model == 'f' ? pow(10.0,u) : u); 
		if (color_only) {
			par = (double *) calloc((size_t)(dk.n_par+1), sizeof(double));
			MLE =  cats_CN_only_stripped2(ts, dk, nm, op, par);
			for (k = 0; k < dk.n_par; k++) params_mle[k] = par[k];
			params_mle[dk.n_par] = 0.0;
			params_mle[dk.n_par+1] = par[dk.n_par];
			free(par);
		} else {
			/* MLE = cats_simplex(start, ts, dk, nm, op, params_mle, cov_mle,0); */
			MLE = brent_angle(ts, dk, nm, op, params_mle, cov_mle,0);
		}

		value[j] = u; mle[j] = MLE; wh[j] = params_mle[dk.n_par]; cn[j] = params_mle[dk.n_par+1]; j++;

		if (op.verbose) {
			fprintf(op.fpout, " %s = %11.6f ", par_names[nm.key_num[0]], nm.model == 'f' ? pow(10.0,u) : u); 
			fprintf(op.fpout, "mle = %13.8f ", MLE);
			fprintf(op.fpout, "wh = %12.6f ",  wh[j-1]*1000.0);
			fprintf(op.fpout, "cn = %12.6f\n", cn[j-1]*1000.0); 
		}

		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.key_num[0]], nm.model == 'f' ? pow(10.0,value[k]) : value[k]); 
				fprintf(op.fpout, "mle = %13.8f ", mle[k]);
				fprintf(op.fpout, "wh = %12.6f ",  wh[k]*1000.0);
				fprintf(op.fpout, "cn = %12.6f\n", cn[k]*1000.0); 
			}
			if (value[k] == xmin) {
				start[0] = wh[k];
				start[1] = cn[k];
			}
		}

		nm.par[0] = (nm.model == 'f' ? pow(10.0,x) : x); 
		MLE = linefit_final(start, ts, dk, nm, op, params_mle, cov_mle);

		i = dk.n_par+2;
		for (j = 0; j < i; j++) params_final[j] = params_mle[j];
		params_final[i] = nm.model == 'f' ? pow(10.0,x) : x;

		for (j = 0; j < i; j++) {
			for (k = 0; k < i; k++) {
				cov_final[j + (i+1)*k] = cov_mle[j + i * k];
			}
		}
		cov_final[i * (i+2)] = 1e-6;


	} else {
		MLE = 0.0;
	}

	free(start);
	free(params_mle);
	free(cov_mle);
	free(params);
	free(cov);
        free(value);
	free(wh);
	free(cn);
	free(mle);

	return(MLE);
}