est_line_wh_only.c

最大似然估计算法

#include "timeseries.h"

double est_line_WH_only(data_kernel dk, options op, double *params_mle, double *cov_mle) {

	int j, k, i, lwork, info;

	int *ipiv;

	double sum, fv, F1, F2, F3, F4, std;
	double second_partial, cross_partial;

	double *dy, *fit, *cov_fit, *data_hat, *residual, *params, *work;
	double *Ioo;

	dy       = (double *) calloc((size_t)  dk.n_data,          sizeof(double));
	fit      = (double *) calloc((size_t)  dk.n_par,           sizeof(double));
	cov_fit  = (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)); 
	params   = (double *) calloc((size_t) (dk.n_par+1),        sizeof(double));
	work     = (double *) calloc((size_t) (dk.n_par+1),        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_fit, 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]);
	}

	/*****************************************************************/
	/* I have a query here, now in general it is believed that       */
	/* std should be sqrt( sum / (n-1) ) and rms is                  */
	/* sqrt( sum / n ). However if you take the MLE equation,        */
	/* for just power law noise (not power law plus white noise)     */
	/* then you can solve for the noise amplitude                    */
	/* However the equation comes out as, if you set the spectral    */
	/* index to zero (white noise), sqrt( sum / n). In Hadleys       */
	/* matlab routines and here we have set the amplitude equal      */
	/* to the std and so we get a subtle but important difference    */
	/* in MLE and noise amplitude to that if                         */
	/* if we arrived at white noise only from the power plus white   */
	/* model which is important when we are testing different        */
	/* "models". So                                                  */
	/* for that reason from now on I am setting (4/10/2002) std      */
	/* equal to sqrt(sum / n). As n becomes large this becomes a     */
	/* mute point but still worth bearing in mind                    */
	/*****************************************************************/

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

	for (j = 0; j < dk.n_par; j++) params_mle[j] = fit[j];

	F2 = MLE_withline_WH(std, dk.n_data, residual);
	fv = -F2;

	/****************************************************************************/
	/* This used to use the Fisher Information Matrix to calculate the full     */
	/* covariance matrix. However since the covariance matrix for the linear    */
	/* parameters is already calculated and the white noise amplitude is always */
	/* uncorrelated with the other parameters then we need only estimate the    */
	/* curvature for the white noise amplitude                                  */
	/****************************************************************************/

	i = dk.n_par + 1;
	for (j = 0; j < i; j++) work[j] = params[j] = params_mle[j];
	Ioo = (double *) calloc((size_t)(i*i), sizeof(double) );

	for (j = 0; j < dk.n_par; j++) {
		for (k = 0; k < dk.n_par; k++) {
			Ioo[j + k * i] = cov_fit[j + k * dk.n_par] * std * std; 
		}
	}

	params[i-1] = work[i-1] + 1e-3 * fabs(work[i-1]);
	F1 = MLE_withline_WH(params[i-1],dk.n_data, residual);
	params[i-1] = work[i-1] - 1e-3 * fabs(work[i-1]);
	F3 = MLE_withline_WH(params[i-1],dk.n_data, residual);

	second_partial = -(F1-2.0*F2+F3)/pow(1e-3*fabs(work[i-1]),2.0);
	Ioo[ i * i - 1] = -1 / second_partial;

	for (j = 0; j < i; j++) {
        	for (k = 0; k < i; k++) cov_mle[j + k * i] = Ioo[j + k * i];
	}

	free(Ioo);
	free(params);
	free(dy);
	free(fit);
	free(cov_fit);
	free(data_hat);
	free(residual);
	free(work);

	return(fv);
}