cats_simplex.c

最大似然估计算法

		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;
			if (op.speed <= 1) {
				cov_scale(C_unit,dk.n_data,ptry,C);
				linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);
				for (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
			}
			mle_try = MLE_withline_CN(ptry, dk.n_data, Eig_values, Eig_vectors, residuals);

			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;
			if (op.speed == 0) {
				cov_scale(C_unit,dk.n_data,ptry,C);
				linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);
				for (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
			}
			mle_try = MLE_withline_CN(ptry, dk.n_data, Eig_values, Eig_vectors, residuals);


			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];
						if (op.speed == 0) {
							cov_scale(C_unit,dk.n_data,params,C);
							linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);
							for (j = 0; j < dk.n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
						}
						mle[i] = MLE_withline_CN(params, dk.n_data, Eig_values, Eig_vectors, residuals);
					}
				}
			}
			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];

		i = dk.n_par + n_params;
		if (op.speed < 3) {
			cov_scale(C_unit,dk.n_data,start,C);
			linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);

			for (j = 0; j < dk.n_par; j++) params_mle[j]          = fit[j];
			for (j = 0; j < n_params; j++) params_mle[j+dk.n_par] = start[j]; 
		} else {

			for (j = 0; j < dk.n_par; j++) params_mle[j]          = 0.0;
			for (j = 0; j < n_params; j++) params_mle[j+dk.n_par] = start[j]; 

		}

		min_mle = MLE_nparam_CN(params_mle,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values,Eig_vectors);

		if (-min_mle < -cn_mle) {

			if (op.fpout != NULL) { 
				fprintf(op.fpout, " cats_simplex : ");
				fprintf(op.fpout, "mle (%f) is less ", -min_mle);
				fprintf(op.fpout, "than cn only mle (%f)\n", -cn_mle); 
			}
			for (j = 0; j < dk.n_par; j++) params_mle[j] = params_cn[j];
			params_mle[dk.n_par] = 0.0;
			params_mle[dk.n_par+1] = params_cn[dk.n_par];

			min_mle = cn_mle;

			if (return_cov) {
				for (j = 0; j < dk.n_par; j++) {
					for (k = 0; k < dk.n_par; k++) {
						cov_mle[j + k * i] = cov_cn[j + k * (dk.n_par + 1)];
					}
				}
				cov_mle[i*i-1] = cov_cn[(dk.n_par+1)*(dk.n_par+1)];
			}
			
		} else if (-min_mle < wh_mle) {

			if (op.fpout != NULL) { 
				fprintf(op.fpout, " cats_simplex : ");
				fprintf(op.fpout, "mle (%f) is less ", -min_mle);
				fprintf(op.fpout, "than wh only mle (%f)\n", wh_mle); 
			}
			for (j = 0; j < dk.n_par; j++) params_mle[j] = params_wh[j];
			params_mle[dk.n_par] = params_wh[dk.n_par];
			params_mle[dk.n_par+1] = 0.0;

			min_mle = -wh_mle;

			if (return_cov) {
				for (j = 0; j < dk.n_par+1; j++) {
					for (k = 0; k < dk.n_par+1; k++) {
						cov_mle[j + k * i] = cov_wh[j + k * (dk.n_par + 1)];
					}
				}
			}
		} else if (return_cov) {

			cov_scale(C_unit,dk.n_data,start,C);
			linefit_all_fast(dk.A,dk.d,C,dk.n_data,dk.n_par,fit,cov_fit,data_hat);

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

			i = dk.n_par + n_params;
			Ioo = (double *) calloc((size_t)(i*i), sizeof(double) );

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

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

			free(params);

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

			for (j = 0; j < dk.n_par; j++)  work[j]          = params[j]          = fit[j];
			for (j = 0; j < n_params; j++)  work[j+dk.n_par] = params[j+dk.n_par] = start[j];

			/*******************************************************************/
			/* 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 < dk.n_par; j++) {
				for (k = 0; k < dk.n_par; k++) Ioo[j+k*i] = -cov_fit[j+k*dk.n_par];
			}

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

			F2 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

			for (j = dk.n_par; j < i; j++) {
				params[j] = work[j] + 1e-3 * fabs(work[j]);
				F1 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);
				params[j] = work[j] - 1e-3 * fabs(work[j]);
				F3 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

				second_partial = -(F1-2.0*F2+F3)/pow(1e-3*fabs(work[j]),2.0);
				Ioo[j + j * i] = second_partial;
				params[j] = work[j];
			}

			for (j = dk.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_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

					params[j] = work[j] - 1e-3 / 2.0 * fabs(work[j]);
					params[k] = work[k] + 1e-3 / 2.0 * fabs(work[k]);
					F2 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

					params[j] = work[j] + 1e-3 / 2.0 * fabs(work[j]);
					params[k] = work[k] - 1e-3 / 2.0 * fabs(work[k]);
					F3 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

					params[j] = work[j] - 1e-3 / 2.0 * fabs(work[j]);
					params[k] = work[k] - 1e-3 / 2.0 * fabs(work[k]);
					F4 = MLE_nparam_CN(params,i,dk.n_par,dk.n_data,dk.d,dk.A,Eig_values, Eig_vectors);

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

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


			free(work);
	
			/**********************/
			/* Calculate inv(Ioo) */
			/**********************/
	
			lwork = i;
			ipiv = (int *) calloc((size_t) i, sizeof(int));
			work = (double *) calloc((size_t) lwork, sizeof(double));
	
			dgetrf_(&i, &i, Ioo, &i, ipiv, &info);
			dgetri_(&i, Ioo, &i, ipiv, work, &lwork, &info);
	
			free(ipiv);
			free(work);
	
			/***************************/
			/* Store -Ioo into cov_mle */
			/***************************/
	
			for (j = 0; j < i; j++) {
				for (k = 0; k < i; k++) cov_mle[j + k * i] = -Ioo[j + k * i];
			}
			free(Ioo);
		}

	} else {
		fprintf(stderr, " cats_simplex : Minimization did not converge!!\n");

		min_mle = 1e10;
		for (j = 0; j < i; j++)        params_mle[j] = 0.0;
		for (j = 0; j < n_params; j++) start[j]      = 0.0;
	}
	

	free(params);
	free(mle);
	free(psum);
	free(ptry);
	for (j = 0; j < n_vert; j++) free(p[j]);
	free(p);
	free(Eig_vectors);
	free(C);
	free(C_unit);
	free(Eig_values);
	free(residuals);
	free(fit);
	free(cov_fit);
	free(data_hat);

	return(-min_mle);
}