white_only.c

来自「最大似然估计算法」· C语言 代码 · 共 68 行

#include "timeseries.h"

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

	int i, j, k, c_id;

	double *fit, *cov, *dy, *data_hat, *residual;
	double wh, dwh, M1, M2, M3;
	double second_partial, work, sum;

	i = dk.n_par;
	dy       = (double *) calloc((size_t) dk.n_data,          sizeof(double));
	fit      = (double *) calloc((size_t) dk.n_par,           sizeof(double));
	cov      = (double *) calloc((size_t)(dk.n_par*dk.n_par), sizeof(double));
	data_hat = (double *) calloc((size_t) dk.n_data,          sizeof(double));
	residual = (double *) calloc((size_t) dk.n_data,          sizeof(double));

	for (j = 0; j < dk.n_data; j++) dy[j] = 1.0;

	linefit_fast(dk.A, dk.d, dy, dk.n_data, dk.n_par, fit, cov, data_hat);

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

	wh = sqrt(sum / (double) dk.n_data);

	M2 = MLE_withline_WH(wh, dk.n_data, residual);

	if (end_flag) {
		dk.MLE[0] = M2;

		work = wh + 1e-3 * fabs(wh);
		M1 = MLE_withline_WH(work, dk.n_data, residual);
		work = wh - 1e-3 * fabs(wh);
		M3 = MLE_withline_WH(work, dk.n_data, residual);

		second_partial = (M1-2.0*M2+M3)/pow(1e-3*fabs(wh),2.0);
		dwh = sqrt(-1.0 / second_partial);

		nm.sigma[0]      = wh;
		nm.dsigma[0]     = dwh; 
		nm.sigma_flag[0] = 2;

		for (j = 0; j < dk.n_par; j++) {
			dk.params[j] = fit[j];
			for (k = 0; k < dk.n_par; k++) {
				dk.covar[j + k * i] = cov[j + k * i] * wh * wh;
			}
		}
	}

	free(dy); 
	free(fit); 
	free(cov); 
	free(data_hat); 
	free(residual);

	return(M2);
}