calculate_fs.c

最大似然估计算法

#include "timeseries.h"

double calculate_fs(double *time, int n_data, int *index, int *f1_count) {

	int i, j, k, n;

	int error, count;
	int year1, year2, leap1, leap2;

	int n_attempt = 9;
	int n_sample  = 7;

	int *bad, *idt;

	double s, dp, min_dt, max_dt, scale, t, ds;
	double mean, median, n1, n2, fs;
	double x1, x2, x3;

	double *dt, *dt2, *tmp;

	double *max_ds, *break_point;

	double sample[] = {0.041666666, 1.0, 2.0, 7.0, 14.0, 30.4375, 365.25};
	double daily[]  = {365.24219, 0.0, 0.0, 0.0, 365.0, 365.25, 365.24219, 365.2564, 365.259635};

	dt     = (double *) calloc((size_t) (n_data-1),  sizeof(double));
	dt2    = (double *) calloc((size_t) (n_data-1),  sizeof(double));
	max_ds = (double *) calloc((size_t) (n_attempt), sizeof(double));
	bad    = (int *)    calloc((size_t) (n_data-1),  sizeof(int));
	idt    = (int *)    calloc((size_t) (n_data-1),  sizeof(int));

	error = 0;
	for (j = 0; j < n_data-1; j++) {

		dt[j] = time[j+1] - time[j];

		year1 = (int)floor(time[j]);
		leap1 = year1%4 == 0 && year1%100 != 0 || year1%400 == 0;
		n1 = 365.0 + (double)leap1;

		year2 = (int)floor(time[j+1]);
		leap2 = year2%4 == 0 && year2%100 != 0 || year2%400 == 0;
		n2 = 365.0 + (double)leap2;

		dt2[j] = 0.5 * ( (n1+n2) / n2 / n1);

		if (dt[j] <= 0.0) bad[error++] = j;

		if (j == 0) min_dt = dt[j];
		else min_dt = dt[j] < min_dt ? dt[j] : min_dt; 
	}
	if (error) {
		fprintf(stderr, " calculate_fs : file is not in ascending");
		fprintf(stderr, " time order or has duplicated points\n");
		for (j = 0; j < error; j++) {
		       	fprintf(stderr, "              : t[%d] - t[%d] = %f\n", bad[j]+1, bad[j], dt[bad[j]]);
		}
		exit(EXIT_FAILURE);
	}
	free(bad);

	/******************************************************************/
	/* First find out the number of decimal places in the time stamps */
	/******************************************************************/

	dp = 0.0;
	do {
		dp += 1.0;
		s = 0.0;
		scale = pow(10.0,dp);
		for (j = 0; j < n_data; j++) {
			t = rint(time[j] * scale) / scale;
			s += ((t-time[j])*(t-time[j]));
		}
		s = sqrt(s / (double) n_data);	
	} while (s > 1e-13);

	/************************************************************************/
	/* The 1e-13 was arrived at partly from trial and error and it is       */
	/* related to the numerical precision of the computer, meaning it       */
	/* can only really distinguish up to 12 decimal places                  */
	/* For a specified decimal place, the differences between the full      */
	/* value and the "clipped" value should form a uniform distribution     */
	/* between -0.5 / 10^dp and 0.5 / 10^dp. Therefore the square root of   */
	/* the sum of the differences squared should have a mean equal to       */
	/* (N / 12 / 10^(2*dp))^0.5. The difficultly is for small time series   */
	/* say less the 250 the error on the mean is quite large and could      */
	/* possibly overlap with the dp level below. A better method will be    */
	/* developed later, its not really critical at this stage, for well     */
	/* defined time series.                                                 */
	/************************************************************************/

	fprintf(stderr, "min_dt = %f dp = %f\n", min_dt, dp);

	/*****************************************************************************/
	/* Try to find out how evenly spaced the file is, and the sampling frequency */
	/*****************************************************************************/

	tmp = (double *) calloc((size_t)(n_data-1), sizeof(double));

	count = 0;
	for (j = 0; j < n_data-1; j++) {
		if (lrint( dt[j] / min_dt) == 1) {
			if (count == 0) max_dt = dt[j];
			else max_dt = dt[j] > max_dt ? dt[j] : max_dt;
			tmp[count] = dt[j];
		       	count++;
		}
	}

	qsort((char *)tmp, (size_t)count, sizeof(double), comp_double_asc);
	if ((int)fmod(count,2) == 1) {
		median = tmp[(count+1)/2 - 1];
	} else {
                median = ( tmp[count/2-1] + tmp[count/2] ) / 2.0;
	}

	mean = 0.0;
	for (j = 0; j < count; j++) mean += tmp[j];
	mean /= (double) count;

	free(tmp);

	fprintf(stderr, "median = %f mean = %f\n", median, mean);

	/***************************************/
	/* Work out what sampling it really is */
	/***************************************/

	break_point = (double *) calloc((size_t)n_sample+1, sizeof(double) );
	break_point[0] = 0;
	scale = pow(10.0,dp);
	for (j = 0; j < n_sample-1; j++) {
		x1 = rint(sample[j] / 365.25 * scale) / scale; 
		x2 = rint(sample[j+1] / 365.25 * scale) / scale;
		x3 = (x2 - x1) * scale / 2.0;
		break_point[j+1] = x1 + x3 / scale;
	}
	break_point[n_sample] = 10.0 * x2;

	for (j = 0; j < n_sample; j++) fprintf(stdout, "break_point %f\n", break_point[j]);

	i = -1;
	n = 0;
	for (j = 0; j < n_sample; j++) {
		if (median >= break_point[j] && median < break_point[j+1]) {
			fprintf(stderr, "sample = %f\n", sample[j]);
		       	i = j;
			n++;
		}
	}
	if (i < 0 || n > 1) {
		fprintf(stderr, " calculate_fs : numerical precision (%f decimal", dp);
		fprintf(stderr, " places) of time samples too small to calculate frequency\n");
		exit(EXIT_FAILURE);
	}

	for (j = 4; j < n_attempt; j++) daily[j]  /= sample[i];

	daily[1] = 1.0 / min_dt;
	daily[2] = 1.0 / median;
	daily[3] = 1.0 / mean;

	for (j = 0; j < n_attempt; j++) fprintf(stderr, "%f\n", daily[j]);

	for (j = 0; j < n_data-1; j++) {
		ds = rint( dt[j] / dt2[j]) - (dt[j] / dt2[j]);
		idt[j] = lrint(dt[j] / dt2[j]); 
		fprintf(stderr, "%f %f %f %d\n", dt[j], dt2[j], ds, idt[j]);
		if (j == 0) max_ds[0] = fabs(ds); 
		else max_ds[0] = fabs(ds) > max_ds[0] ? fabs(ds) : max_ds[0];
	}
	fprintf(stderr, "max_ds[0] = %f\n", max_ds[0]);

	for (k = 1; k < n_attempt; k++) {
		for (j = 0; j < n_data-1; j++) {
			ds = rint( dt[j] * daily[k]) - (dt[j] * daily[k]);
			if (j == 0) max_ds[k] = fabs(ds); 
			else max_ds[k] = fabs(ds) > max_ds[k] ? fabs(ds) : max_ds[k];
		}
		fprintf(stderr, "max_ds[%d] = %f\n", k, max_ds[k]);
	}

	for (j = 0; j < n_attempt; j++) fprintf(stderr, "%f\n", max_ds[j] / daily[j]);

	i = 0;
	for (k = 1; k < n_attempt; k++) i = max_ds[k] < max_ds[i] ? k : i;

	fprintf(stderr, "Choice %d %f\n", i, daily[i]);
	count = 0;
	index[0] = 0;
	for (j = 0; j < n_data-1; j++) {
	       	if (i > 0) idt[j] = lrint( dt[j] * daily[i]);
		if (idt[j] == 1) count++;
		index[j+1] = index[j] + idt[j];
	}

	fs = daily[i] / 365.24219 / 24.0 / 3600; 
	
	(*f1_count) = count;

	free(dt);
	free(dt2);
	free(max_ds);
	free(idt);
	free(break_point);

	return(fs);
}