brent_wrap.c

最大似然估计算法

#include "timeseries.h"

void brent_wrap(time_series ts, data_kernel dk, noise_model nm, int color_only, options op) {

	int i, j, k, c_id;

	double mle;

	double *params_final, *cov_final, *start_value;

	i = dk.n_par + 3;
	params_final  = (double *) calloc((size_t) i,    sizeof(double));
	cov_final     = (double *) calloc((size_t)(i*i), sizeof(double));
	start_value   = (double *) calloc((size_t) 3,    sizeof(double));

	if (nm.model == 'f') {
		start_value[0] = log10(M_PI / 4.0 / ts.time_span);
		start_value[2] = log10(2.0 * M_PI * ts.fs * sec_per_year);
		start_value[1] = 0.5 * (start_value[0] + start_value[2]);
	} else if (nm.model == 'p') {
		start_value[0] = -3;
		start_value[1] = -1;
		start_value[2] =  1;
	}

	c_id = ts.c_id[ts.current_series];

	mle =  brent(start_value, ts, dk, nm, op, color_only, params_final, cov_final);

	fprintf(op.fpout, " Number of Parameters = %d\n", i);
	fprintf(op.fpout, "+%s MLE = %.4f\n", comp_names[c_id], mle);

	for (j = 0; j < dk.n_par; j++) {
		switch (dk.name_id[j]) {
			case 2:
			case 3:
				fprintf(op.fpout, "+%s %13d %s", comp_names[c_id], dk.alt_id[j], par_names[dk.name_id[j]]); 
				break;
			case 4:
				fprintf(op.fpout, "+%s %13.8f %s", comp_names[c_id], ts.t[dk.alt_id[j]], par_names[dk.name_id[j]]);
				break;
			case 0:
			case 1:
			default:
				fprintf(op.fpout, "+%s               %s", comp_names[c_id], par_names[dk.name_id[j]]);
				break;
		}
		fprintf(op.fpout, " %10.4f +- %10.4f\n", params_final[j] * 1000.0, sqrt(cov_final[j + j * i])*1000.0); 
	}
	k = dk.n_par;
	fprintf(op.fpout, "+%s               %s", comp_names[c_id], par_names[5]);
	fprintf(op.fpout, " %10.4f +- %10.4f\n", params_final[k] * 1000.0, sqrt(cov_final[k + k * i])*1000.0); 
	k = dk.n_par+1;
	fprintf(op.fpout, "+%s               %s", comp_names[c_id], par_names[nm.sigma_id]);
	fprintf(op.fpout, " %10.4f +- %10.4f\n", params_final[k] * 1000.0, sqrt(cov_final[k + k * i])*1000.0); 
	fprintf(op.fpout, "+%s               %s", comp_names[c_id], par_names[nm.names[1]]);
	fprintf(op.fpout, " %10.4f\n", params_final[dk.n_par+2]); 


	fprintf(op.fpout, "\n+COVAR\n");
	for (j = 0; j < i; j++) {
		fprintf(op.fpout, "XX");
		for (k = 0; k < i; k++) fprintf(op.fpout," %10.4f", cov_final[j + k * i] * 1000.0 * 1000.0);
		fprintf(op.fpout, "\n");
	}
	fprintf(op.fpout, "-COVAR\n\n");

	free(start_value);
	free(params_final);
	free(cov_final);
}