general_inner_simplex.c

最大似然估计算法

#include "timeseries.h" 

double general_inner_simplex(time_series ts, data_kernel dk, noise_model *nm, int n_models, options op, int end_flag) {
		
	int    i, j, k, l, m;
	int    n_vert, n_params, ilo, ihi, inhi, info;
	int    nfun_eval, got_answer, n_squared;

	double **p, *mle, *psum, *params, *ptry, *start;
	double **Cunits, *Cfull, *unit, *fit, *cov_fit;

	double mle_try, mle_save, min_mle;
	double fac, fac1, fac2, r;
	double md1, md2, md3, dist, scale;

	clock_t time1, time2;


	/*************************************************/
	/* Just check that no models have fixed sigmas   */ 
	/* Unfortunately you cant fix a sigma when using */
        /* the angle method                              */
	/*************************************************/

	for (j = 0; j < n_models; j++) {
		if (nm[j].sigma_flag[0] != 0) {
			fprintf(stderr, "Cannot have fixed covariance in the stochastic model\n");
			fprintf(stderr, "exiting....\n");
			exit(EXIT_FAILURE);
		}
	}
	n_vert   = n_models;
	n_params = n_vert-1;

	/***********************************/
	/* Set the initial starting values */
	/***********************************/

	ptry   = (double *)  calloc((size_t) n_params, sizeof(double));
	psum   = (double *)  calloc((size_t) n_params, sizeof(double));
	start  = (double *)  calloc((size_t) n_params, sizeof(double));
	params = (double *)  calloc((size_t) n_params, sizeof(double));
	unit   = (double *)  calloc((size_t) n_vert,   sizeof(double));
	mle    = (double *)  calloc((size_t) n_vert,   sizeof(double));

	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));

	p        = (double **) calloc((size_t)(n_params*n_vert), sizeof(double *));
	for (j = 0; j < n_vert; j++) p[j] = (double *) calloc((size_t) n_params, sizeof(double));

	for (j = 0; j < n_params; j++) start[j] = 45.0;
	
	for (k = 0; k < n_vert; k++) {
		for (j = 0; j < n_params; j++) {
			p[k][j] = start[j];
			if (k == j+1) p[k][j] = op.delta * start[j];
		}
	}

	for (i = 0; i < n_models; i++) {
			create_covariance(nm[i], ts, op.cov_method, Cunits[i], 1);
	}

	if (op.verbose) fprintf(op.fpout, " general_inner_simplex : starting simplex\n");
	if (op.verbose) fprintf(op.fpout, " general_inner_simplex : n_vert = %d\n", n_vert);

	/*************************************/
	/* setup starting values for simplex */
	/*************************************/

	for (k = 0; k < n_vert; k++) {
		for (j = 0; j < n_params; j++) params[j] = p[k][j];

		for (l = 0; l < n_squared; l++) Cfull[l] = 0.0;

		sigma_from_rphi(n_vert,params,unit);

		if (op.verbose) {
			for (l = 0; l < n_vert; l++)   fprintf(op.fpout, " unit[%d] = %f\n", l, unit[l]);
			for (l = 0; l < n_params; l++) fprintf(op.fpout, "  phi[%d] = %f\n", l, params[l]);
			fprintf(stderr, "\n");
		}

		for (i = 0; i < n_models; i++) {

			scale = unit[i]*unit[i];

			time1 = clock();
			for (l = 0; l < n_squared; l++) Cfull[l] += scale*Cunits[i][l];
			time2 = clock();


			if (op.verbose > 1) {
				fprintf(op.fpout," Time to scale covariance matrix ");
				fprintf(op.fpout,"%f seconds\n", ((double) (time2-time1)) / CLOCKS_PER_SEC);
			}
		}
				
		mle[k] = general_MLE(dk, Cfull, op, fit, cov_fit, &r, 1);
	/* if (op.verbose > 0) fprintf(op.fpout, "START : %12.8f %.8f\n", params[0], mle[k]); */
		for (j = 0 ; j < n_params; j++) if (params[j] < 0.0 || params[j] > 90.0) mle[k] -= 1e5; 
		mle[k] = -mle[k];
	}

	for (j = 0; j < n_params; j++) {
		psum[j] = 0.0;
		for (k = 0; k < n_vert; k++) psum[j] += p[k][j];
	}
		
	/*******************************/
	/* Begin the simplex algorithm */
	/*******************************/

	nfun_eval = 0;
	got_answer = 0;
	while ((double) nfun_eval <= op.tol[0]) {

		if (op.verbose) {
			for (k = 0; k < n_vert; k++) {
				fprintf(op.fpout, "%2d %3d %15.8f ", k, nfun_eval,  -mle[k]);
				for (j = 0; j < n_params; j++) fprintf(op.fpout, "%12.8f ", p[k][j]);
				fprintf(op.fpout, "\n");
			}
			fprintf(op.fpout, "\n");
			fflush(op.fpout);
		}
			
		ilo = 0;

		ihi = mle[0] > mle[1] ? (inhi=1,0) : (inhi=0,1);

		for (i = 0; i < n_vert; i++) {
			if (mle[i] <= mle[ilo]) ilo = i;
			if (mle[i] > mle[ihi]) {
				inhi = ihi;
				ihi = i;
			} else if (mle[i] > mle[inhi] && i != ihi) inhi = i;
		}

		/**************************************************/
		/* Check that this stopping criteria is good..... */
		/* It is not the same as the one in NR!           */
		/* Yes it its better but does need tuning         */
		/**************************************************/

		md1 = md2 = md3 = 0.0;
		for (k = 1; k < n_vert; k++) {
			for (j = 0; j < n_params; j++) {
				dist = fabs(p[k][j] - p[0][j]) / p[0][j]; 
				md1 = dist > md1 ? dist : md1;
				dist = fabs(p[k][j] - p[0][j]);
				md2 = dist > md2 ? dist : md2;
			}
			dist = fabs((mle[k] - mle[0]) / mle[0]);
			md3 = dist > md3 ? dist : md3;
		}

		if ((md1 <= op.tol[1] || md2 <= op.tol[3]) && (md3 <= op.tol[2])) {
			if (op.verbose) {
				fprintf(op.fpout, " Stopping criteria met.\n");
				fprintf(op.fpout, " Stopping criteria :(%g <= %g ||  ", md1, op.tol[1]);
				fprintf(op.fpout, " %g <= %g) &&  ", md2, op.tol[3]);
				fprintf(op.fpout, " (%g <= %g)\n", md3, op.tol[2]);
			}
			got_answer = 1;
			break;
		}

		nfun_eval += 2;

		/************************************************************/
		/* first extrapolate by a factor -1 through the face of the */
		/* simplex across from the high point                       */
		/************************************************************/

		fac  = -1.0;
		fac1 = (1.0 - fac) / (double) n_params;
		fac2 = fac1 - fac;

		for (j = 0 ; j < n_params; j++) ptry[j] = psum[j]*fac1 - p[ihi][j]*fac2;

		sigma_from_rphi(n_vert,ptry,unit);

		for (l = 0; l < n_squared; l++) Cfull[l] = 0.0;

		for (i = 0; i < n_models; i++) {
			scale = unit[i]*unit[i];
			for (l = 0; l < n_squared; l++) Cfull[l] += scale*Cunits[i][l];
		}
				
		mle_try = general_MLE(dk, Cfull, op, fit, cov_fit, &r, 1);
		/* if (op.verbose > 0) fprintf(op.fpout, "EXTRA : %12.8f %.8f\n", ptry[0], mle_try); */
		for (j = 0 ; j < n_params; j++) if (ptry[j] < 0.0 || ptry[j] > 90.0) mle_try -= 1e5; 
		mle_try = -mle_try;

		if (mle_try < mle[ihi]) {
			mle[ihi] = mle_try;
			for (j = 0; j < n_params; j++) {
				psum[j] += ptry[j]-p[ihi][j];
				p[ihi][j] = ptry[j];
			}
		}

		/************************************************************/

		if (mle_try < mle[ilo]) {

			fac = 2.0;
			fac1 = (1.0 - fac) / (double) n_params;
			fac2 = fac1 - fac;

			for (j = 0 ; j < n_params; j++) ptry[j] = psum[j]*fac1 - p[ihi][j]*fac2;

			sigma_from_rphi(n_vert,ptry,unit);

			for (l = 0; l < n_squared; l++) Cfull[l] = 0.0;

			for (i = 0; i < n_models; i++) {
				scale = unit[i]*unit[i];
				for (l = 0; l < n_squared; l++) Cfull[l] += scale*Cunits[i][l];
			}

			mle_try = general_MLE(dk, Cfull, op, fit, cov_fit, &r, 1);
			/* if (op.verbose > 0) fprintf(op.fpout, "OTHER1 : %12.8f %.8f\n", ptry[0], mle_try); */
			for (j = 0; j < n_params; j++) if (ptry[j] < 0.0 || ptry[j] > 90.0) mle_try -= 1e5; 
			mle_try = -mle_try;

			if (mle_try < mle[ihi]) {
				mle[ihi] = mle_try;
				for (j = 0; j < n_params; j++) {
					psum[j] += ptry[j]-p[ihi][j];
					p[ihi][j] = ptry[j];
				}
			}
		} else if (mle_try >= mle[inhi]) {
			mle_save = mle[ihi];

			fac = 0.5;
			fac1 = (1.0 - fac) / (double) n_params;
			fac2 = fac1 - fac;
	
			for (j = 0 ; j < n_params; j++) ptry[j] = psum[j]*fac1 - p[ihi][j]*fac2;
			for (l = 0; l < n_squared; l++) Cfull[l] = 0.0;

			sigma_from_rphi(n_vert,ptry,unit);

			for (i = 0; i < n_models; i++) {
				scale = unit[i]*unit[i];
				for (l = 0; l < n_squared; l++) Cfull[l] += scale*Cunits[i][l];
			}
				
			mle_try = general_MLE(dk, Cfull, op, fit, cov_fit, &r, 1);
			/* if (op.verbose > 0) fprintf(op.fpout, "OTHER2 : %12.8f %.8f\n", ptry[0], mle_try); */
			for (j = 0 ; j < n_params; j++) if (ptry[j] < 0.0 || ptry[j] > 90.0) mle_try -= 1e5; 
			mle_try = -mle_try;

			if (mle_try < mle[ihi]) {
				mle[ihi] = mle_try;
				for (j = 0; j < n_params; j++) {
					psum[j] += ptry[j]-p[ihi][j];
					p[ihi][j] = ptry[j];
				}
			}
	
			if (mle_try >= mle_save) {
				for (i = 0; i < n_vert; i++) {
					if (i != ilo) {
						for (j = 0; j < n_params; j++) {
							p[i][j] = psum[j] = 0.5 * (p[i][j] + p[ilo][j]);
						}
						for (j = 0; j < n_params; j++) params[j] = p[i][j];
						for (l = 0; l < n_squared; l++) Cfull[l] = 0.0;

						sigma_from_rphi(n_vert,params,unit);

						for (i = 0; i < n_models; i++) {
							scale = unit[i]*unit[i];
							for (l = 0; l < n_squared; l++) Cfull[l] += scale*Cunits[i][l];
						}
				
						mle_try = general_MLE(dk, Cfull, op, fit, cov_fit, &r, 1);
						/* if (op.verbose > 0) fprintf(op.fpout, "OTHER3 : %12.8f %.8f\n", params[0], mle_try); */
						for (j = 0 ; j < n_params; j++) if (params[j] < 0.0 || params[j] > 0.0) mle_try -= 1e5; 
						mle_try = -mle_try;
					}
				}
			}
			nfun_eval += n_params;
			for (j = 0; j < n_params; j++) {
				psum[j] = 0.0;
				for (k = 0; k < n_vert; k++) psum[j] += p[k][j];
			}
		} else --nfun_eval;
	}

	/**************************/
	/* The end of the simplex */
	/**************************/

	if (got_answer) {

		for (j = 0; j < n_params;  j++) start[j] = p[ilo][j];
		for (l = 0; l < n_squared; l++) Cfull[l] = 0.0;

		sigma_from_rphi(n_vert,start,unit);

		for (i = 0; i < n_models; i++) {
			scale = unit[i]*unit[i];
			for (l = 0; l < n_squared; l++) Cfull[l] += scale*Cunits[i][l];
		}
				
		min_mle = general_MLE(dk, Cfull, op, fit, cov_fit, &r, 1);
		/* if (op.verbose > 0) fprintf(op.fpout, "END : %12.8f %.8f\n", start[0], min_mle); */

		if (end_flag) {

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

			dk.MLE[0] = min_mle;

			calculate_fisher(Cunits, dk, nm, n_models, fit, cov_fit); 

     			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];
				}
			}
		}

	} else {

		fprintf(op.fpout, " Stopping criteria not met.\n");
		fprintf(op.fpout, " Stopping criteria :(%f <= %f ||  ", md1, op.tol[1]);
		fprintf(op.fpout, " %f <= %f) &&  ", md2, op.tol[3]);
		fprintf(op.fpout, " (%f <= %f)\n", md3, op.tol[2]);

		fprintf(stderr, " general_inner_simplex : Could not converge on a solution\n");

		min_mle = -1e10;
		for (j = 0; j < n_params; j++) start[j] = 0.0;
	}
	
	free(ptry);
	free(start);
	free(params);
	free(mle);
	free(psum);
	free(unit);
	free(Cfull);
	free(fit);
	free(cov_fit);

	for (j = 0; j < n_vert; j++) free(p[j]);
	free(p);
	for (j = 0; j < n_models; j++) free(Cunits[j]);
	free(Cunits);

	return(min_mle);
}