mle_withline.c

最大似然估计算法

#include "timeseries.h"

double MLE_withline_CN(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]);

out = -out;
if (flag) out += 1e6;

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

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

	int j, k, i;

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

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

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

	out = -out;

	return(out);
}

double MLE_withline_WH(double param, int n_data, double *residual) {

int j;

double b, out;

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

out = -0.5 * (double) n_data * 2.0 * log(param) - 0.5 * b / (param * param) - 
	(double) n_data / 2.0 * log(2.0 * M_PI);

return(out);

}