cats_simplex_stripped.c

最大似然估计算法

#include "timeseries.h" 

/************************************************************************/
/* Stripped down simplex_CN that doesnt calculate the covariance matrix */
/************************************************************************/

double cats_simplex_stripped(double *start, time_series ts, data_kernel dk, noise_model nm, options op, double *params_mle) {
		
	int    i, j, k;
	int    n_vert, ilo, ihi, inhi, n_params, info;
	int    nfun_eval, liwork, lwork, got_answer;

	int    *iwork;

	double **p, *mle, *psum, *params, *ptry;
	double *fit, *cov_fit, *data_hat, *work;
	double *Eig_vectors, *Eig_values, *C_unit, *C;
	double *residuals;
	double *cov_wh, *params_wh, *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;

	time_t time1, time2;
	struct tm *loctime;

	/*****************************************************************/
	/* 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;
	work        = (double *) calloc((size_t) lwork,                sizeof(double));

	liwork = 30 + 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                                                 */
	/*************************************************************************/

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

	time1 = time(NULL);
	loctime = localtime(&time1);

	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); 

	free(iwork);
	free(work);

	time2 = time(NULL);
	loctime = localtime(&time2);

	if (op.fpout != NULL) {
		fprintf(op.fpout, " Time taken to create covariance matrix and compute eigen value and vectors : ");
		fprintf(op.fpout, "%lg seconds\n", difftime(time2,time1) );
	}

	/***********************************************************/
	/* 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));
	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_stripped(dk, C_unit, Eig_values, Eig_vectors, params_cn, op.speed);

	fprintf(stdout, "wh_only = %.8f (%.4f) pl_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 {
		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.fpout != NULL) {
			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);

		if (mle_try < mle[ihi]) {
			mle[ihi] = mle_try;
			for (j = 0; j < n_params; j++) {
				psum[j] += ptry[j]-p[ihi][j];
				p[ihi][j] = ptry[j];
			}
		}

		/************************************************************/

		if (mle_try <= mle[ilo]) {

			fac = 2.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);

			if (mle_try < mle[ihi]) {
				mle[ihi] = mle_try;
				for (j = 0; j < n_params; j++) {
					psum[j] += ptry[j]-p[ihi][j];
					p[ihi][j] = ptry[j];
				}
			}
		} else if (mle_try >= mle[inhi]) {
			mle_save = mle[ihi];

			fac = 0.5;
			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 == 0) {
				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);


			if (mle_try < mle[ihi]) {
				mle[ihi] = mle_try;
				for (j = 0; j < n_params; j++) {
					psum[j] += ptry[j]-p[ihi][j];
					p[ihi][j] = ptry[j];
				}
			}
	
			if (mle_try >= mle_save) {
				for (i = 0; i < n_vert; i++) {
					if (i != ilo) {
						for (j = 0; j < n_params; j++) {
							p[i][j] = psum[j] = 0.5 * (p[i][j] + p[ilo][j]);
						}
						for (j = 0; j < n_params; j++) params[j] = p[i][j];
						if (op.speed == 0) {
							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[i] = MLE_withline_CN(params, dk.n_data, Eig_values, Eig_vectors, residuals);
					}
				}
			}
			nfun_eval += n_params;
			for (j = 0; j < n_params; j++) {
				psum[j] = 0.0;
				for (k = 0; k < n_vert; k++) psum[j] += p[k][j];
			}
		} else --nfun_eval;
	}

	/**************************/
	/* The end of the simplex */
	/**************************/

	if (got_answer) {

		for (j = 0; j < n_params; j++) start[j] = p[ilo][j];

		if (op.speed < 3) {
			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);

			i = dk.n_par + n_params;
			for (j = 0; j < dk.n_par; j++) params_mle[j]          = fit[j];
			for (j = 0; j < n_params; j++) params_mle[j+dk.n_par] = start[j]; 
		} else {

			for (j = 0; j < dk.n_par; j++) params_mle[j]          = 0.0;
			for (j = 0; j < n_params; j++) params_mle[j+dk.n_par] = start[j]; 

		}

		min_mle = MLE_nparam_CN(params_mle,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values,Eig_vectors);

		if (wh_mle < cn_mle) {
			if (-min_mle < cn_mle) {
				if (op.fpout != NULL) { 
					fprintf(op.fpout, " cats_simplex : ");
					fprintf(op.fpout, "MLE (%f) is less ", -min_mle);
					fprintf(op.fpout, "than cn only MLE (%f)\n", cn_mle); 
				}
				for (j = 0; j < dk.n_par; j++) params_mle[j] = params_cn[j];
				params_mle[dk.n_par] = 0.0;
				params_mle[dk.n_par+1] = params_cn[dk.n_par];
				min_mle = -cn_mle;
			}
		} else {
			if (-min_mle < wh_mle) {
				if (op.fpout != NULL) { 
					fprintf(op.fpout, " cats_simplex : ");
					fprintf(op.fpout, "MLE (%f) is less ", -min_mle);
					fprintf(op.fpout, "than wh only MLE (%f)\n", wh_mle); 
				}
				for (j = 0; j < dk.n_par; j++) params_mle[j] = params_wh[j];
				params_mle[dk.n_par] = params_wh[dk.n_par];
				params_mle[dk.n_par+1] = 0.0;
				min_mle = -wh_mle;
			}
		}

	} else {
		fprintf(stderr, " cats_simplex : Minimization did not converge!!\n");

		min_mle = 1e10;
		for (j = 0; j < i; j++) params_mle[j] = 0.0;
		for (j = 0; j < n_params; j++) start[j] = 0.0;
	}
	
	free(params);
	free(mle);
	free(psum);
	free(ptry);
	for (j = 0; j < n_vert; j++) free(p[j]);
	free(p);
	free(Eig_vectors);
	free(C);
	free(C_unit);
	free(Eig_values);
	free(residuals);
	free(fit);
	free(cov_fit);
	free(data_hat);

	return(-min_mle);
}