brent_color_white.c

最大似然估计算法

	if (got_answer) {

		if (op.verbose) fprintf(op.fpout, "\n Finished : \n\n");
		for (k = 0; k < j; k++) {
			if (op.verbose) {
				fprintf(op.fpout, " Angle = %9.6f ", value[k]); 
				fprintf(op.fpout, "mle = %13.8f ", mle[k]);
				fprintf(op.fpout, "wh = %12.6f ",  wh[k]);
				fprintf(op.fpout, "cn = %12.6f\n", cn[k]); 
			}
			if (value[k] == xmin) {
				start[0] = wh[k];
				start[1] = cn[k];
			}
		}

		if (op.speed < 3) {
			cov_scale(C_unit,n_data,start,C);
			linefit_all_fast(dk.A,dk.d,C,n_data,n_par,fit,cov_fit,data_hat);
			for (j = 0; j < n_data; j++) residuals[j] = dk.d[j] - data_hat[j];
		} else {
			for (j = 0; j < n_data; j++) residuals[j] = dk.d[j];
		}

		i = n_par + n_params;

		min_mle = color_white_mle(start, n_data, Eig_values, Eig_vectors, residuals);

		if (min_mle < cn_mle) {

			/***************************************************************/
			/* Check whether min_mle is lower than coloured_noise only mle */
			/***************************************************************/

			if (op.fpout != NULL) {
				fprintf(op.fpout, " brent_color_white : ");
				fprintf(op.fpout, "mle (%f) is less ", min_mle);
				fprintf(op.fpout, "than color only mle (%f)\n", cn_mle);
			}       

			min_mle = cn_mle;
			nm[wh_flag].sigma[0]  = 0.0;
			nm[cn_flag].sigma[0] = cn_sigma;

			if (end_flag) {

				b = cn_sigma + 1e-3 * fabs(cn_sigma);
				F1 = color_only_mle(&b, n_data, Eig_values, Eig_vectors, res_cn, 0);
	
	       			b = cn_sigma - 1e-3 * fabs(cn_sigma);
	       			F3 = color_only_mle(&b, n_data, Eig_values, Eig_vectors, res_cn, 0);

	       			second_partial = (F1-2.0*cn_mle+F3)/pow(1e-3*fabs(cn_sigma),2.0);

	       			nm[cn_flag].dsigma[0] = sqrt(-1.0 / second_partial);
				nm[wh_flag].dsigma[0] = 0.0;

	       			nm[cn_flag].sigma_flag[0] = 2;
	       			nm[wh_flag].sigma_flag[0] = 2;

				for (j = 0; j < n_par; j++) {
					dk.params[j] = params_cn[j];
					for (k = 0; k < n_par; k++) {
						dk.covar[j + k * n_par] = cov_cn[j + k * n_par] * cn_sigma * cn_sigma;
					}
				}
			}

		} else if (min_mle < wh_mle) {

			/************************************************************/
			/* Check whether min_mle is lower than white_noise only mle */
			/************************************************************/

			if (op.fpout != NULL) {
				fprintf(op.fpout, " brent_color_white : ");
				fprintf(op.fpout, "mle (%f) is less ", min_mle);
				fprintf(op.fpout, "than wh only mle (%f)\n", wh_mle);
			}

			min_mle = wh_mle;
			nm[wh_flag].sigma[0] = wh_sigma;
			nm[cn_flag].sigma[0] = 0.0;

			if (end_flag) {

				b = wh_sigma + 1e-3 * fabs(wh_sigma);
				F1 = MLE_withline_WH(b, n_data, res_wh);

				b = wh_sigma - 1e-3 * fabs(wh_sigma);
				F3 = MLE_withline_WH(b, n_data, res_wh);

	       			second_partial = (F1-2.0*wh_mle+F3)/pow(1e-3*fabs(wh_sigma),2.0);

	       			nm[wh_flag].dsigma[0] = sqrt(-1.0 / second_partial);
	       			nm[cn_flag].dsigma[0] = 0.0;

	       			nm[cn_flag].sigma_flag[0] = 2;
	       			nm[wh_flag].sigma_flag[0] = 2;

				for (j = 0; j < n_par; j++) {
					dk.params[j] = params_wh[j];
					for (k = 0; k < n_par; k++) {
						dk.covar[j + k * n_par] = cov_wh[j + k * n_par] * wh_sigma * wh_sigma;
					}
				}
			}

		/**************************************************/
		/* Okay the minimisation was the best so carry on */
		/**************************************************/

		} else {


			nm[wh_flag].sigma[0] = start[0];
			nm[cn_flag].sigma[0] = start[1];

			if (end_flag) {

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

				lwork = n_par;
				ipiv = (int *)    calloc((size_t) n_par, sizeof(int));
				work = (double *) calloc((size_t) lwork, 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);

				i = n_par + n_params;
				Ioo    = (double *) calloc((size_t)(i*i), sizeof(double) );
				work   = (double *) calloc( (size_t) i,   sizeof(double) );
				params = (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_params; j++)  work[j+n_par] = params[j+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 < 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 */
				/*************************************/
	
				F2 = color_white_mle_res(params,i,n_par,n_data,dk.d,dk.A,Eig_values, Eig_vectors);
	
				for (j = n_par; j < i; j++) {
					params[j] = work[j] + 1e-3 * fabs(work[j]);
					F1 = color_white_mle_res(params,i,n_par,n_data,dk.d,dk.A,Eig_values, Eig_vectors);

					params[j] = work[j] - 1e-3 * fabs(work[j]);
					F3 = color_white_mle_res(params,i,n_par,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 = 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 = color_white_mle_res(params,i,n_par,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 = color_white_mle_res(params,i,n_par,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 = color_white_mle_res(params,i,n_par,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 = color_white_mle_res(params,i,n_par,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);

				/******************************************************/
				/* Now we dont care about the full covariance anymore */  
				/* When is it ever used? so we just need              */
				/* the covaraince of the linear-parameters and        */
				/* the variance of the model sigma's                  */
				/******************************************************/

				j = n_par;
	       			nm[wh_flag].dsigma[0] = sqrt(-Ioo[j+j*i]);


				j = n_par+1;
	       			nm[cn_flag].dsigma[0] = sqrt(-Ioo[j+j*i]);

	       			nm[cn_flag].sigma_flag[0] = 2;
	       			nm[wh_flag].sigma_flag[0] = 2;

				for (j = 0; j < n_par; j++) {
					dk.params[j] = params[j];
					for (k = 0; k < n_par; k++) {
						dk.covar[j + k * n_par] = -Ioo[j + k * i];
					}
				}

				free(params);
				free(Ioo);
			} /* end of end_flag bit */
		} /* end of mle checking */
	} else {

		min_mle = 1e10;
	}

	free(mle);
	free(residuals);
	free(res_wh);
	free(res_cn);
	free(start_value);
	free(Eig_vectors);
	free(Eig_values);
	free(C);
	free(C_unit);
	free(fit);
	free(cov_fit);
	free(data_hat);

	free(cn);
	free(start);
	free(radius);
	free(value);
	free(wh);

	free(cov_wh);
	free(cov_cn);
	free(params_wh);
	free(params_cn);
	free(dy);

	dk.MLE[0] = min_mle;
	return(min_mle);
}