cats_cn_only.c

最大似然估计算法

#include "timeseries.h"

double cats_CN_only(data_kernel dk, double *C_unit, double *Eig_values, double *Eig_vectors, double *params_mle, double *cov_mle, int speed) {

	int i, j, k, info;
	int n_params, n_data, n_par, lwork;

	int *ipiv;

	double a, b, P, N, alpha, beta, min_mle;
	double F1, F2, F3, F4;
	double second_partial, cross_partial;

	double *Ioo, *work;
	double *fit, *cov_fit, *data_hat, *params;
	double *work1, *work2, *work3;
	double *residuals;

	/******************************************************/
	/* This is for one color noise only so n_params = 1   */
	/* with fixed parameters                              */
	/******************************************************/

	n_params = 1;
	n_data = dk.n_data;
	n_par  = dk.n_par;

	/**********************************************************************************/
	/* First of all generate all the relevant matrices for the rest of the subroutine */
	/**********************************************************************************/

	residuals   = (double *) calloc((size_t) n_data,          sizeof(double));
	fit         = (double *) calloc((size_t) n_par,           sizeof(double));
	cov_fit     = (double *) calloc((size_t)(n_par*n_par),    sizeof(double));
	data_hat    = (double *) calloc((size_t) n_data,          sizeof(double));

	params      = (double *) calloc((size_t)(n_par+n_params), sizeof(double));

	work1       = (double *) calloc( (size_t)(n_data*n_data), sizeof(double));
	work2       = (double *) calloc( (size_t) n_data,         sizeof(double));
	work3       = (double *) calloc( (size_t) n_data,         sizeof(double));

	/**********************************************/
	/* Now estimate the parameters and covariance */
	/**********************************************/

	if (speed < 3) {
		linefit_all_fast(dk.A,dk.d,C_unit,n_data,dk.n_par,fit,cov_fit,data_hat);
		for (j = 0; j < n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
	} else {
		for (j = 0; j < n_data; j++) residuals[j] = dk.d[j];
		for (j = 0; j < dk.n_par; j++) params[j] = 0.0;
	}

	/**************************************************/
	/* Calculate the amplitude of the power-law noise */
	/**************************************************/

	alpha = 1.0;
	beta  = 0.0;

	N = (double) n_data;

	for (j = 0; j < n_data; j++) {
		for (k = 0; k < n_data; k++) {
			work1[j + k * n_data] = Eig_vectors[k + j * n_data] / Eig_values[j];
		}
	}

	i = 1;
	dgemm_("T","N",&i,&n_data,&n_data,&alpha,residuals,&n_data,Eig_vectors,&n_data,&beta,work2,&i);
	dgemm_("N","N",&i,&n_data,&n_data,&alpha,work2,&i,work1,&n_data,&beta,work3,&i);

	P = 0.0;
	for (j = 0; j < n_data; j++) P += work3[j] * residuals[j];

	b = sqrt( P / N);

	/**********************/
	/* Free some matrices */
	/**********************/

	free(work1);
	free(work2);
	free(work3);

	/************************************/
	/* Store the params into params_mle */
	/************************************/

	i = n_par + n_params;

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

	for (j = 0; j < n_par; j++) {
		for (k = 0; k < n_par; k++) cov_fit[j + k * n_par] *= (b * b);
	}

	/************************************************************/
	/* Now form the covariance matrix                           */
	/* We dont need to do the linear parts as they are already  */
	/* calculated in cov_fit. The amplitude of the CN noise     */
	/* the is no correlation with the linear parts so there is  */
	/* no need to calculate the cross partials                  */
	/************************************************************/

	/*********************************************/
	/* Invert cov_fit and substitute it into Ioo */
	/*********************************************/

	Ioo  = (double *) calloc((size_t)(i*i),  sizeof(double));

	lwork = n_par;
	ipiv = (int *)    calloc((size_t) n_par, sizeof(int));
	work = (double *) calloc((size_t) lwork, sizeof(double));

	dgetrf_(&n_par, &n_par, cov_fit, &n_par, ipiv, &info);
	dgetri_(&n_par, cov_fit, &n_par, ipiv, work, &lwork, &info);

	free(ipiv);
	free(work);

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

	/***************************************************************/
	/* Now calculate the second_partial for the CN noise parameter */
	/***************************************************************/

	F2 = MLE_nparam_CN_only(params_mle, i, n_par, n_data, dk.d, dk.A, Eig_values, Eig_vectors);

	work = (double *) calloc((size_t) (i), sizeof(double));
	for (j = 0; j < i; j++) work[j] = params_mle[j];

	work[n_par] = params_mle[n_par] + 1e-3 * fabs(params_mle[n_par]);
	F1 = MLE_nparam_CN_only(work, i, n_par, n_data, dk.d, dk.A, Eig_values, Eig_vectors);

	work[n_par] = params_mle[n_par] - 1e-3 * fabs(params_mle[n_par]);
	F3 = MLE_nparam_CN_only(work, i, n_par, n_data, dk.d, dk.A, Eig_values, Eig_vectors);

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

	free(work);

	/**********************/
        /* Calculate inv(Ioo) */
        /**********************/

        lwork = i;
        ipiv = (int    *) calloc((size_t) i,     sizeof(int));
        work = (double *) calloc((size_t) lwork, sizeof(double));

        dgetrf_(&i, &i, Ioo, &i, ipiv, &info);
        dgetri_(&i, Ioo, &i, ipiv, work, &lwork, &info);
        free(ipiv);
        free(work);

	/***************************/
        /* Store -Ioo into cov_mle */
        /***************************/

        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(residuals);
	free(fit);
	free(cov_fit);
	free(data_hat);

	return(F2);
}

/*******************************************************************************************/
/* Stripped down version of color_only_mle that doesnt calculate the covariance at the end */
/*******************************************************************************************/

double cats_CN_only_stripped(data_kernel dk, double *C_unit, double *Eig_values, double *Eig_vectors, double *params_mle, int speed) {

	int i, j, k, info;
	int n_params, lwork, n_data, n_par;

	double a, b, P, N, alpha, beta, min_mle;
	double fv, F2;

	double *fit, *cov_fit, *data_hat, *params;
	double *work1, *work2, *work3;
	double *residuals;


	/*****************************************************/
	/* This is for a coloured noise only so n_params = 1 */
	/*****************************************************/

	n_params = 1;
	n_data = dk.n_data;
	n_par = dk.n_par;

	/**********************************************************************************/
	/* First of all generate all the relevant matrices for the rest of the subroutine */
	/**********************************************************************************/

	residuals = (double *) calloc((size_t)  n_data,           sizeof(double));
	fit       = (double *) calloc((size_t)  n_par,            sizeof(double));
	cov_fit   = (double *) calloc((size_t) (n_par*n_par),     sizeof(double));
	data_hat  = (double *) calloc((size_t)  n_data,           sizeof(double));

	params    = (double *) calloc((size_t) (n_par+n_params),  sizeof(double));

	work1     = (double *) calloc((size_t) (n_data * n_data), sizeof(double) );
	work2     = (double *) calloc((size_t)  n_data,           sizeof(double) );
	work3     = (double *) calloc((size_t)  n_data,           sizeof(double) );

	/**********************************************/
	/* Now estimate the parameters and covariance */
	/**********************************************/

	if (speed < 3) {
		linefit_all_fast(dk.A,dk.d,C_unit,n_data,dk.n_par,fit,cov_fit,data_hat);
		for (j = 0; j < n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
	} else {
		for (j = 0; j < n_data; j++) residuals[j] = dk.d[j];
		for (j = 0; j < dk.n_par; j++) params[j] = 0.0;
	}

	/*************************************************/
	/* Calculate the amplitude of the coloured noise */
	/*************************************************/

	alpha = 1.0;
	beta  = 0.0;

	N = (double) n_data;

	for (j = 0; j < n_data; j++) {
		for (k = 0; k < n_data; k++) {
			work1[j + k * n_data] = Eig_vectors[k + j * n_data] / Eig_values[j];
		}
	}

	i = 1;
	dgemm_("T","N",&i,&n_data,&n_data,&alpha,residuals,&n_data,Eig_vectors,&n_data,&beta,work2,&i);
	dgemm_("N","N",&i,&n_data,&n_data,&alpha,work2,&i,work1,&n_data,&beta,work3,&i);

	P = 0.0;
	for (j = 0; j < n_data; j++) P += work3[j] * residuals[j];

	b = sqrt( P / N);

	/**********************/
	/* Free some matrices */
	/**********************/

	free(work1);
	free(work2);
	free(work3);

	/************************************/
	/* Store the params into params_mle */
	/************************************/

	i = n_par + n_params;

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

	F2 = MLE_nparam_CN_only(params_mle, i, n_par, n_data, dk.d, dk.A, Eig_values, Eig_vectors);
	fv = -F2;

	free(params);
	free(residuals);
	free(fit);
	free(cov_fit);
	free(data_hat);

	return(fv);
}

double cats_CN_only_stripped2(time_series ts, data_kernel dk, noise_model nm, options op, double *params_mle) {

        int j, k, i, lwork, liwork;
        int info, c_id;
        
        int *iwork;
        
        double cn_mle;
        
        double *Eig_vectors, *C_unit, *Eig_values;
        double *work;
        
        /*******************************************************************************/
        /* First of all generate all the relevant matrices for the rest of the simplex */
        /*******************************************************************************/
        
        Eig_vectors = (double *) calloc((size_t)(dk.n_data*dk.n_data), sizeof(double));
        Eig_values  = (double *) calloc((size_t) dk.n_data,            sizeof(double));
        C_unit      = (double *) calloc((size_t)(dk.n_data*dk.n_data), sizeof(double));
        
        /* lwork = 100 + 6 * dk.n_data + 2 * dk.n_data * dk.n_data; 
        liwork = 30 + 5 * dk.n_data; */

	lwork = 35 * dk.n_data;
	liwork = 5 * dk.n_data;

        
        work        = (double *) calloc((size_t) lwork,                sizeof(double));
        iwork       = (int *)    calloc((size_t) liwork,               sizeof(int));
        
        /*************************************************************************/
        /* Create the unit power law covariance matrix and find the eigen values */
        /* and the eigen vectors                                                 */
        /*************************************************************************/

        create_covariance(nm, ts, op.cov_method, Eig_vectors, 1);

        dlacpy_("Full", &dk.n_data, &dk.n_data, Eig_vectors, &dk.n_data, C_unit, &dk.n_data);
        /* dsyevd_("V","L",&dk.n_data,Eig_vectors,&dk.n_data,Eig_values,work,&lwork,iwork,&liwork,&info); */
	dsyev_("V","L",&dk.n_data,Eig_vectors,&dk.n_data,Eig_values,work,&lwork,&info);


        free(iwork);
        free(work);

        i = dk.n_par + 1;

        cn_mle = cats_CN_only_stripped(dk, C_unit, Eig_values, Eig_vectors, params_mle, op.speed);

	free(Eig_vectors);
	free(Eig_values);
	free(C_unit);

	return(cn_mle);
}