create_a_and_d.c

最大似然估计算法

#include "timeseries.h"

data_kernel create_A_and_d(time_series ts, double *periods, int n_periods, int component, int trend) {

	int j,k,n,l;
	int n_par;
	int used_off;
	int offsetStart;
	int nip;

	int *timeIndex;

	double t0;

	data_kernel dk;

	timeIndex = (int *) calloc( (size_t)ts.n_offsets, sizeof(int) );
	used_off = sort_offsets(ts, component, timeIndex);

	nip = trend == 1 ? 2 : 1;

	n_par = nip + used_off + 2 * n_periods;

	dk.name_id = (int *) calloc((size_t)n_par,            sizeof(int));
	dk.alt_id  = (int *) calloc((size_t)n_par,            sizeof(int));

	dk.d      = (double *) calloc((size_t)ts.n_data,           sizeof(double));
	dk.dd     = (double *) calloc((size_t)ts.n_data,           sizeof(double));
	dk.A      = (double *) calloc((size_t)(ts.n_data * n_par), sizeof(double));
	dk.params = (double *) calloc((size_t)n_par,               sizeof(double));
	dk.covar  = (double *) calloc((size_t)(n_par*n_par),       sizeof(double));
	dk.MLE    = (double *) calloc((size_t)1,                   sizeof(double));

	for (j = 0; j < ts.n_data; j++) {
		 dk.d[j]  = ts.data[(j*ts.n_series)+component];
		 dk.dd[j] = ts.formal_error[(j*ts.n_series)+component];
	}

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

	t0 = floor(ts.t[0]);

	offsetStart = nip + n_periods * 2;

	/*******************************/
	/* Form the major name tag ids */
	/*******************************/

	dk.name_id[0] = 0; 
	if (trend == 1) dk.name_id[1] = 1;

	for (l = 0; l < n_periods; l++) {
		dk.name_id[l*2+nip] = 2;
		dk.name_id[l*2+nip+1] = 3;
	}
	for (l = 0; l < used_off; l++) dk.name_id[offsetStart+l] = 4;

	/*****************************/
	/* Form the sub name tag ids */
	/*****************************/
	dk.alt_id[0] = -1;
	if (trend == 1) dk.alt_id[1] = -1;

	for (l = 0; l < n_periods; l++) {
		dk.alt_id[l*2+nip]   = l;
		dk.alt_id[l*2+nip+1] = l;
	}
	for (l = 0; l < used_off; l++) dk.alt_id[offsetStart+l] = timeIndex[l];
	

	for (k = 0; k < ts.n_data; k++) {
		dk.A[k] = 1.0;
		if (trend == 1) dk.A[k + ts.n_data] = ts.t[k] - t0;

		for (l = 0; l  < n_periods; l++) {
			dk.A[k + (l*2+nip)*ts.n_data]   = sin(2.0 * M_PI * (ts.t[k]-t0) / periods[l]);
			dk.A[k + (l*2+nip+1)*ts.n_data] = cos(2.0 * M_PI * (ts.t[k]-t0) / periods[l]);

		}

		for (j = 0; j < used_off; j++) {
			if (k >= timeIndex[j]) 	dk.A[k+(offsetStart+j)*ts.n_data] = 1.0;
			else 			dk.A[k+(offsetStart+j)*ts.n_data] = 0.0;
		}
	}
	free(timeIndex);
	return(dk);
}