emp_wrapper.c

最大似然估计算法

#include "timeseries.h"

void emp_wrapper(noise_model *nm, int n_models, time_series ts, data_kernel dk, options op) {

	int i, j, k, l, m;
	int error;

	double sig_wh, sig_cn, min_cn;

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


	if (n_models == 1 && nm[0].model == 'w') {
		white_only(ts, dk, nm[0], op, 1);
	} else if (n_models == 2) {
		if (nm[0].model == 'w' || nm[1].model == 'w') {
			k = (nm[0].model == 'w' ? 1 : 0);
			l = (nm[0].model == 'w' ? 0 : 1);

			cats_empirical(ts, dk, nm[k], op, &sig_wh, &sig_cn, &min_cn);

			nm[k].sigma[0] = sig_cn;
			nm[l].sigma[0] = sig_wh;

			if (op.speed == 0) {
				nm[l].sigma_flag[0] = 1;
				nm[k].sigma_flag[0] = 1;

				wls_wrapper(nm,n_models,ts,dk,op);

			} else if (op.speed == 1 || op.speed == 2) {
				white_only(ts, dk, nm[l], op, 1);
				dk.MLE[0] = 0.0;
				
				for (j = 0; j < dk.n_par; j++) {
					for (i = 0; i < dk.n_par; i++) {
						dk.covar[j + i * dk.n_par] = 0.0;
					}
				}
				nm[l].sigma[0] = sig_wh;
				nm[l].dsigma[0] = 0.0;
				nm[k].sigma_flag[0] = 2;
			} else {
				dk.MLE[0] = 0.0;
				for (j = 0; j < dk.n_par; j++) {
					dk.params[j] = 0.0;
					for (i = 0; i < dk.n_par; i++) {
						dk.covar[j + i * dk.n_par] = 0.0;
					}
				}
				nm[k].sigma_flag[0] = 2;
				nm[l].sigma_flag[0] = 2;
			}
		}
	} else {
		fprintf(stderr, " Error : empirical : desired stochastic model can not be estimated empirically, use MLE or PSD methods\n");
		exit(EXIT_FAILURE);
	}
}