color_only_mle.c

最大似然估计算法

#include "timeseries.h"

double color_only_mle(double *sigma, int n_data, double *Eig_values, double *Eig_vectors, double *residuals, int get_sigma) {

	int j, k, i;

	double alpha, beta, N, a, P, out, sig;

	double *work1, *work2, *work3;

	alpha = 1.0;
	beta  = 0.0;

	N = (double) n_data;

	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) );

	for (a = 0.0, j = 0; j < n_data; j++) a += log(Eig_values[j]);

	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);

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

	free(work1);
	free(work2);
	free(work3);
	
	if (get_sigma) {
		sig = sqrt(P / N);
		*sigma = sig;	
	}
	else sig = *sigma;

	out = -N * log(2.0 * M_PI) / 2.0 - a / 2.0 - N * log(sig*sig) / 2.0 - P / 2.0 / sig / sig;

	return(out);
}