cats_simplex.c

最大似然估计算法

#include "timeseries.h" 

double cats_simplex(double *start, time_series ts, data_kernel dk, noise_model nm, options op, double *params_mle, double *cov_mle, int return_cov) {
		
	int    i, j, k, m, ia;
	int    n_vert, ilo, ihi, inhi, n_params, info;
	int    nfun_eval, liwork, lwork, got_answer;
	int    il, iu;

	int    *ifail, *iwork, *ipiv;

	double **p, *mle, *psum, *params, *ptry;
	double *fit, *cov_fit, *data_hat, *work;
	double *Eig_vectors, *Eig_values, *C_unit, *C, *Ctmp;
	double *residuals, *Ioo;
	double *cov_wh, *params_wh, *cov_cn, *params_cn;

	double mle_try, mle_save, min_mle;
	double fac, fac1, fac2;
	double md1, md2, md3, dist;
	double F2, b, wh_mle, cn_mle, kappa;
	double F, dummy1, dummy2, dummy3;
	double F1, F3, F4, cross_partial, second_partial;
	double abstol, a;
	double vl, vu, xo;

	clock_t time1, time2;

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

	n_params = 2;
	n_vert = n_params+1;

	/*******************************************************************************/
	/* 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 = 100 + 6 * dk.n_data + 2 * dk.n_data * dk.n_data; */
	lwork = 35 * dk.n_data;
	work        = (double *) calloc((size_t) lwork,                sizeof(double));

	/* liwork = 30 + 5 * dk.n_data; */
	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);
	/* dsyevd_("V","L",&dk.n_data,Eig_vectors,&dk.n_data,Eig_values,work,&lwork,iwork,&liwork,&info);  */
	/* abstol = 2.0 * dlamch_("S");
	fprintf(stderr, "abstol = %e\n", abstol);
	dsyevx_("V","A","L",&dk.n_data,Ctmp,&dk.n_data,&vl,&vu,&il,&iu,&abstol,&m,Eig_values,Eig_vectors,&dk.n_data,work,&lwork,iwork,ifail,&info);  */
	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) cn_only = %.8f (%.4f)\n", params_wh[dk.n_par], wh_mle, params_cn[dk.n_par], -cn_mle);

	/***********************************/
	/* Set the initial starting values */
	/***********************************/

	params = (double *)  calloc((size_t) n_params, sizeof(double));
	mle    = (double *)  calloc((size_t) n_vert,   sizeof(double));
	psum   = (double *)  calloc((size_t) n_params, sizeof(double));
	ptry   = (double *)  calloc((size_t) n_params, 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));

	if (op.input_method == 1) {
		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];
			}
		}
	} else if (op.input_method == 2) {
		b = params_wh[dk.n_par];
		start[0] = b;
		F = 1.0 / ts.time_span / sec_per_year;
		switch (nm.model) {
			case 'b':
				start[1] = b;
				break;
			case 'f':
				kappa = nm.par[0];
				start[1] = b * sqrt(sec_per_year) * sqrt(kappa*kappa + 4.0 * M_PI * M_PI * F * F) / sqrt(2.0);
				break;
			case 'p':
				kappa = nm.par[0];
				start[1] = b * pow(ts.fs,kappa/4.0) / pow(2.0*M_PI,kappa/2.0) / pow(sec_per_year,kappa/4.0), pow(F,kappa/2.0);
				break;
			default:
				fprintf(stderr, " cats_simplex : unknown covariance model : %s\n", nm.model);
				exit(EXIT_FAILURE);
		}
		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];
			}
		}
	} else if (op.input_method == 3) {
		cats_empirical(ts, dk, nm, op, &dummy1, &dummy2, &dummy3); 
		start[0] = dummy1;
		start[1] = dummy2;
		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];
			}
		}
	} else if (op.input_method == 4) {
		p[0][0] = params_wh[dk.n_par]; 
		p[0][1] = 1e-10;
		p[1][0] = 0.0;
		p[1][1] = params_cn[dk.n_par];
		p[2][0] = 0.75 * p[0][0];
		p[2][1] = 0.75 * p[1][1];

	} else if (op.input_method == 5) {
		xo = params_cn[dk.n_par] * params_wh[dk.n_par] / (params_cn[dk.n_par] + params_wh[dk.n_par]);

		p[0][0] = xo;
		p[0][1] = xo;
		p[1][0] = 0.25 * xo;
		p[1][1] = 0.25 * xo;

		if (-cn_mle > wh_mle) {
			p[2][0] = 0.25 * xo;
			p[2][1] = params_cn[dk.n_par] - 0.25 * params_cn[dk.n_par]*params_cn[dk.n_par] / (params_cn[dk.n_par] +  params_wh[dk.n_par]);
		} else {
			p[2][0] = params_wh[dk.n_par] - 0.25 * params_wh[dk.n_par]*params_wh[dk.n_par] / (params_cn[dk.n_par] +  params_wh[dk.n_par]);
			p[2][1] = 0.25 * xo;

		}
	} else {
		fprintf(stderr, " cats_simplex : unknown method for creating starting parameters : %d\n", op.input_method);
		exit(EXIT_FAILURE);

	}


	/*******************************/
	/* Begin the simplex algorithm */
	/*******************************/

	if (op.speed == 3) {
		for (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j];
	} else if (op.speed == 2) {
		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 (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
	}

	for (k = 0; k < n_vert; k++) {
		for (j = 0; j < n_params; j++) params[j] = p[k][j];
		if (op.speed < 2) {
			cov_scale(C_unit,dk.n_data,params,C);
			linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);
			for (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
		}
		mle[k] = MLE_withline_CN(params, dk.n_data, Eig_values, Eig_vectors, residuals);
	}

	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;

		if (op.speed <= 1) {
			cov_scale(C_unit,dk.n_data,ptry,C);
			linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);
			for (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
		}

		mle_try = MLE_withline_CN(ptry, dk.n_data, Eig_values, Eig_vectors, residuals);