mle_nparam.c

最大似然估计算法

#include "timeseries.h"

double MLE_nparam_CN(double *all_params, int n_params, int n_lpar, int n_data, double *d, double *A, double *Eig_values, double *Eig_vectors){

	int i, j;

	double out, alpha, beta;
	double *residual, *params, *fit;

	alpha = 1.0;
	beta  = 0.0;

	i = n_params - n_lpar;
	residual = (double *) calloc( (size_t) n_data,            sizeof(double) );
	fit      = (double *) calloc( (size_t) n_lpar,            sizeof(double) );
	params   = (double *) calloc( (size_t) i,                 sizeof(double) );

	/***********************/
	/* Get residuals first */
	/***********************/

	for (j = 0; j < n_lpar; j++) fit[j]    = all_params[j];
	for (j = 0; j < i; j++)      params[j] = all_params[j+n_lpar];

	i = 1;
	dgemm_("N","N",&n_data,&i,&n_lpar,&alpha,A,&n_data,fit,&n_lpar,&beta,residual,&n_data);
	for (j = 0; j < n_data; j++) residual[j] = d[j] - residual[j];

	out = MLE_withline_CN(params, n_data, Eig_values, Eig_vectors, residual); 

	free(residual);
	free(fit);
	free(params);
	
	return(out);
}

double MLE_nparam_CN_only(double *all_params, int n_params, int n_lpar, int n_data, double *d, double *A, double *Eig_values, double *Eig_vectors){

	int i, j;

	double out, alpha, beta, param;
	double *residual, *fit;

	alpha = 1.0;
	beta  = 0.0;

	i = n_params - n_lpar;
	residual = (double *) calloc( (size_t) n_data,            sizeof(double) );
	fit      = (double *) calloc( (size_t) n_lpar,            sizeof(double) );

	/***********************/
	/* Get residuals first */
	/***********************/

	for (j = 0; j < n_lpar; j++) fit[j] = all_params[j];
	param = all_params[n_lpar];

	i = 1;
	dgemm_("N","N",&n_data,&i,&n_lpar,&alpha,A,&n_data,fit,&n_lpar,&beta,residual,&n_data);
	for (j = 0; j < n_data; j++) residual[j] = d[j] - residual[j];

	out = MLE_withline_CN_only(param, n_data, Eig_values, Eig_vectors, residual); 

	free(residual);
	free(fit);
	
	return(out);

}

double MLE_nparam_WH(double *all_params, int n_lpar, int n_data, double *d, double *A){

	int i, j;
	
	double alpha, beta, param, out;
	double *fit, *residual;

	alpha = 1.0;
	beta  = 0.0;

	/***********************/
	/* Get residuals first */
	/***********************/

	fit      = (double *) calloc((size_t) n_lpar, sizeof(double));
	residual = (double *) calloc((size_t) n_data, sizeof(double));

	for (j = 0; j < n_lpar; j++) fit[j] = all_params[j];
	param = all_params[n_lpar];

	i = 1;
	dgemm_("N","N",&n_data,&i,&n_lpar,&alpha,A,&n_data,fit,&n_lpar,&beta,residual,&n_data);
	for (j = 0; j < n_data; j++) residual[j] = d[j] - residual[j];

	out = MLE_withline_WH(param, n_data, residual);

	free(fit);
	free(residual);

	return(out);

}