color_white_mle.c

最大似然估计算法

#include "timeseries.h"

double color_white_mle(double *params, int nrow, double *Eig_values, double *Eig_vectors, double *residual){

int i, j, k;
int flag;

double a, a1, a2, b, sum, out, alpha, beta;
double *temp, *work1, *work2, *work3;

alpha = 1.0;
beta  = 0.0;

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

flag = 0;
for (j = 0; j < 2; j++) {
	if (params[j] < 0.0) { 
		flag = 1;
		params[j] = fabs(params[j]);
	}
}

/* a = (params[0] * params[0]) / (params[1] * params[1]); */
a1 = params[1] * params[1];
a2 = params[0] * params[0];

for (sum = 0.0, j = 0; j < nrow; j++) {
	temp[j] = a1 / (a2 + a1 * Eig_values[j]);
	sum += log(params[0]*params[0] + params[1]*params[1]*Eig_values[j]);
}

for (j = 0; j < nrow; j++) {
	for (k = 0; k < nrow; k++) {
		work1[j + k * nrow] = temp[j] * Eig_vectors[k + j * nrow];
	}
}

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

b = 0.0;
for (j = 0; j < nrow; j++) b += work3[j] * residual[j];

out = -0.5 * (double) nrow * log(2.0 * M_PI) - 0.5 * sum - 0.5 * b / (params[1] * params[1]);

if (flag) out -= 1e6;

free(temp);
free(work1);
free(work2);
free(work3);
	
return(out);
}