brent_angle.c

最大似然估计算法

#include "timeseries.h"

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

double brent_angle(time_series ts, data_kernel dk, noise_model nm, options op, double *params_mle, double *cov_mle, int return_cov) {

	int j, jj, k, i, iter;
	int n_params, lwork, liwork, info;

	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;
	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.01;

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

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

	clock_t time1, time2;

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

        n_params = 2;

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

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

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

        lwork = 35 * dk.n_data;
        work        = (double *) calloc((size_t) lwork,                sizeof(double));

        liwork = 5 * dk.n_data;
        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, ts, op.cov_method, Eig_vectors, 1);

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

        if (op.verbose) fprintf(op.fpout, " work(1)  = %f\n", work[0]);
        if (op.verbose) fprintf(op.fpout, " info     = %d\n", info);

        free(iwork);
        free(work);

        time2 = clock();

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

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

        i = dk.n_par + 1;
        cov_wh    = (double *) calloc((size_t)(i*i), sizeof(double));
        cov_cn    = (double *) calloc((size_t)(i*i), sizeof(double));
        params_wh = (double *) calloc((size_t) i,    sizeof(double));
        params_cn = (double *) calloc((size_t) i,    sizeof(double));

        wh_mle = est_line_WH_only(dk, op, params_wh, cov_wh);

        cn_mle = cats_CN_only(dk, C_unit, Eig_values, Eig_vectors, params_cn, cov_cn, op.speed);

        if (op.verbose) {
		fprintf(op.fpout, " wh_only = %.8f (%.4f),",  params_wh[dk.n_par],  wh_mle); 
		fprintf(op.fpout, " cn_only = %.8f (%.4f)\n", params_cn[dk.n_par], -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));
	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     = (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;

	if (op.speed == 3) {
		for (k = 0; k < dk.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,dk.n_data,start,C);
		linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);
		for (k = 0; k < dk.n_data; k++) residuals[k] = dk.d[k] - data_hat[k];
	}

	MLE = exact_radius(dk.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]*1000.0);
		fprintf(op.fpout, "wh = %12.6f ",  wh[j-1]*1000.0);
		fprintf(op.fpout, "cn = %12.6f\n", cn[j-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 angle = %9.6f\n", u);

		if (op.speed == 3) {
			for (k = 0; k < dk.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,dk.n_data,start,C);
			linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);
			for (k = 0; k < dk.n_data; k++) residuals[k] = dk.d[k] - data_hat[k];
		}

		MLE = exact_radius(dk.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]*1000.0);
			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;
			}