Method runMaximumLikelihood

src/backend/worksheet/plots/cartesian/XYFitCurve.cpp:2132–2312

}

void XYFitCurvePrivate::runMaximumLikelihood(const AbstractColumn* tmpXDataColumn, const double norm) {
	const size_t n = tmpXDataColumn->rowCount();

	fitResult.available = true;
	fitResult.valid = true;
	fitResult.status = i18n("Success"); // can it fail in any way?

	const unsigned int np = fitData.paramNames.size(); // number of fit parameters
	fitResult.dof = n - np;
	fitResult.paramValues.resize(np);
	fitResult.errorValues.resize(np);
	fitResult.tdist_tValues.resize(np);
	fitResult.tdist_pValues.resize(np);
	fitResult.marginValues.resize(np);
	fitResult.correlationMatrix.resize(np * (np + 1) / 2);

	DEBUG(Q_FUNC_INFO << ", DISTRIBUTION: " << fitData.modelType)
	fitResult.paramValues[0] = norm; // A - normalization
	// TODO: parameter values (error, etc.)
	// TODO: currently all values are used (data range not changeable)
	const double alpha = 1.0 - fitData.confidenceInterval / 100.;
	const auto& statistics = ((Column*)tmpXDataColumn)->statistics();
	const double mean = statistics.arithmeticMean;
	const double var = statistics.variance;
	const double median = statistics.median;
	const double madmed = statistics.meanDeviationAroundMedian;
	const double iqr = statistics.iqr;
	switch (fitData.modelType) { // only these are supported
	case nsl_sf_stats_gaussian: {
		const double sigma = std::sqrt(var);
		const double mu = mean;
		fitResult.paramValues[1] = sigma;
		fitResult.paramValues[2] = mu;
		DEBUG(Q_FUNC_INFO << ", mu = " << mu << ", sigma = " << sigma)

		fitResult.errorValues[2] = sigma / std::sqrt(n);
		double margin = nsl_stats_tdist_margin(alpha, fitResult.dof, fitResult.errorValues.at(2));
		// DEBUG("z = " << nsl_stats_tdist_z(alpha, fitResult.dof))
		fitResult.marginValues[2] = margin;

		fitResult.errorValues[1] = sigma * sigma / std::sqrt(2 * n);
		margin = nsl_stats_tdist_margin(alpha, fitResult.dof, fitResult.errorValues.at(1));
		// WARN("sigma CONFIDENCE INTERVAL: " << fitResult.paramValues[1] - margin << " .. " << fitResult.paramValues[1] + margin)
		fitResult.marginValues[1] = margin;

		// normalization for spreadsheet or curve
		if (dataSourceType != XYAnalysisCurve::DataSourceType::Histogram)
			fitResult.paramValues[0] /=
				(gsl_sf_erf((tmpXDataColumn->maximum() - mu) / sigma) - gsl_sf_erf((tmpXDataColumn->minimum() - mu) / sigma)) / (2. * std::sqrt(2.));
		break;
	}
	case nsl_sf_stats_exponential: {
		const double mu = tmpXDataColumn->minimum();
		const double lambda = 1. / (mean - mu); // 1/(<x>-\mu)
		fitResult.paramValues[1] = lambda * (1 - 1. / (n - 1)); // unbiased
		fitResult.paramValues[2] = mu;

		fitResult.errorValues[1] = lambda / std::sqrt(n);