cats_simplex_wrap.c

最大似然估计算法

#include "timeseries.h"

void cats_simplex_wrap(time_series ts, data_kernel dk, noise_model nm, noise_model wn, options op) {

	int i, j, k, c_id;

	double mle;

	double *params_mle, *cov_mle, *start;

	i = dk.n_par + 2;
	params_mle = (double *) calloc((size_t) i,    sizeof(double));
	cov_mle    = (double *) calloc((size_t)(i*i), sizeof(double));
	start      = (double *) calloc((size_t) 2,    sizeof(double));

	k = op.input_method;
	/* if (nm.got_sigma && wn.got_sigma) {
		start[0] = wn.sigma[wn.current_field];
		start[1] = nm.sigma[nm.current_field];
		op.input_method = 1; 
	} */

	c_id = ts.c_id[ts.current_series];

	/* mle = cats_simplex(start, ts, dk, nm, op, params_mle, cov_mle, 1); */
	mle = brent_angle(ts, dk, nm, op, params_mle, cov_mle, 1);
	op.input_method = k;

	fprintf(op.fpout, " Number of Parameters = %d\n", i);
	fprintf(op.fpout, "+%s MLE = %.8f\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_mle[j] * 1000.0, sqrt(cov_mle[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_mle[k] * 1000.0, sqrt(cov_mle[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_mle[k] * 1000.0, sqrt(cov_mle[k + k * i])*1000.0); 

	if (nm.n_pvec > 0) {
		fprintf(op.fpout, "+%s FIXED         %s", comp_names[c_id], par_names[nm.names[1]]);
		fprintf(op.fpout, " %11.6f\n", nm.pvec[0]);
	}


	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_mle[j + k * i] * 1000.0 * 1000.0);
		fprintf(op.fpout, "\n");
	}
	fprintf(op.fpout, "-COVAR\n\n");

	free(start);
	free(params_mle);
	free(cov_mle);
}