linefit_all_fast3.c

最大似然估计算法

#include "timeseries.h"

double linefit_all_fast3(data_kernel dk, double *C, int nrow, double *fit, double *cov_fit, double *data_hat, options op) {

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

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

clock_t time1, time2;

ncol = dk.n_par;

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

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

time2 = clock();

if (op.verbose > 0) {
	fprintf(op.fpout, " Time taken to create cholesky of covariance matrix : ");
	fprintf(op.fpout, "%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] = dk.A[j + k * nrow];
	work[j + ncol * nrow] = dk.d[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,dk.A,&nrow,fit,&ncol,&beta,data_hat,&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] = dk.d[j]-data_hat[j];

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

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

sig = sqrt( P / (double) nrow);

free(ipiv);
free(work);

N = (double) nrow;

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

time2 = clock();

if (op.verbose > 1) {
	fprintf(op.fpout, " linefit_all_fast : %13.6f %13.6f %13.6f %13.6f\n", -N * log(2.0 * M_PI) / 2.0, - a, - N * log(sig*sig) / 2.0, - P / 2.0 / sig / sig);
	fprintf(op.fpout, " linefit_all_fast : mle : %f\n", out);


	fprintf(op.fpout, " Time taken to do final mle part : ");
	fprintf(op.fpout, "%f seconds\n", ((double) (time2-time1)) / CLOCKS_PER_SEC);
}

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

return(out);

}