color_only.c

最大似然估计算法

#include "timeseries.h"

double color_only(time_series ts, data_kernel dk, noise_model nm, options op, int end_flag) {

	int j, k, i, l, lwork, liwork, info;
	int c_id, n_data, n_par;

	int *iwork, *ipiv;

	double a, b, P, N, alpha, beta, min_mle;
        double M1, M2, M3;
        double second_partial;

	double sigma, dsigma;

	double *Eig_vectors, *C_unit, *Eig_values;
        double *fit, *cov, *data_hat;
        double *work1;
        double *residuals;

	clock_t time1, time2;

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

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

	l = 0;	
	Eig_vectors = (double *) calloc((size_t)(n_data*n_data),  sizeof(double));
	Eig_values  = (double *) calloc((size_t) n_data,          sizeof(double));
	C_unit      = (double *) calloc((size_t)(n_data*n_data),  sizeof(double));

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

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

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

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

	time1 = clock();

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

	time2 = clock();

	if (op.verbose) {
		fprintf(op.fpout, " Time taken to create covariance matrix %f seconds", ((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,work1,&lwork,&info);

	time1 = clock();

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

	free(iwork);
	free(work1);

	/**********************************************/
	/* Now estimate the parameters and covariance */
        /**********************************************/


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

        /**************************************************/
        /* Calculate the amplitude of the power-law noise */
        /**************************************************/

	time1 = clock();

	M2 = color_only_mle(&sigma, n_data, Eig_values, Eig_vectors, residuals, 1); 

	time2 = clock();

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

	switch (nm.model) {
		case 'b':
			if (nm.pvec[0] <= 0.0) M2 = -1e10;
			if (nm.pvec[1] >= 1.0) M2 = -1e10;
			if (nm.pvec[2] <= 0.0) M2 = -1e10;
			break;
	}

	dk.MLE[0] = M2;

	nm.sigma[0]      = sigma;

	if (end_flag) {
		b = sigma + 1e-3 * fabs(sigma);
		M1 = color_only_mle(&b, n_data, Eig_values, Eig_vectors, residuals, 0); 

		b = sigma - 1e-3 * fabs(sigma);
		M3 = color_only_mle(&b, n_data, Eig_values, Eig_vectors, residuals, 0); 

		second_partial = (M1-2.0*M2+M3)/pow(1e-3*fabs(sigma),2.0);

		dsigma = sqrt(-1.0 / second_partial);

		nm.dsigma[0]     = dsigma; 
		nm.sigma_flag[0] = 2;

		i = dk.n_par;
		for (j = 0; j < n_par; j++) {
			dk.params[j] = fit[j];
			for (k = 0; k < n_par; k++) {
				dk.covar[j + k * n_par] = cov[j + k * n_par] * sigma * sigma;
			}
		}
	
	}
	free(Eig_vectors);
	free(Eig_values);
	free(C_unit);

	free(residuals);
	free(fit);
	free(cov);
	free(data_hat);

	return(M2);
}