power_law_cov.c

最大似然估计算法

#include "timeseries.h"

void power_law_cov(double *time, int n, double Delta_t, double sigma, double spectral_index, double *C) {

	/*****************************************************************/
	/*    power_law_cov                                              */
	/*                                                               */
	/*    Compute the covariance matrix for any power-law process    */
	/*    using fractional differencing (see fractional_diff.c)      */
	/*                                                               */
	/*    time           time vector of data points (years)          */
	/*    sigma          standard deviation of noise (meters)        */
	/*    spectral_index                                             */
	/*                                                               */
	/* Hadley Johnson, [e-mail], 28-May-1998 matlab version */
	/* Simon  Williams, [e-mail], 08-Nov-2000 c-version       */
	/*                                                               */
	/* added Delta_t to calculate first offset instead of assuming   */
	/* 1 / 365.24219 so it is better for data other than daily       */
	/* for example monthly data in sea-level records. Could use the  */
	/* median delta instead but that requires sorting the data which */
	/* I didnt want to do. Also could have specified/calculated      */
	/* Delta_t in the calling subroutines but that means changing    */
	/* a lot of routines. May do this later when I have more time    */
	/* S. Williams 18 Dec 2002                                       */
	/*****************************************************************/

	int i, j, k;

	double *dt, *c, *filter;

	dt = (double *) calloc((size_t)n, sizeof(double) );
	c  = (double *) calloc((size_t)n, sizeof(double) );

	filter = (double *)calloc((size_t)(n*n), sizeof(double));

	dt[0] = Delta_t;
	for (j = 0; j < n-1; j++) { 
		dt[j+1] = time[j+1] - time[j];
	}

	for (j = 0; j < n; j++) c[j] = sigma * pow(dt[j],-spectral_index / 4.0);

	if (spectral_index > 0.0) c[0] = 0.0;

	fractional_diff(filter,n,spectral_index / 2.0);

	for (j = 0; j < n; j++) {
		for (k = 0; k < n; k++){
			filter[j + k * n] *= c[k];
		}
	}

	for (j = 0; j < n; j++) {
		for (k = 0; k < n; k++) {
			C[j + k * n] = 0.0;
			for (i = 0; i < n; i++) {
				C[j + k * n] += filter[j + i * n] * filter[k + i * n];
			}
		}
	}
	free(dt);
	free(c);
	free(filter);

}