cats_empirical.c

最大似然估计算法

#include "timeseries.h"

void cats_empirical(time_series ts, data_kernel dk, noise_model nm, options op, double *sig_wh, double *sig_cn, double *min_cn) {

	int i, j, k, n;
	int N, idt;

	double T, Dk, c, gamma, degree, dt;
	double sum, mean, std, sig_dx;
	double Fs, tmp, wh_mle;
	double spectral_index, rms;

	double *b, *cc, *dd, *r;
	double *params_wh, *cov_wh;

	Fs = ts.fs; 

	dd = (double *) calloc((size_t) (dk.n_data-1), sizeof(double));
	r  = (double *) calloc((size_t)  dk.n_data,    sizeof(double));

	/*********************/
	/* Get the residuals */
	/*********************/

	i = dk.n_par + 1;
	params_wh = (double *) calloc((size_t) i,    sizeof(double));
	cov_wh    = (double *) calloc((size_t)(i*i), sizeof(double));

	wh_mle = est_line_WH_only(dk, op, params_wh, cov_wh);
	for (j = 0; j < dk.n_data; j++) {
		r[j] = 0.0;
		for (k = 0; k < dk.n_par; k++) r[j] += dk.A[j + k * dk.n_data] * params_wh[k];
		r[j] -= dk.d[j];
	}
	rms = params_wh[dk.n_par];

	free(params_wh);
	free(cov_wh);

	/*****************************/
	/* Get the daily differences */
	/*****************************/

	n = 0;
	for (j = 1; j < dk.n_data; j++) {
		idt = ts.index[j] - ts.index[j-1];
		if (idt == 1) {
			dd[n] = (r[j] - r[j-1]);
			n++;
		}
	}

	sum = 0.0;
	for (j = 0; j < n; j++) sum += dd[j];
	mean = sum / (double) n;

	sum = 0.0;
	for (j = 0; j < n; j++) sum += (dd[j]-mean)*(dd[j]-mean);
	std = sqrt(sum /  (double) (n-1));
	sig_dx = std;

	sum = 0.0;
	for (j = 0; j < dk.n_data; j++) sum += r[j];
	mean = sum /  (double) dk.n_data;

	sum = 0.0;
	for (j = 0; j < dk.n_data; j++) sum += (r[j]-mean)*(r[j]-mean);
	std = sqrt(sum / (double) (dk.n_data-1));

	T = ts.time_span * sec_per_year;
	N = ts.n_full; 

	if (nm.model == 'p') {

		spectral_index = nm.pvec[0];

		Dk = 2.0 * pow(2.0 * M_PI,spectral_index) * pow(sec_per_year,spectral_index / 2.0);

		if (spectral_index < -1.0) {
			gamma = 0.5718 * spectral_index * spectral_index + 1.0826 * spectral_index + 1.5838;
		} else {
        		gamma = -0.0681 * spectral_index + 1.0;
		}
		if (spectral_index == -1.0) {
        		c = Dk * (log(Fs / 2.0) - log(1.0/T)) / sqrt(Fs);
		} else {
        		c = Dk * ( pow(Fs / 2.0,spectral_index+1.0) - pow(1.0/T,spectral_index+1.0) ) /
 				pow(Fs,spectral_index / 2.0 + 1.0) / (spectral_index + 1.0);
		}
	
		c *= gamma;
	
		degree = (spectral_index + 2.0) / 2.0;

		dt = 1.0 / Fs / sec_per_year;

		b  = (double *)calloc((size_t)N, sizeof(double));
		cc = (double *)calloc((size_t)N, sizeof(double));
	
		b[0] = 1.0;
		for (j = 0; j < N-1; j++)  b[j+1] = ((double) j - degree) / ((double) (j+1));
	
		cc[0] = b[0];
		for (j = 1; j < N; j++) cc[j] = cc[j-1] * b[j];
	
		sum = 0.0;
		for (j = 0; j < N; j++) sum += cc[j] *cc[j];
	
		free(b);
		free(cc);
		free(dd);
		free(r);
	
		/*********************************************************************/
		/* Now find the minimum amount of PL noise given the RMS noise value */
		/*********************************************************************/
	
		n = (int) ceil( log((double)N) / log(2.0) );	
		k = (int) pow(2.0, n);
	
		(*min_cn) = sqrt(2.0 * pow((double) k, spectral_index) * rms * rms * pow(Fs,-spectral_index/2.0) / Dk);
	
		/**********************************/
		/* Now estimate sig_wh and sig_cn */
		/**********************************/
	
		tmp = (sig_dx*sig_dx * c - std * std * pow(dt,-spectral_index/2.0) * sum);
		tmp /= (2 * c - pow(dt,-spectral_index/2.0) * sum);
	
		if (tmp >= 0.0) {
			(*sig_wh) = sqrt( tmp );
			if (std > sqrt(tmp))  (*sig_cn) = sqrt((std*std - (*sig_wh)*(*sig_wh)) / c);
			else (*sig_cn) = (*min_cn);
		} else { 
			/* white_noise is  very small */
			(*sig_cn) = sqrt((std*std) / c);
			/*****************************************************************************************/
			/* instead of giving a value of zero for sig_wh give a value one order of magniute below */
			/* the possible level based on sig_dx                                                    */
			/*****************************************************************************************/
	
			(*sig_wh) = sqrt( std * std * sqrt(dt) * sum / 20.0 / c );
		}
	} else {
		fprintf(stderr, " cats_empirical : not derived for model (%c) yet...\n", nm.model);
		exit(EXIT_FAILURE);
	}	
	
	
}