linefit_all_fast2.c

最大似然估计算法

#include "timeseries.h"

double linefit_all_fast2(double *A, double *data, double *C, int nrow, int ncol, double *fit, double *cov_fit, double *residuals, double *sigma, int get_sigma) {

int i, j, k, info, nrhs, lwork, minrc;
int *ipiv;

double alpha, beta, a, P, sig, N, out;
double *W, *W2, *Q, *R, *t1, *new_A, *new_d, *work;

clock_t time1, time2;

lwork = (nrow * (ncol+1));

W     = (double *) calloc((size_t) (nrow*nrow), sizeof(double));
W2    = (double *) calloc((size_t) (nrow*nrow), sizeof(double));
work  = (double *) calloc((size_t) lwork,       sizeof(double));

ipiv = (int *) calloc((size_t) nrow, sizeof(int));

dlacpy_("Full",&nrow,&nrow,C,&nrow,W,&nrow); 

/***************/
/* C = chol(C) */
/***************/
 
time1 = clock();

dpotrf_("L",&nrow,W,&nrow,&info);

for (a = 0.0, j = 0; j < nrow; j++) {
	a += log(W[j + j * nrow]);
	for (k = j+1; k < nrow; k++) W[j + k * nrow] = 0.0;
}
a = a * 2.0;

dlacpy_("Full",&nrow,&nrow,W,&nrow,W2,&nrow); 

time2 = clock();

fprintf(stderr, " Time taken to create cholesky of covariance matrix : ");
fprintf(stderr, "%f seconds\n", ((double) (time2-time1)) / CLOCKS_PER_SEC);

nrhs = ncol + 1;
for (j = 0; j < nrow; j++) {
	for (k = 0; k < ncol; k++) work[j + k * nrow] = A[j + k * nrow];
	work[j + ncol * nrow] = data[j];
}


/**************************/
/* New_d = C \ data       */
/* New_A = C \ A          */
/* both computed together */
/**************************/

dgesv_(&nrow, &nrhs, W, &nrow, ipiv, work, &nrow, &info);

new_A = (double *) calloc((size_t) (nrow*ncol), sizeof(double));
new_d = (double *) calloc((size_t) nrow,        sizeof(double));

for (j = 0; j < nrow; j++) {
	for (k = 0; k < ncol; k++) new_A[j + k * nrow] = work[j + k * nrow];
	new_d[j] = work[j + ncol * nrow];
}

/*********************/
/* QR factorization  */
/* [Q,R] = qr(new_A) */
/*********************/

minrc = nrow < ncol ? nrow : ncol;

Q = (double *) calloc((size_t) (nrow*minrc), sizeof(double));
R = (double *) calloc((size_t) (minrc*ncol), sizeof(double));

qr(new_A, nrow, ncol, Q, R);
free(new_A);

/********************/
/* Calculate inv(R) */
/********************/

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

dgetrf_(&ncol, &ncol, R, &ncol, ipiv, &info);
dgetri_(&ncol, R, &ncol, ipiv, work, &lwork, &info);
free(ipiv);
free(work);

/********************/
/* t1 = inv(R) * Q' */
/********************/

alpha = 1.0;
beta  = 0.0;
t1 = (double *) calloc( (size_t) (ncol*nrow), sizeof(double));
dgemm_("N","T",&ncol,&nrow,&ncol,&alpha,R,&ncol,Q,&nrow,&beta,t1,&ncol);
free(Q);

/**********************/
/* cov_fit = t1 * t1' */
/**********************/

dgemm_("N","T",&ncol, &ncol,&nrow,&alpha,t1,&ncol,t1,&ncol,&beta,cov_fit,&ncol);

/********************/
/* fit = t1 * new_d */
/********************/

i = 1;
dgemm_("N","N",&ncol,&i,&nrow,&alpha,t1,&ncol,new_d,&nrow,&beta,fit,&ncol);

/**********************/
/* data_hat = A * fit */
/**********************/

dgemm_("N","N",&nrow,&i,&ncol,&alpha,A,&nrow,fit,&ncol,&beta,residuals,&nrow);

/********/
/* Done */
/********/

time1 = clock();

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

nrhs = 1;
for (j = 0; j < nrow; j++) {
	work[j]      = data[j] - residuals[j];
	residuals[j] = work[j];
}

dgesv_(&nrow, &nrhs, W2, &nrow, ipiv, work, &nrow, &info);

for (P = 0.0, j = 0; j < nrow; j++) P += work[j]*work[j];

free(ipiv);
free(work);

N = (double) nrow;

if (get_sigma) {
	sig = sqrt(P / N);
	*sigma = sig;
}
else sig = *sigma;
                                                                                                                                           
fprintf(stderr, " linefit_all_fast2 : %13.6f %13.6f %13.6f %13.6f\n", -N * log(2.0 * M_PI) / 2.0, - a / 2.0, - N * log(sig*sig) / 2.0, - P / 2.0 / sig / sig);

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

fprintf(stderr, " linefit_all_fast2 : mle : %f\n", out);

time2 = clock();

fprintf(stderr, " Time taken to do final MLE part : ");
fprintf(stderr, "%f seconds\n", ((double) (time2-time1)) / CLOCKS_PER_SEC);

free(new_d);
free(R);
free(t1);
free(W);
free(W2);

return(out);

}