linefit_fast.c

最大似然估计算法

#include "timeseries.h"

void linefit_fast(double *A, double *d, double *dy, int nrow, int ncol, double *fit, double *cov_fit, double *data_hat) {

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

double alpha, beta;
double *Q, *R, *t1, *new_A, *new_d, *work;

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++) {
	new_d[j] = d[j] / dy[j];
	for (k = 0; k < ncol; k++) new_A[j + k * nrow] = A[j + k * nrow] / dy[j];
}

/*********************/
/* 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) */
/********************/

ipiv = (int *) calloc((size_t) ncol, sizeof(int));
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,data_hat,&nrow);

free(new_d);
free(R);
free(t1);

}