wls_wrapper.c

最大似然估计算法

#include "timeseries.h"
                                                                                                                                                  
void wls_wrapper(noise_model *nm, int n_models, time_series ts, data_kernel dk, options op) {

int i, j, k;
int n_squared, error;

double **Cunits, *Cfull;
double *fit, *cov_fit, *r, scale, mle;

/********************************************************/
/* For this every parameter must be fixed so check this */
/* and exit if it is not true                           */
/********************************************************/

for (j = 0, error = 0; j < n_models; j++) {
	if (nm[j].sigma_flag[0] != 1) error++;
	for (k = 0; k < nm[j].n_pvec; k++) {
		if (nm[j].pvec_flag[k] != 1) error++;
	}
}
if (error) {
	fprintf(stderr, " Error : wls_wrapper : Not all parameters are fixed. Exiting...\n");
	exit(EXIT_FAILURE);
}

n_squared  = dk.n_data * dk.n_data;
                                                                                                                                                  
Cunits   = (double **) calloc((size_t) n_models*n_squared, sizeof(double *));
for (j = 0; j < n_models; j++) Cunits[j] = (double *) calloc((size_t) n_squared, sizeof(double));
                                                                                                                                                  
Cfull    = (double *)  calloc((size_t) n_squared, sizeof(double));
                                                                                                                                                  
fit      = (double *) calloc((size_t) dk.n_par,            sizeof(double));
cov_fit  = (double *) calloc((size_t) (dk.n_par*dk.n_par), sizeof(double));
r        = (double *) calloc((size_t) 1,                   sizeof(double));

for (j = 0; j < n_models; j++) {
	create_covariance(nm[j], ts, op.cov_method, Cunits[j], 1);
	scale = nm[j].sigma[0] * nm[j].sigma[0];
	for (k = 0; k < n_squared; k++) Cfull[k] += scale*Cunits[j][k];
}

mle = general_MLE(dk, Cfull, op, fit, cov_fit, r, 0);

if (op.verbose) {
	fprintf(op.fpout, " Weighted Least Squares : MLE would be maximised (at the current sigma ratios) by scaling by %f\n", r[0]);
}

dk.MLE[0] = mle;

for (j = 0; j < dk.n_par; j++) {
	dk.params[j] = fit[j];
	for (k = 0; k < dk.n_par; k++) {
		dk.covar[j + k * dk.n_par] = cov_fit[j + k * dk.n_par];
	}
}

free(fit);
free(cov_fit);
free(r);
free(Cfull);
for (j = 0; j < n_models; j++) free(Cunits[j]);
free(Cunits);

}