calculate_fisher.c

最大似然估计算法

#include "timeseries.h"

void calculate_fisher(double **Cunits, data_kernel dk, noise_model *nm, int n_models, double *fit, double *cov_fit) {

int i, j, k, l, n_par, n_squared;
int lwork, info;

int *ipiv, *iwork;

double F1, F2, F3, F4, second_partial, cross_partial;

double *Ioo, *sigma, *work, *params;

n_squared = dk.n_data * dk.n_data;

sigma = (double *) calloc((size_t) n_models,  sizeof(double));

for (i = 0; i < n_models; i++) sigma[i] = nm[i].sigma[0];

/*********************************************/
/* Invert cov_fit and substitute it into Ioo */
/*********************************************/

n_par = dk.n_par;

i = n_par + n_models;

Ioo = (double *) calloc((size_t)(i*i), sizeof(double) );

lwork = n_par;
ipiv = (int *)    calloc((size_t) n_par, sizeof(int));
work = (double *) calloc((size_t) n_par, sizeof(double));

dgetrf_(&n_par, &n_par, cov_fit, &n_par, ipiv, &info);
dgetri_(&n_par, cov_fit, &n_par, ipiv, work, &lwork, &info);

free(ipiv);
free(work);

/*******************************************************************/
/* Now calculate an estimate of the covariance matrix of the MLE   */
/* based on the Fisher Information Matrix for the non-linear parts */
/*******************************************************************/

for (j = 0; j < n_par; j++) {
	for (k = 0; k < n_par; k++) Ioo[j+k*i] = -cov_fit[j+k*n_par];
}

/*************************************/
/* Now for the non-linear parameters */
/*************************************/


params = (double *) calloc( (size_t) i, sizeof(double) );
work   = (double *) calloc( (size_t) i, sizeof(double) );

for (j = 0; j < n_par; j++)     work[j]       = params[j]       = fit[j];
for (j = 0; j < n_models; j++)  work[j+n_par] = params[j+n_par] = sigma[j];

F2 = MLE_fixed_param(Cunits,params,n_models,dk);

for (j = n_par; j < i; j++) {

	params[j] = work[j] + 1e-3 * fabs(work[j]);
	F1 = MLE_fixed_param(Cunits,params,n_models,dk);

	params[j] = work[j] - 1e-3 * fabs(work[j]);
	F3 = MLE_fixed_param(Cunits,params,n_models,dk);
	
	second_partial = (F1-2.0*F2+F3)/pow(1e-3*fabs(work[j]),2.0);

	/* fprintf(stderr, "F1 = %f F2 = %f F3 = %f sp = %f\n", F1, F2, F3, second_partial); */

	params[j] = work[j];
	
	Ioo[j + j * i] = second_partial;

}

for (j = n_par; j < i; j++) {
	for (k = 0; k < j; k++) {
		params[j] = work[j] + 1e-3 / 2.0 * fabs(work[j]);
		params[k] = work[k] + 1e-3 / 2.0 * fabs(work[k]);
		F1 = MLE_fixed_param(Cunits,params,n_models,dk);

		params[j] = work[j] - 1e-3 / 2.0 * fabs(work[j]);
		params[k] = work[k] + 1e-3 / 2.0 * fabs(work[k]);
		F2 = MLE_fixed_param(Cunits,params,n_models,dk);

		params[j] = work[j] + 1e-3 / 2.0 * fabs(work[j]);
		params[k] = work[k] - 1e-3 / 2.0 * fabs(work[k]);
		F3 = MLE_fixed_param(Cunits,params,n_models,dk);

		params[j] = work[j] - 1e-3 / 2.0 * fabs(work[j]);
		params[k] = work[k] - 1e-3 / 2.0 * fabs(work[k]);
		F4 = MLE_fixed_param(Cunits,params,n_models,dk);

		cross_partial = (F1-F2-F3+F4) / (1e-3 * fabs(work[j]) * 1e-3 * fabs(work[k]));

	/*	fprintf(stderr, "F1 = %f F2 = %f F3 = %f F4 = %f cp = %f\n", F1, F2, F3, F4, cross_partial); */
		Ioo[j + k * i] = Ioo[k + j * i] = cross_partial;

		params[j] = work[j];
		params[k] = work[k];
	}
}

/* for (j = 0; j < i; j++) {
	for (k = 0; k < i; k++) {
		fprintf(stderr, "%9.4f ", Ioo[j + k * i]);
	}
	fprintf(stderr, "\n");
}
fprintf(stderr, "\n"); */

free(work);
free(params);

/**********************/
/* Calculate inv(Ioo) */
/**********************/

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

dgetrf_(&i, &i, Ioo, &i, ipiv, &info);
dgetri_(&i, Ioo, &i, ipiv, work, &lwork, &info);

free(ipiv);
free(work);

/* for (j = 0; j < i; j++) {
	for (k = 0; k < i; k++) {
		fprintf(stderr, "%9.4f ", -Ioo[j + k * i]);
	}
	fprintf(stderr, "\n");
}
fprintf(stderr, "\n"); */

for (j = 0; j < n_par; j++) {
	for (k = 0; k < n_par; k++) {
		cov_fit[j + k * n_par] = -Ioo[j + k * i];
	}
}


for (j = 0; j < n_models; j++) {
	k = j + n_par;
	nm[j].dsigma[0] = sqrt(-Ioo[k + k * i]); 
	nm[j].sigma_flag[0] = 2;
}

free(sigma);
free(Ioo);

}