brent_angle.c

最大似然估计算法

		}
	}

	if (got_answer) {

		if (op.verbose) fprintf(op.fpout, " 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]*1000.0);
				fprintf(op.fpout, "cn = %12.6f\n", cn[k]*1000.0); 
			}
			if (value[k] == xmin) {
				start[0] = wh[k];
				start[1] = cn[k];
			}
		}

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

		i = dk.n_par + n_params;

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

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

			if (op.fpout != NULL) {
                                fprintf(op.fpout, " brent_angle : ");
                                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) {

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

                        if (op.fpout != NULL) {
                                fprintf(op.fpout, " brent_angle : ");
                                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)];
                                        }
                                }
                        }

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

                } 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); */

                        i = dk.n_par + n_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(params);
                        free(Ioo);
                }
	} else {
		fprintf(stderr, " brent_angle : 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(mle);
	free(residuals);
	free(start_value);
	free(Eig_vectors);
	free(C);
	free(C_unit);
	free(Eig_values);
	free(fit);
	free(cov_fit);
	free(data_hat);

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

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

	return(-min_mle);
}