mle_fixed_param.c

来自「最大似然估计算法」· C语言 代码 · 共 68 行

#include "timeseries.h"

double MLE_fixed_param(double **Cunits, double *params, int n_models, data_kernel dk) {

int i, j, k, l, info, nrhs;
int nrow, nsquared, ncol;

int *ipiv;

double alpha, beta, P, N, out, s2, a;
double *Cfull, *work, *fit;

nrow = dk.n_data;
ncol = dk.n_par;
nsquared = nrow * nrow;

fit   = (double *) calloc((size_t) ncol,     sizeof(double));
Cfull = (double *) calloc((size_t) nsquared, sizeof(double));

for (j = 0; j < ncol; j++) fit[j] = params[j];


for (j = 0; j < n_models; j++) {
	s2 = params[ncol+j] * params[ncol+j];
	for (l = 0; l < nsquared; l++) Cfull[l] += s2*Cunits[j][l];
}

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

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

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

alpha = 1.0;
beta  = 0.0;

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


dgemm_("N","N",&nrow,&i,&ncol,&alpha,dk.A,&nrow,fit,&ncol,&beta,work,&nrow);

nrhs = 1;
for (j = 0; j < nrow; j++) work[j] = dk.d[j] - work[j];

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

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

free(ipiv);
free(work);
free(fit);
free(Cfull);

N = (double) nrow;

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

return(out);
}