linefit_final.c

最大似然估计算法

#include "timeseries.h"

double linefit_final(double *start, time_series ts, data_kernel dk, noise_model nm, options op, double *params_mle, double *cov_mle) {

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

	int *iwork, *ipiv;

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

	double min_mle, wh_mle, cn_mle, b;
	double F1, F2, F3, F4, second_partial, cross_partial;

	/******************************************************************/
	/* This is for a power-law 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 = 100 + 6 * dk.n_data + 2 * dk.n_data * dk.n_data;
	liwork = 30 + 5 * dk.n_data; */
	lwork = 35 * dk.n_data;
	liwork = 5 * dk.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                                                 */
	/*************************************************************************/

	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);  */
	dsyev_("V","L",&dk.n_data,Eig_vectors,&dk.n_data,Eig_values,work,&lwork,&info);

	free(iwork);
	free(work);

	/***************************************************************/
	/* While this is a "mimimum" we still need to check it against */
	/* the white noise only case and the pl only case              */
	/***************************************************************/

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

	cn_mle = -cn_mle;

	/*****************************************************************/
	/* Now calculate an estimate of the covariance matrix of the MLE */
	/*****************************************************************/

	i = dk.n_par + n_params;

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

	if (op.speed == 3) {
		for (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j];
		for (j = 0; j < dk.n_par; j++) fit[j] = 0.0;
	} else {
		for (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
	}

	F2 = MLE_withline_CN(start, dk.n_data, Eig_values, Eig_vectors, residuals);
	min_mle = -F2;

	/* fprintf(stderr, " min_mle = %.8f wh_mle = %.8f cn_mle = %.8f\n", min_mle, wh_mle, cn_mle); */
	/* fprintf(stderr, " min_mle-wh_mle = %.8f min_mle-cn_mle = %.8f\n", min_mle-wh_mle, min_mle-cn_mle); */
	if (abs(min_mle-wh_mle) < 1e-8) {

		start[0] = params_wh[dk.n_par];
		start[1] = 0.0;

		for (j = 0; j < dk.n_par; j++) {
			params_mle[j] = params_wh[j];
			for (k = 0; k < dk.n_par; k++) cov_mle[j+k*i] = cov_wh[j+k*(dk.n_par+1)];
		}
		params_mle[dk.n_par] = params_wh[dk.n_par];
		params_mle[i-1]    = 0.0;

		cov_mle[dk.n_par + dk.n_par * i] = cov_wh[dk.n_par * (dk.n_par + 2)];
		min_mle = wh_mle;

	} else if (abs(min_mle-cn_mle) < 1e-8) {

		start[0] = 0.0;
		start[1] = params_cn[dk.n_par];
	
		for (j = 0; j < dk.n_par; j++) {
			params_mle[j] = params_cn[j];
			for (k = 0; k < dk.n_par; k++) cov_mle[j+k*i] = cov_cn[j+k*(dk.n_par+1)];
		}
		params_mle[i-1]    = params_cn[dk.n_par];
		params_mle[dk.n_par] = 0.0;
		cov_mle[i * i - 1] = cov_cn[dk.n_par * (dk.n_par + 2)];
		min_mle = cn_mle;
	} else {

		
		/*********************************************/
		/* Invert cov_fit and substitute it into Ioo */
		/*********************************************/

		Ioo = (double *) calloc((size_t)(i*i), sizeof(double) ); 
	
		lwork = dk.n_par;

		ipiv = (int *)    calloc((size_t) lwork, sizeof(int));
		work = (double *) calloc((size_t) lwork, sizeof(double));
	
		dgetrf_(&dk.n_par, &dk.n_par, cov_fit, &dk.n_par, ipiv, &info);
		dgetri_(&dk.n_par, cov_fit, &dk.n_par, ipiv, work, &lwork, &info);

		free(ipiv);
		free(work);
	
		for (j = 0; j < dk.n_par; j++) {
			for (k = 0; k < dk.n_par; k++) Ioo[j+k*i] = -cov_fit[j+k*dk.n_par];
		}

		work   = (double *) calloc( (size_t) i, sizeof(double) );
		params = (double *) calloc( (size_t) i, sizeof(double) );

		for (j = 0; j < dk.n_par; j++) work[j]          = params[j]          = fit[j];
		for (j = 0; j < n_params; j++) work[j+dk.n_par] = params[j+dk.n_par] = start[j]; 
	
		/******************************************************************************/
		/* For the linear parts we already have the covariance matrix so we just      */
		/* need to obtain the parts of the covariance matrix for the non-linear parts */
		/* which is obtained using the Fisher information matrix                      */
		/******************************************************************************/

		for (j = dk.n_par; j < i; j++) {

			params[j] = work[j] + 1e-3 * fabs(work[j]);
			F1 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

			params[j] = work[j] - 1e-3 * fabs(work[j]);
			F3 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

			second_partial = -(F1-2.0*F2+F3)/pow(1e-3*fabs(work[j]),2.0);

			Ioo[j + j * i] = second_partial;
			params[j] = work[j];
		}
			
		for (j = dk.n_par; j < i; j++) {
			for (k = 0; k < j; k++) {
				params[j] = work[j] + 1e-3 / 2.0 * fabs(work[j]);
				params[k] = work[k] + 1e-3 / 2.0 * fabs(work[k]);
				F1 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

				params[j] = work[j] - 1e-3 / 2.0 * fabs(work[j]);
				params[k] = work[k] + 1e-3 / 2.0 * fabs(work[k]);
				F2 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

				params[j] = work[j] + 1e-3 / 2.0 * fabs(work[j]);
				params[k] = work[k] - 1e-3 / 2.0 * fabs(work[k]);
				F3 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

				params[j] = work[j] - 1e-3 / 2.0 * fabs(work[j]);
				params[k] = work[k] - 1e-3 / 2.0 * fabs(work[k]);
				F4 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

				cross_partial = -(F1-F2-F3+F4) / (1e-3 * fabs(work[j]) * 1e-3 * fabs(work[k]));
				Ioo[j + k * i] = Ioo[k + j * i] = cross_partial;
			
				params[j] = work[j];
				params[k] = work[k];
			}
		}

		/************************************/
		/* Store the params into params_mle */
		/************************************/
	
		for (j = 0; j < i; j++) params_mle[j] = work[j];
	
		/**********************/
		/* Calculate inv(Ioo) */
		/**********************/
	
		lwork = i;
		free(work);

		ipiv = (int *)    calloc((size_t) lwork, sizeof(int));
		work = (double *) calloc((size_t) lwork, sizeof(double));
	
		dgetrf_(&i, &i, Ioo, &i, ipiv, &info);
		dgetri_(&i, Ioo, &i, ipiv, work, &lwork, &info);

		free(ipiv);
		free(work);
	
		/***************************/
		/* Store -Ioo into cov_mle */
		/***************************/
	
		for (j = 0; j < i; j++) {
			for (k = 0; k < i; k++) cov_mle[j + k * i] = -Ioo[j + k * i];
		}

		free(Ioo);
		free(params);

	}

	free(Eig_vectors);
	free(Eig_values);
	free(C);
	free(C_unit);
	free(residuals);
	free(fit);
	free(cov_fit);
	free(data_hat);

	free(cov_wh);
	free(params_wh);
	free(cov_cn);
	free(params_cn);

	return(min_mle);
}