brent_color_white.c

最大似然估计算法

#include "timeseries.h"

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

double brent_color_white(time_series ts, data_kernel dk, noise_model *nm, int n_models, options op, int end_flag) {

	int j, jj, k, i, iter;
	int n_params, lwork, liwork, info, n_data;
	int wh_flag, cn_flag, n_par;

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

	int *iwork, *ipiv;

	double a, b, d, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, x, xm;
	double wh_mle, MLE, xmin;
	double cn_mle, min_mle, wh_sigma, cn_sigma, sum;
	double F1, F2, F3, F4, second_partial, cross_partial;

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

	double *Eig_vectors, *Eig_values, *C_unit, *C;
	double *residuals, *res_wh, *res_cn,  *Ioo, *dy;
	double *cov_wh, *params_wh, *cov_cn, *params_cn;
	double *fit, *cov_fit, *data_hat, *work, *radius;

	double *start, *params, *start_value;
	double *value, *wh, *cn, *mle;

	clock_t time1, time2;

	/*****************************************************************/
	/* This is for a coloured noise plus white noise so n_params = 2 */
	/*****************************************************************/

	n_params = 2;

	wh_flag = (nm[0].model == 'w' ? 0 : 1);
	cn_flag = (nm[0].model == 'w' ? 1 : 0);

	n_data = dk.n_data;
	n_par  = dk.n_par;

	/*******************************************************************************/
	/* First of all generate all the relevant matrices for the rest of the simplex */
	/*******************************************************************************/


	Eig_vectors = (double *) calloc((size_t)(n_data*n_data), sizeof(double));
	C           = (double *) calloc((size_t)(n_data*n_data), sizeof(double));
	C_unit      = (double *) calloc((size_t)(n_data*n_data), sizeof(double));
	Eig_values  = (double *) calloc((size_t) n_data,         sizeof(double));
	residuals   = (double *) calloc((size_t) n_data,         sizeof(double));
	res_wh      = (double *) calloc((size_t) n_data,         sizeof(double));
	res_cn      = (double *) calloc((size_t) n_data,         sizeof(double));

	fit         = (double *) calloc((size_t) n_par,          sizeof(double));
	cov_fit     = (double *) calloc((size_t)(n_par*n_par),   sizeof(double));
	data_hat    = (double *) calloc((size_t) n_data,         sizeof(double));

	lwork = 35 * n_data;
	liwork = 5 * n_data;

	work        = (double *) calloc((size_t) lwork,          sizeof(double));
	iwork       = (int *)    calloc((size_t) liwork,         sizeof(int));


	/*************************************************************************/
	/* Create the unit power law covariance matrix and find the eigen values */
	/* and the eigen vectors                                                 */
	/*************************************************************************/

	time1 = clock();

	create_covariance(nm[cn_flag], ts, op.cov_method, Eig_vectors, 1);

	time2 = clock();

	if (op.verbose) {
		fprintf(op.fpout, " Time taken to create covariance matrix : ");
		fprintf(op.fpout, "%f seconds\n", ((double) (time2-time1)) / CLOCKS_PER_SEC);
	}

	dlacpy_("Full", &n_data, &n_data, Eig_vectors, &n_data, C_unit, &n_data);
	dsyev_("V","L",&n_data,Eig_vectors,&n_data,Eig_values,work,&lwork,&info);


	free(iwork);
	free(work);

	time1 = clock();

	if (op.verbose) {
		fprintf(op.fpout, " Time taken to compute eigen value and vectors : ");
		fprintf(op.fpout, "%lg seconds\n",((double) (time1-time2)) / CLOCKS_PER_SEC );
	}

	/***********************************************************/
	/* Find the white noise only and coloured noise only mle's */
	/***********************************************************/

	dy = (double *) calloc((size_t) n_data, sizeof(double));

	params_wh = (double *) calloc((size_t) n_par,        sizeof(double));
	cov_wh    = (double *) calloc((size_t)(n_par*n_par), sizeof(double));
	params_cn = (double *) calloc((size_t) n_par,        sizeof(double));
	cov_cn    = (double *) calloc((size_t)(n_par*n_par), sizeof(double));

	for (j = 0; j < dk.n_data; j++) dy[j] = 1.0;

	linefit_fast(dk.A, dk.d, dy, n_data, n_par, params_wh, cov_wh, data_hat);

	if (op.speed == 3) {
		for (j = 0; j < n_data; j++) data_hat[j]  = 0.0;
		for (j = 0; j < n_par; j++)  params_wh[j] = 0.0;
	}

	for (sum = 0.0, j = 0; j < n_data; j++) {
		res_wh[j] = dk.d[j] - data_hat[j];
		sum += (res_wh[j]*res_wh[j]);
	}

	wh_sigma = sqrt(sum / (double) n_data);

	wh_mle = MLE_withline_WH(wh_sigma, n_data, res_wh);


	if (op.speed < 3) {
		linefit_all_fast(dk.A,dk.d,C_unit,n_data,n_par,params_cn,cov_cn,data_hat);
		for (j = 0; j < n_data; j++) res_cn[j] = dk.d[j] - data_hat[j];
	} else {
		for (j = 0; j < n_data; j++) res_cn[j] = dk.d[j];
		for (j = 0; j < n_par; j++) params_cn[j] = 0.0;
	}

	cn_mle = color_only_mle(&cn_sigma, n_data, Eig_values, Eig_vectors, res_cn, 1);

	if (op.verbose) {
		fprintf(op.fpout, " wh_only = %.8f (%.4f),",  wh_sigma, wh_mle); 
		fprintf(op.fpout, " cn_only = %.8f (%.4f)\n", cn_sigma, cn_mle);
	}


	/***********************************************************************/
	/* Now Begin the Brent one-dimensional minimisation based on the angle */
	/***********************************************************************/

	start_value  = (double *) calloc((size_t)(3),       sizeof(double));
	
	start_value[0] = 0.0;
	start_value[1] = 45.0;
	start_value[2] = 90.0;

	if (op.verbose) {
		fprintf(op.fpout, " Starting a one-dimensional minimisation : initial angle %5.2f\n", 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)(2),       sizeof(double));
	radius     = (double *) calloc((size_t)(1),       sizeof(double));
	cn         = (double *) calloc((size_t)(200),     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));

	j = 0;

	if (op.speed == 3) {
		for (k = 0; k < n_data; k++) residuals[k] = dk.d[k];
	} else {
		start[0] = tan(x*M_PI/180.0);
		start[1] = 1.0;
		cov_scale(C_unit,n_data,start,C);
		linefit_all_fast(dk.A,dk.d,C,n_data,n_par,fit,cov_fit,data_hat);
		for (k = 0; k < n_data; k++) residuals[k] = dk.d[k] - data_hat[k];
	}

	MLE = exact_radius(n_data, Eig_values, Eig_vectors, residuals, x*M_PI/180.0, radius); 

	value[j] = x; mle[j] = MLE; wh[j] = radius[0] * sin(x*M_PI/180.0); cn[j] = radius[0] * cos(x*M_PI/180.0); j++;
	
	if (op.verbose) {
		fprintf(op.fpout, " Angle = %9.6f ", x); 
		fprintf(op.fpout, "mle = %13.8f ", MLE);
		fprintf(op.fpout, "radius = %12.6f ", radius[0]);
		fprintf(op.fpout, "wh = %12.6f ",  wh[j-1]);
		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 angle = %9.6f\n", u);

		if (op.speed == 3) {
			for (k = 0; k < n_data; k++) residuals[k] = dk.d[k];
		} else if (op.speed <= 1) {
			start[0] = tan(u*M_PI/180.0);
			start[1] = 1.0;
			cov_scale(C_unit,n_data,start,C);
			linefit_all_fast(dk.A,dk.d,C,n_data,n_par,fit,cov_fit,data_hat);
			for (k = 0; k < n_data; k++) residuals[k] = dk.d[k] - data_hat[k];
		}

		MLE = exact_radius(n_data, Eig_values, Eig_vectors, residuals, u*M_PI/180.0, radius); 

		value[j] = u; mle[j] = MLE; wh[j] = radius[0] * sin(u*M_PI/180.0); cn[j] = radius[0] * cos(u*M_PI/180.0); j++;

		if (op.verbose) {
			fprintf(op.fpout, " angle = %9.6f ", u); 
			fprintf(op.fpout, "mle = %13.8f ", MLE);
			fprintf(op.fpout, "radius = %12.6f ", radius[0]);
			fprintf(op.fpout, "wh = %12.6f ",  wh[j-1]);
			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;
			}
		}
	}