diff --git a/Framework/CurveFitting/CMakeLists.txt b/Framework/CurveFitting/CMakeLists.txt
index 88ca5fabe5d005b2a18387f90c80d30dcaebab98..57f22e7576c99d51c7daaf1c3f0dd972f4d17ef7 100644
--- a/Framework/CurveFitting/CMakeLists.txt
+++ b/Framework/CurveFitting/CMakeLists.txt
@@ -19,6 +19,7 @@ set(SRC_FILES
src/Algorithms/PawleyFit.cpp
src/Algorithms/PlotPeakByLogValue.cpp
src/Algorithms/PlotPeakByLogValueHelper.cpp
+ src/Algorithms/ProfileChiSquared1D.cpp
src/Algorithms/QENSFitSequential.cpp
src/Algorithms/QENSFitSimultaneous.cpp
src/Algorithms/QENSFitUtilities.cpp
@@ -187,6 +188,7 @@ set(INC_FILES
inc/MantidCurveFitting/Algorithms/PawleyFit.h
inc/MantidCurveFitting/Algorithms/PlotPeakByLogValue.h
inc/MantidCurveFitting/Algorithms/PlotPeakByLogValueHelper.h
+ inc/MantidCurveFitting/Algorithms/ProfileChiSquared1D.h
inc/MantidCurveFitting/Algorithms/QENSFitSequential.h
inc/MantidCurveFitting/Algorithms/QENSFitSimultaneous.h
inc/MantidCurveFitting/Algorithms/QENSFitUtilities.h
@@ -351,6 +353,7 @@ set(TEST_FILES
Algorithms/PawleyFitTest.h
Algorithms/PlotPeakByLogValueTest.h
Algorithms/PlotPeakByLogValueHelperTest.h
+ Algorithms/ProfileChiSquared1DTest.h
Algorithms/QENSFitSequentialTest.h
Algorithms/QENSFitSimultaneousTest.h
Algorithms/RefinePowderInstrumentParameters3Test.h
diff --git a/Framework/CurveFitting/inc/MantidCurveFitting/Algorithms/CalculateChiSquared.h b/Framework/CurveFitting/inc/MantidCurveFitting/Algorithms/CalculateChiSquared.h
index bdf47c12734e9a43646f7db517acf3e9986c0b5e..31ea00716bcee331032d846844ca802ff92fa618 100644
--- a/Framework/CurveFitting/inc/MantidCurveFitting/Algorithms/CalculateChiSquared.h
+++ b/Framework/CurveFitting/inc/MantidCurveFitting/Algorithms/CalculateChiSquared.h
@@ -6,6 +6,7 @@
// SPDX - License - Identifier: GPL - 3.0 +
#pragma once
+#include "MantidAPI/IFunction_fwd.h"
#include "MantidCurveFitting/IFittingAlgorithm.h"
#include "MantidKernel/System.h"
@@ -14,29 +15,24 @@ namespace CurveFitting {
namespace Algorithms {
/**
-
Calculate chi squared for a function and a data set in a workspace.
- Optionally outputs slices of the chi^2 along the parameter axes
- and estimates the standard deviations.
*/
class MANTID_CURVEFITTING_DLL CalculateChiSquared : public IFittingAlgorithm {
public:
const std::string name() const override;
int version() const override;
const std::vector<std::string> seeAlso() const override {
- return {"CalculateCostFunction", "Fit"};
+ return {"CalculateCostFunction", "Fit", "ProfileChiSquared1D"};
}
const std::string summary() const override;
+ static void calcChiSquared(const API::IFunction &fun, size_t nParams,
+ const API::FunctionDomain &domain,
+ API::FunctionValues &values, double &chiSquared,
+ double &chiSquaredWeighted, double &dof);
private:
void initConcrete() override;
void execConcrete() override;
- void estimateErrors();
- void unfixParameters();
- void refixParameters();
-
- /// Cache indices of fixed parameters
- std::vector<size_t> m_fixedParameters;
};
} // namespace Algorithms
diff --git a/Framework/CurveFitting/inc/MantidCurveFitting/Algorithms/ProfileChiSquared1D.h b/Framework/CurveFitting/inc/MantidCurveFitting/Algorithms/ProfileChiSquared1D.h
new file mode 100644
index 0000000000000000000000000000000000000000..762995aff93a7e6e16f97258e7950bb2e8cb116a
--- /dev/null
+++ b/Framework/CurveFitting/inc/MantidCurveFitting/Algorithms/ProfileChiSquared1D.h
@@ -0,0 +1,47 @@
+// Mantid Repository : https://github.com/mantidproject/mantid
+//
+// Copyright © 2021 ISIS Rutherford Appleton Laboratory UKRI,
+// NScD Oak Ridge National Laboratory, European Spallation Source,
+// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
+// SPDX - License - Identifier: GPL - 3.0 +
+#pragma once
+
+#include "MantidCurveFitting/Functions/ChebfunBase.h"
+#include "MantidCurveFitting/GSLMatrix.h"
+#include "MantidCurveFitting/IFittingAlgorithm.h"
+#include "MantidKernel/System.h"
+
+namespace Mantid {
+namespace CurveFitting {
+namespace Algorithms {
+
+/**
+Profiles chi2 about its minimum to find parameter errors
+*/
+class MANTID_CURVEFITTING_DLL ProfileChiSquared1D : public IFittingAlgorithm {
+public:
+ ProfileChiSquared1D();
+ const std::string name() const override;
+ int version() const override;
+ const std::vector<std::string> seeAlso() const override {
+ return {"CalculateChiSquared", "Fit"};
+ }
+ const std::string summary() const override;
+
+private:
+ void initConcrete() override;
+ void execConcrete() override;
+ void unfixParameters();
+ void refixParameters();
+ GSLMatrix getCovarianceMatrix();
+ std::tuple<double, double>
+ getChiSquaredRoots(const Functions::ChebfunBase_sptr &approximation,
+ std::vector<double> &coeffs, double qvalue, double rBound,
+ double lBound) const;
+ /// Cache indices of fixed parameters
+ std::vector<size_t> m_fixedParameters;
+};
+
+} // namespace Algorithms
+} // namespace CurveFitting
+} // namespace Mantid
diff --git a/Framework/CurveFitting/src/Algorithms/CalculateChiSquared.cpp b/Framework/CurveFitting/src/Algorithms/CalculateChiSquared.cpp
index 290ed48e4923513773e12b95e32526b2d45ac834..26812a161aa0390a24c5bf576ec56aa41cfa76ec 100644
--- a/Framework/CurveFitting/src/Algorithms/CalculateChiSquared.cpp
+++ b/Framework/CurveFitting/src/Algorithms/CalculateChiSquared.cpp
@@ -5,12 +5,6 @@
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidCurveFitting/Algorithms/CalculateChiSquared.h"
-#include "MantidAPI/Column.h"
-#include "MantidAPI/ITableWorkspace.h"
-#include "MantidAPI/TableRow.h"
-#include "MantidAPI/WorkspaceFactory.h"
-#include "MantidCurveFitting/Functions/ChebfunBase.h"
-#include "MantidCurveFitting/GSLJacobian.h"
namespace Mantid {
namespace CurveFitting {
@@ -18,51 +12,10 @@ namespace Algorithms {
using namespace Kernel;
using namespace API;
-using namespace Functions;
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(CalculateChiSquared)
-//----------------------------------------------------------------------------------------------
-namespace {
-
-/// Caclculate the chi squared, weighted chi squared and the number of degrees
-/// of freedom.
-/// @param domain :: Function's domain.
-/// @param nParams :: Number of free fitting parameters.
-/// @param values :: Functin's values.
-/// @param chi0 :: Chi squared at the minimum.
-/// @param sigma2 :: Estimated variance of the fitted data.
-void calcChiSquared(const API::IFunction &fun, size_t nParams,
- const API::FunctionDomain &domain,
- API::FunctionValues &values, double &chiSquared,
- double &chiSquaredWeighted, double &dof) {
-
- // Calculate function values.
- fun.function(domain, values);
-
- // Calculate the chi squared.
- chiSquared = 0.0;
- chiSquaredWeighted = 0.0;
- dof = -static_cast<double>(nParams);
- for (size_t i = 0; i < values.size(); ++i) {
- auto weight = values.getFitWeight(i);
- if (weight > 0.0) {
- double tmp = values.getFitData(i) - values.getCalculated(i);
- chiSquared += tmp * tmp;
- tmp *= weight;
- chiSquaredWeighted += tmp * tmp;
- dof += 1.0;
- }
- }
- if (dof <= 0.0) {
- dof = 1.0;
- }
-}
-} // namespace
-
-//----------------------------------------------------------------------------------------------
-
/// Algorithms name for identification. @see Algorithm::name
const std::string CalculateChiSquared::name() const {
return "CalculateChiSquared";
@@ -101,7 +54,6 @@ void CalculateChiSquared::initConcrete() {
"Output value of weighted chi squared divided by the "
"number of data points).",
Direction::Output);
- declareProperty("Output", "", "A base name for output workspaces.");
declareProperty("Weighted", false,
"Option to use the weighted chi squared "
"in error estimation. Default is false.");
@@ -172,512 +124,43 @@ void CalculateChiSquared::execConcrete() {
// Store the result.
setProperty("ChiSquaredDividedByDOF", chiSquaredDOF);
setProperty("ChiSquaredWeightedDividedByDOF", chiSquaredWeightedDOF);
-
- std::string baseName = getProperty("Output");
- if (!baseName.empty()) {
- estimateErrors();
- }
}
//----------------------------------------------------------------------------------------------
-namespace {
-/// Calculate the negative logarithm of the probability density function (PDF):
-/// if a = getDiff(...) then exp(-a) is a value of the PDF.
-/// @param domain :: Function's domain.
+/// Caclculate the chi squared, weighted chi squared and the number of degrees
+/// of freedom.
+/// @param fun :: Function's domain.
/// @param nParams :: Number of free fitting parameters.
-/// @param values :: Functin's values.
-/// @param chi0 :: Chi squared at the minimum.
-/// @param sigma2 :: Estimated variance of the fitted data.
-double getDiff(const API::IFunction &fun, size_t nParams,
- const API::FunctionDomain &domain, API::FunctionValues &values,
- double chi0, double sigma2) {
- double chiSquared = 0.0;
- double chiSquaredWeighted = 0.0;
- double dof = 0;
- calcChiSquared(fun, nParams, domain, values, chiSquared, chiSquaredWeighted,
- dof);
- double res = 0.0;
- if (sigma2 > 0) {
- res = (chiSquared - chi0) / 2 / sigma2;
- } else {
- res = (chiSquaredWeighted - chi0) / 2;
- }
- return res;
-}
-
-/// Helper class to calculate the chi squared along a direction in the parameter
-/// space.
-class ChiSlice {
-public:
- /// Constructor.
- /// @param f :: The fitting function
- /// @param dir :: A normalised direction vector in the parameter space.
- /// @param domain :: Function's domain.
- /// @param values :: Functin's values.
- /// @param chi0 :: Chi squared at the minimum.
- /// @param sigma2 :: Estimated variance of the fitted data.
- ChiSlice(IFunction &f, const GSLVector &dir,
- const API::FunctionDomain &domain, API::FunctionValues &values,
- double chi0, double sigma2)
- : m_function(f), m_direction(dir), m_domain(domain), m_values(values),
- m_chi0(chi0), m_sigma2(sigma2) {}
- /// Calculate the value of chi squared along the chosen direction at a
- /// distance from
- /// the minimum point.
- /// @param p :: A distance from the minimum.
- double operator()(double p) {
- std::vector<double> par0(m_function.nParams());
- for (size_t ip = 0; ip < m_function.nParams(); ++ip) {
- par0[ip] = m_function.getParameter(ip);
- m_function.setParameter(ip, par0[ip] + p * m_direction[ip]);
- }
- double res = getDiff(m_function, m_function.nParams(), m_domain, m_values,
- m_chi0, m_sigma2);
- for (size_t ip = 0; ip < m_function.nParams(); ++ip) {
- m_function.setParameter(ip, par0[ip]);
- }
- return res;
- }
-
- /// Make an approximation for this slice on an interval.
- /// @param lBound :: The left bound of the approximation interval.
- /// @param rBound :: The right bound of the approximation interval.
- /// @param P :: Output vector with approximation parameters.
- /// @param A :: Output vector with approximation parameters.
- ChebfunBase_sptr makeApprox(double lBound, double rBound,
- std::vector<double> &P, std::vector<double> &A,
- bool &ok) {
-
- auto base = ChebfunBase::bestFitAnyTolerance(lBound, rBound, *this, P, A,
- 1.0, 1e-4, 129);
- ok = bool(base);
- if (!base) {
- base = std::make_shared<ChebfunBase>(10, lBound, rBound, 1e-4);
- P = base->fit(*this);
- A = base->calcA(P);
- }
- return base;
- }
-
- /// Fiind a displacement in the parameter space from the initial point
- /// to a point where the PDF drops significantly.
- /// @param shift :: Initial shift form par0 value.
- double findBound(double shift) {
- double bound0 = 0;
- double diff0 = (*this)(0);
- double bound = shift;
- bool canDecrease = true;
- for (size_t i = 0; i < 100; ++i) {
- double diff = (*this)(bound);
-
- bool isIncreasing = fabs(bound) > fabs(bound0) && diff > diff0;
- if (canDecrease) {
- if (isIncreasing)
- canDecrease = false;
- } else {
- if (!isIncreasing) {
- bound = bound0;
- break;
- }
- }
-
- bound0 = bound;
- diff0 = diff;
-
- if (diff > 3.0) {
- if (diff < 4.0) {
- break;
- }
- // diff is too large
- bound *= 0.75;
- } else {
- // diff is too small
- bound *= 2;
- }
- }
- return bound;
- }
-
-private:
- /// The fitting function
- IFunction &m_function;
- /// The direction in the parameter space
- GSLVector m_direction;
- /// The domain
- const API::FunctionDomain &m_domain;
- /// The values
- API::FunctionValues &m_values;
- /// The chi squared at the minimum
- double m_chi0;
- /// The data variance.
- double m_sigma2;
-};
-} // namespace
-
-//----------------------------------------------------------------------------------------------
-/// Examine the chi squared as a function of fitting parameters and estimate
-/// errors for each parameter.
-void CalculateChiSquared::estimateErrors() {
- // Number of fiting parameters
- auto nParams = m_function->nParams();
- // Create an output table for displaying slices of the chi squared and
- // the probabilitydensity function
- auto pdfTable = API::WorkspaceFactory::Instance().createTable();
-
- std::string baseName = getProperty("Output");
- if (baseName.empty()) {
- baseName = "CalculateChiSquared";
- }
- declareProperty(
- std::make_unique<API::WorkspaceProperty<API::ITableWorkspace>>(
- "PDFs", "", Kernel::Direction::Output),
- "The name of the TableWorkspace in which to store the "
- "pdfs of fit parameters");
- setPropertyValue("PDFs", baseName + "_pdf");
- setProperty("PDFs", pdfTable);
-
- // Create an output table for displaying the parameter errors.
- auto errorsTable = API::WorkspaceFactory::Instance().createTable();
- auto nameColumn = errorsTable->addColumn("str", "Parameter");
- auto valueColumn = errorsTable->addColumn("double", "Value");
- auto minValueColumn = errorsTable->addColumn("double", "Value at Min");
- auto leftErrColumn = errorsTable->addColumn("double", "Left Error");
- auto rightErrColumn = errorsTable->addColumn("double", "Right Error");
- auto quadraticErrColumn = errorsTable->addColumn("double", "Quadratic Error");
- auto chiMinColumn = errorsTable->addColumn("double", "Chi2 Min");
- errorsTable->setRowCount(nParams);
- declareProperty(
- std::make_unique<API::WorkspaceProperty<API::ITableWorkspace>>(
- "Errors", "", Kernel::Direction::Output),
- "The name of the TableWorkspace in which to store the "
- "values and errors of fit parameters");
- setPropertyValue("Errors", baseName + "_errors");
- setProperty("Errors", errorsTable);
-
- // Calculate initial values
- double chiSquared = 0.0;
- double chiSquaredWeighted = 0.0;
- double dof = 0;
- API::FunctionDomain_sptr domain;
- API::FunctionValues_sptr values;
- m_domainCreator->createDomain(domain, values);
- calcChiSquared(*m_function, nParams, *domain, *values, chiSquared,
- chiSquaredWeighted, dof);
- // Value of chi squared for current parameters in m_function
- double chi0 = chiSquared;
- // Fit data variance
- double sigma2 = chiSquared / dof;
- bool useWeighted = getProperty("Weighted");
+/// @param domain :: Function's domain
+/// @param values :: Function's values
+/// @param chiSquared :: unweighted chi squared
+/// @param chiSquaredWeighted :: weighted chi squared
+/// @param dof :: degrees of freedom
+void CalculateChiSquared::calcChiSquared(
+ const API::IFunction &fun, size_t nParams,
+ const API::FunctionDomain &domain, API::FunctionValues &values,
+ double &chiSquared, double &chiSquaredWeighted, double &dof) {
- if (useWeighted) {
- chi0 = chiSquaredWeighted;
- sigma2 = 0.0;
- }
-
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << "chi0=" << chi0 << '\n';
- g_log.debug() << "sigma2=" << sigma2 << '\n';
- g_log.debug() << "dof=" << dof << '\n';
- }
-
- // Parameter bounds that define a volume in the parameter
- // space within which the chi squared is being examined.
- GSLVector lBounds(nParams);
- GSLVector rBounds(nParams);
-
- // Number of points in lines for plotting
- size_t n = 100;
- pdfTable->setRowCount(n);
- const double fac = 1e-4;
-
- // Loop over each parameter
- for (size_t ip = 0; ip < nParams; ++ip) {
-
- // Add columns for the parameter to the pdf table.
- auto parName = m_function->parameterName(ip);
- nameColumn->read(ip, parName);
- // Parameter values
- auto col1 = pdfTable->addColumn("double", parName);
- col1->setPlotType(1);
- // Chi squared values
- auto col2 = pdfTable->addColumn("double", parName + "_chi2");
- col2->setPlotType(2);
- // PDF values
- auto col3 = pdfTable->addColumn("double", parName + "_pdf");
- col3->setPlotType(2);
-
- double par0 = m_function->getParameter(ip);
- double shift = fabs(par0 * fac);
- if (shift == 0.0) {
- shift = fac;
- }
-
- // Make a slice along this parameter
- GSLVector dir(nParams);
- dir.zero();
- dir[ip] = 1.0;
- ChiSlice slice(*m_function, dir, *domain, *values, chi0, sigma2);
-
- // Find the bounds withn which the PDF is significantly above zero.
- // The bounds are defined relative to par0:
- // par0 + lBound is the lowest value of the parameter (lBound <= 0)
- // par0 + rBound is the highest value of the parameter (rBound >= 0)
- double lBound = slice.findBound(-shift);
- double rBound = slice.findBound(shift);
- lBounds[ip] = lBound;
- rBounds[ip] = rBound;
-
- // Approximate the slice with a polynomial.
- // P is a vector of values of the polynomial at special points.
- // A is a vector of Chebyshev expansion coefficients.
- // The polynomial is defined on interval [lBound, rBound]
- // The value of the polynomial at 0 == chi squared at par0
- std::vector<double> P, A;
- bool ok = true;
- auto base = slice.makeApprox(lBound, rBound, P, A, ok);
- if (!ok) {
- g_log.warning() << "Approximation failed for parameter " << ip << '\n';
- }
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << "Parameter " << ip << '\n';
- g_log.debug() << "Slice approximated by polynomial of order "
- << base->size() - 1;
- g_log.debug() << " between " << lBound << " and " << rBound << '\n';
- }
-
- // Write n slice points into the output table.
- double dp = (rBound - lBound) / static_cast<double>(n);
- for (size_t i = 0; i < n; ++i) {
- double par = lBound + dp * static_cast<double>(i);
- double chi = base->eval(par, P);
- col1->fromDouble(i, par0 + par);
- col2->fromDouble(i, chi);
- }
-
- // Check if par0 is a minimum point of the chi squared
- std::vector<double> AD;
- // Calculate the derivative polynomial.
- // AD are the Chebyshev expansion of the derivative.
- base->derivative(A, AD);
- // Find the roots of the derivative polynomial
- std::vector<double> minima = base->roots(AD);
- if (minima.empty()) {
- minima.emplace_back(par0);
- }
-
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << "Minima: ";
- }
-
- // If only 1 extremum is found assume (without checking) that it's a
- // minimum.
- // If there are more than 1, find the one with the smallest chi^2.
- double chiMin = std::numeric_limits<double>::max();
- double parMin = par0;
- for (double minimum : minima) {
- double value = base->eval(minimum, P);
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << minimum << " (" << value << ") ";
- }
- if (value < chiMin) {
- chiMin = value;
- parMin = minimum;
- }
- }
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << '\n';
- g_log.debug() << "Smallest minimum at " << parMin << " is " << chiMin
- << '\n';
- }
-
- // Points of intersections with line chi^2 = 1/2 give an estimate of
- // the standard deviation of this parameter if it's uncorrelated with the
- // others.
- A[0] -= 0.5; // Now A are the coefficients of the original polynomial
- // shifted down by 1/2.
- std::vector<double> roots = base->roots(A);
- std::sort(roots.begin(), roots.end());
-
- if (roots.empty()) {
- // Something went wrong; use the whole interval.
- roots.resize(2);
- roots[0] = lBound;
- roots[1] = rBound;
- } else if (roots.size() == 1) {
- // Only one root found; use a bound for the other root.
- if (roots.front() < 0) {
- roots.emplace_back(rBound);
- } else {
- roots.insert(roots.begin(), lBound);
- }
- } else if (roots.size() > 2) {
- // More than 2 roots; use the smallest and the biggest
- auto smallest = roots.front();
- auto biggest = roots.back();
- roots.resize(2);
- roots[0] = smallest;
- roots[1] = biggest;
- }
-
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << "Roots: ";
- for (double root : roots) {
- g_log.debug() << root << ' ';
- }
- g_log.debug() << '\n';
- }
-
- // Output parameter info to the table.
- valueColumn->fromDouble(ip, par0);
- minValueColumn->fromDouble(ip, par0 + parMin);
- leftErrColumn->fromDouble(ip, roots[0] - parMin);
- rightErrColumn->fromDouble(ip, roots[1] - parMin);
- chiMinColumn->fromDouble(ip, chiMin);
-
- // Output the PDF
- for (size_t i = 0; i < n; ++i) {
- double chi = col2->toDouble(i);
- col3->fromDouble(i, exp(-chi + chiMin));
- }
-
- // make sure function parameters don't change.
- m_function->setParameter(ip, par0);
- }
-
- // Improve estimates for standard deviations.
- // If parameters are correlated the found deviations
- // most likely underestimate the true values.
- unfixParameters();
- GSLJacobian J(*m_function, values->size());
- m_function->functionDeriv(*domain, J);
- refixParameters();
- // Calculate the hessian at the current point.
- GSLMatrix H;
- if (useWeighted) {
- H.resize(nParams, nParams);
- for (size_t i = 0; i < nParams; ++i) {
- for (size_t j = i; j < nParams; ++j) {
- double h = 0.0;
- for (size_t k = 0; k < values->size(); ++k) {
- double w = values->getFitWeight(k);
- h += J.get(k, i) * J.get(k, j) * w * w;
- }
- H.set(i, j, h);
- if (i != j) {
- H.set(j, i, h);
- }
- }
- }
- } else {
- H = J.matrix().tr() * J.matrix();
- }
- // Square roots of the diagonals of the covariance matrix give
- // the standard deviations in the quadratic approximation of the chi^2.
- GSLMatrix V(H);
- if (!useWeighted) {
- V *= 1. / sigma2;
- }
- V.invert();
- // In a non-quadratic asymmetric case the following procedure can give a
- // better result:
- // Find the direction in which the chi^2 changes slowest and the positive and
- // negative deviations in that direction. The change in a parameter at those
- // points can be a better estimate for the standard deviation.
- GSLVector v(nParams);
- GSLMatrix Q(nParams, nParams);
- // One of the eigenvectors of the hessian is the direction of the slowest
- // change.
- H.eigenSystem(v, Q);
-
- // Loop over the eigenvectors
- for (size_t i = 0; i < nParams; ++i) {
- auto dir = Q.copyColumn(i);
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << "Direction " << i << '\n';
- g_log.debug() << dir << '\n';
- }
- // Make a slice in that direction
- ChiSlice slice(*m_function, dir, *domain, *values, chi0, sigma2);
- double rBound0 = dir.dot(rBounds);
- double lBound0 = dir.dot(lBounds);
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << "lBound " << lBound0 << '\n';
- g_log.debug() << "rBound " << rBound0 << '\n';
- }
- double lBound = slice.findBound(lBound0);
- double rBound = slice.findBound(rBound0);
- std::vector<double> P, A;
- // Use a polynomial approximation
- bool ok = true;
- auto base = slice.makeApprox(lBound, rBound, P, A, ok);
- if (!ok) {
- g_log.warning() << "Approximation failed in direction " << i << '\n';
- }
- // Find the deviation points where the chi^2 = 1/2
- A[0] -= 0.5;
- std::vector<double> roots = base->roots(A);
- std::sort(roots.begin(), roots.end());
- // Sort out the roots
- auto nRoots = roots.size();
- if (nRoots == 0) {
- roots.resize(2, 0.0);
- } else if (nRoots == 1) {
- if (roots.front() > 0.0) {
- roots.insert(roots.begin(), 0.0);
- } else {
- roots.emplace_back(0.0);
- }
- } else if (nRoots > 2) {
- roots[1] = roots.back();
- roots.resize(2);
- }
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << "Roots " << roots[0] << " (" << slice(roots[0]) << ") "
- << roots[1] << " (" << slice(roots[1]) << ") \n";
- }
- // Loop over the parameters and see if there deviations along
- // this direction is greater than any previous value.
- for (size_t ip = 0; ip < nParams; ++ip) {
- auto lError = roots.front() * dir[ip];
- auto rError = roots.back() * dir[ip];
- if (lError > rError) {
- std::swap(lError, rError);
- }
- if (lError < leftErrColumn->toDouble(ip)) {
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << " left for " << ip << ' ' << lError << ' '
- << leftErrColumn->toDouble(ip) << '\n';
- }
- leftErrColumn->fromDouble(ip, lError);
- }
- if (rError > rightErrColumn->toDouble(ip)) {
- if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
- g_log.debug() << " right for " << ip << ' ' << rError << ' '
- << rightErrColumn->toDouble(ip) << '\n';
- }
- rightErrColumn->fromDouble(ip, rError);
- }
- }
- // Output the quadratic estimate for comparrison.
- quadraticErrColumn->fromDouble(i, sqrt(V.get(i, i)));
- }
-}
+ // Calculate function values.
+ fun.function(domain, values);
-/// Temporary unfix any fixed parameters.
-void CalculateChiSquared::unfixParameters() {
- for (size_t i = 0; i < m_function->nParams(); ++i) {
- if (!m_function->isActive(i)) {
- m_function->unfix(i);
- m_fixedParameters.emplace_back(i);
+ // Calculate the chi squared.
+ chiSquared = 0.0;
+ chiSquaredWeighted = 0.0;
+ dof = -static_cast<double>(nParams);
+ for (size_t i = 0; i < values.size(); ++i) {
+ auto weight = values.getFitWeight(i);
+ if (weight > 0.0) {
+ double tmp = values.getFitData(i) - values.getCalculated(i);
+ chiSquared += tmp * tmp;
+ tmp *= weight;
+ chiSquaredWeighted += tmp * tmp;
+ dof += 1.0;
}
}
-}
-
-/// Restore the "fixed" status of previously unfixed paramters.
-void CalculateChiSquared::refixParameters() {
- for (auto &fixedParameter : m_fixedParameters) {
- m_function->fix(fixedParameter);
+ if (dof <= 0.0) {
+ dof = 1.0;
}
}
diff --git a/Framework/CurveFitting/src/Algorithms/ProfileChiSquared1D.cpp b/Framework/CurveFitting/src/Algorithms/ProfileChiSquared1D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d4a624d8b01f9d084da6e6aba55a88600daebd5
--- /dev/null
+++ b/Framework/CurveFitting/src/Algorithms/ProfileChiSquared1D.cpp
@@ -0,0 +1,491 @@
+// Mantid Repository : https://github.com/mantidproject/mantid
+//
+// Copyright © 2021 ISIS Rutherford Appleton Laboratory UKRI,
+// NScD Oak Ridge National Laboratory, European Spallation Source,
+// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
+// SPDX - License - Identifier: GPL - 3.0 +
+#include "MantidCurveFitting/Algorithms/ProfileChiSquared1D.h"
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidAPI/Column.h"
+#include "MantidAPI/IFunction.h"
+#include "MantidAPI/ITableWorkspace.h"
+#include "MantidAPI/MatrixWorkspace.h"
+#include "MantidAPI/WorkspaceFactory.h"
+#include "MantidCurveFitting/Algorithms/CalculateChiSquared.h"
+#include "MantidCurveFitting/Algorithms/Fit.h"
+#include "MantidCurveFitting/GSLJacobian.h"
+
+#include <boost/math/distributions/chi_squared.hpp>
+
+namespace {
+// The maximum difference of chi squared to search for
+// 10.8276 covers 99.9% of the distrubition
+constexpr double MAXCHISQUAREDIFFERENCE = 10.8276;
+
+/// Calculate the change in chi2
+/// @param domain :: Function's domain.
+/// @param nParams :: Number of free fitting parameters.
+/// @param values :: Functin's values.
+/// @param chi0 :: Chi squared at the minimum.
+double getDiff(const Mantid::API::IFunction &fun, size_t nParams,
+ const Mantid::API::FunctionDomain &domain,
+ Mantid::API::FunctionValues &values, double chi0) {
+ double chiSquared = 0.0;
+ double chiSquaredWeighted = 0.0;
+ double dof = 0;
+ Mantid::CurveFitting::Algorithms::CalculateChiSquared::calcChiSquared(
+ fun, nParams, domain, values, chiSquared, chiSquaredWeighted, dof);
+ return chiSquaredWeighted - chi0;
+}
+
+} // namespace
+
+namespace Mantid {
+namespace CurveFitting {
+namespace Algorithms {
+
+DECLARE_ALGORITHM(ProfileChiSquared1D)
+
+using namespace Mantid::API;
+/// Helper class to calculate the chi squared along a direction in the parameter
+/// space.
+class ChiSlice {
+public:
+ /// Constructor.
+ /// @param inputFunction :: The fitting function
+ /// @param fixedParameterIndex :: index of the parameter which is fixed
+ /// @param inputWS :: The input workspace (used for fit algorithm)
+ /// @param workspaceIndex :: Workspace index (used for fit algorithm)
+ /// @param domain :: Function's domain.
+ /// @param values :: Functin's values.
+ /// @param chi0 :: Chi squared at the minimum.
+ /// @param freeParameters :: Parameters which are free in the function.
+ ChiSlice(IFunction_sptr inputFunction, int fixedParameterIndex,
+ API::MatrixWorkspace_sptr inputWS, int workspaceIndex,
+ const API::FunctionDomain &domain, API::FunctionValues &values,
+ double chi0, std::vector<int> &freeParameters)
+ : m_fixedParameterIndex(fixedParameterIndex), m_domain(domain),
+ m_values(values), m_chi0(chi0), m_function(inputFunction),
+ m_ws(inputWS), m_workspaceIndex(workspaceIndex),
+ m_freeParameters(freeParameters) {
+ // create a fitting algorithm based on least squares (which is the default)
+ m_fitalg = AlgorithmFactory::Instance().create("Fit", -1);
+ m_fitalg->setChild(true);
+ }
+ /// Calculate the value of chi squared along the chosen direction at a
+ /// distance from
+ /// the minimum point.
+ /// @param p :: A distance from the minimum.
+ double operator()(double p) {
+ m_fitalg->initialize();
+ m_fitalg->setProperty("Function", m_function);
+ m_fitalg->setProperty("InputWorkspace", m_ws);
+ m_fitalg->setProperty("WorkspaceIndex", m_workspaceIndex);
+ IFunction_sptr function = m_fitalg->getProperty("Function");
+ std::vector<double> originalParamValues(function->nParams());
+ for (auto ip = 0u; ip < function->nParams(); ++ip) {
+ originalParamValues[ip] = function->getParameter(ip);
+ }
+ function->setParameter(m_fixedParameterIndex,
+ originalParamValues[m_fixedParameterIndex] + p);
+ function->fix(m_fixedParameterIndex);
+
+ // re run the fit to minimze the unfixed parameters
+ m_fitalg->execute();
+ // find change in chi 2
+ // num free parameters is the number of global free parameters - the 1 we've
+ // just fixed
+ int numFreeParameters = static_cast<int>(m_freeParameters.size() - 1);
+ double res =
+ getDiff(*function, numFreeParameters, m_domain, m_values, m_chi0);
+ // reset fit to original values
+ for (auto ip = 0u; ip < function->nParams(); ++ip) {
+ function->setParameter(ip, originalParamValues[ip]);
+ }
+ function->unfix(m_fixedParameterIndex);
+ return res;
+ }
+
+ /// Make an approximation for this slice on an interval.
+ /// @param lBound :: The left bound of the approximation interval.
+ /// @param rBound :: The right bound of the approximation interval.
+ /// @param P :: Output vector with approximation parameters.
+ /// @param A :: Output vector with approximation parameters.
+ Functions::ChebfunBase_sptr makeApprox(double lBound, double rBound,
+ std::vector<double> &P,
+ std::vector<double> &A) {
+
+ auto base = Functions::ChebfunBase::bestFitAnyTolerance(
+ lBound, rBound, *this, P, A, 1.0, 1e-4, 129);
+ if (!base) {
+ base = std::make_shared<Functions::ChebfunBase>(10, lBound, rBound, 1e-4);
+ P = base->fit(*this);
+ A = base->calcA(P);
+ }
+ return base;
+ }
+
+ /// Fiind a displacement in the parameter space from the initial point
+ /// to a point where the PDF drops significantly.
+ /// @param shift :: Initial shift form par0 value.
+ double findBound(double shift) {
+ double bound0 = 0;
+ double diff0 = (*this)(0);
+ double bound = shift;
+ bool canDecrease = true;
+ for (size_t i = 0; i < 100; ++i) {
+ double diff = (*this)(bound);
+
+ bool isIncreasing = fabs(bound) > fabs(bound0) && diff > diff0;
+ if (canDecrease) {
+ if (isIncreasing)
+ canDecrease = false;
+ } else {
+ if (!isIncreasing) {
+ bound = bound0;
+ break;
+ }
+ }
+
+ bound0 = bound;
+ diff0 = diff;
+
+ if (diff > MAXCHISQUAREDIFFERENCE - 1) {
+ if (diff < MAXCHISQUAREDIFFERENCE) {
+ break;
+ }
+ // diff is too large
+ bound *= 0.75;
+ } else {
+ // diff is too small
+ bound *= 2;
+ }
+ }
+ return bound;
+ }
+
+private:
+ // Fixed parameter index
+ int m_fixedParameterIndex;
+ /// The domain
+ const API::FunctionDomain &m_domain;
+ /// The values
+ API::FunctionValues &m_values;
+ /// The chi squared at the minimum
+ double m_chi0;
+ // fitting algorithm
+ IAlgorithm_sptr m_fitalg;
+ // Input workspace and function
+ IFunction_sptr m_function;
+ MatrixWorkspace_sptr m_ws;
+ int m_workspaceIndex;
+ // Vector of free parameter indices
+ std::vector<int> m_freeParameters;
+}; // namespace Algorithms
+
+/// Default constructor
+ProfileChiSquared1D::ProfileChiSquared1D() : IFittingAlgorithm() {}
+
+const std::string ProfileChiSquared1D::name() const {
+ return "ProfileChiSquared1D";
+}
+
+int ProfileChiSquared1D::version() const { return 1; }
+
+const std::string ProfileChiSquared1D::summary() const {
+ return "Profiles chi squared about its minimum to obtain parameter errors "
+ "for the input function.";
+}
+
+void ProfileChiSquared1D::initConcrete() {
+ declareProperty("Output", "", "A base name for output workspaces.");
+}
+
+void ProfileChiSquared1D::execConcrete() {
+ // Number of fiting parameters
+ auto nParams = m_function->nParams();
+ // Create an output table for displaying slices of the chi squared and
+ // the probabilitydensity function
+ auto pdfTable = API::WorkspaceFactory::Instance().createTable();
+
+ // Sigma confidence levels, could be an input but for now look for 1 sigma
+ // (68%) 2 sigma (95) and 3(99%) error bounds chi2 disturbution has 1 degree
+ // of freedom if we are changing 1 parameter at a time
+ boost::math::chi_squared chi2Dist(1);
+ std::array<double, 3> sigmas = {1, 2, 3};
+ std::array<double, 3> qvalues;
+ for (size_t i = 0; i < sigmas.size(); i++) {
+ double pvalue = std::erf(sigmas[i] / sqrt(2));
+ // find chi2 quanitile for given p value
+ qvalues[i] = boost::math::quantile(chi2Dist, pvalue);
+ }
+
+ // Find number of free parameter, should be >= 2
+ std::vector<int> freeParameters;
+ for (size_t ip = 0; ip < nParams; ++ip) {
+ if (m_function->isActive(ip)) {
+ freeParameters.push_back(static_cast<int>(ip));
+ }
+ }
+
+ if (freeParameters.size() < 2) {
+ throw std::invalid_argument("Function must have 2 or more free parameters");
+ }
+
+ std::string baseName = getProperty("Output");
+ Workspace_sptr ws = getProperty("InputWorkspace");
+ int workspaceIndex = getProperty("WorkspaceIndex");
+ MatrixWorkspace_sptr inputws = std::dynamic_pointer_cast<MatrixWorkspace>(ws);
+ if (baseName.empty()) {
+ baseName = "ProfileChiSquared1D";
+ }
+ declareProperty(
+ std::make_unique<API::WorkspaceProperty<API::ITableWorkspace>>(
+ "PDFs", "", Kernel::Direction::Output),
+ "The name of the TableWorkspace in which to store the "
+ "pdfs of fit parameters");
+ setPropertyValue("PDFs", baseName + "_pdf");
+ setProperty("PDFs", pdfTable);
+
+ // Create an output table for displaying the parameter errors.
+ auto errorsTable = API::WorkspaceFactory::Instance().createTable();
+ auto nameColumn = errorsTable->addColumn("str", "Parameter");
+ auto valueColumn = errorsTable->addColumn("double", "Value");
+ auto minValueColumn = errorsTable->addColumn("double", "Value at Min");
+ auto leftErrColumn = errorsTable->addColumn("double", "Left Error (1-sigma)");
+ auto rightErrColumn =
+ errorsTable->addColumn("double", "Right Error (1-sigma)");
+ auto leftErrColumn_2 =
+ errorsTable->addColumn("double", "Left Error (2-sigma)");
+ auto rightErrColumn_2 =
+ errorsTable->addColumn("double", "Right Error (2-sigma )");
+ auto leftErrColumn_3 =
+ errorsTable->addColumn("double", "Left Error (3-sigma)");
+ auto rightErrColumn_3 =
+ errorsTable->addColumn("double", "Right Error (3-sigma )");
+ auto quadraticErrColumn =
+ errorsTable->addColumn("double", "Quadratic Error (1-sigma)");
+ errorsTable->setRowCount(freeParameters.size());
+ declareProperty(
+ std::make_unique<API::WorkspaceProperty<API::ITableWorkspace>>(
+ "Errors", "", Kernel::Direction::Output),
+ "The name of the TableWorkspace in which to store the "
+ "values and errors of fit parameters");
+ setPropertyValue("Errors", baseName + "_errors");
+ setProperty("Errors", errorsTable);
+
+ // Calculate initial values
+ double chiSquared = 0.0;
+ double chiSquaredWeighted = 0.0;
+ double dof = 0;
+ API::FunctionDomain_sptr domain;
+ API::FunctionValues_sptr values;
+ m_domainCreator->createDomain(domain, values);
+ CalculateChiSquared::calcChiSquared(*m_function, nParams, *domain, *values,
+ chiSquared, chiSquaredWeighted, dof);
+ // Value of chi squared for current parameters in m_function
+ double chi0 = chiSquaredWeighted;
+
+ // Parameter bounds that define a volume in the parameter
+ // space within which the chi squared is being examined.
+ GSLVector lBounds(nParams);
+ GSLVector rBounds(nParams);
+
+ // Number of points in lines for plotting
+ size_t n = 100;
+ pdfTable->setRowCount(n);
+ const double fac = 1e-4;
+
+ for (auto p = 0u; p < freeParameters.size(); ++p) {
+ int row = p;
+ int ip = freeParameters[p];
+ // Add columns for the parameter to the pdf table.
+ auto parName = m_function->parameterName(ip);
+ nameColumn->read(row, parName);
+ // Parameter values
+ auto col1 = pdfTable->addColumn("double", parName);
+ col1->setPlotType(1);
+ // Chi squared values
+ auto col2 = pdfTable->addColumn("double", parName + "_chi2");
+ col2->setPlotType(2);
+ // PDF values
+ auto col3 = pdfTable->addColumn("double", parName + "_pdf");
+ col3->setPlotType(2);
+
+ double par0 = m_function->getParameter(ip);
+ double shift = fabs(par0 * fac);
+ if (shift == 0.0) {
+ shift = fac;
+ }
+
+ // Make a slice along this parameter
+ ChiSlice slice(m_function, ip, inputws, workspaceIndex, *domain, *values,
+ chi0, freeParameters);
+
+ // Find the bounds withn which the PDF is significantly above zero.
+ // The bounds are defined relative to par0:
+ // par0 + lBound is the lowest value of the parameter (lBound <= 0)
+ // par0 + rBound is the highest value of the parameter (rBound >= 0)
+ double lBound = slice.findBound(-shift);
+ double rBound = slice.findBound(shift);
+ lBounds[ip] = lBound;
+ rBounds[ip] = rBound;
+
+ // Approximate the slice with a polynomial.
+ // P is a vector of values of the polynomial at special points.
+ // A is a vector of Chebyshev expansion coefficients.
+ // The polynomial is defined on interval [lBound, rBound]
+ // The value of the polynomial at 0 == chi squared at par0
+ std::vector<double> P, A;
+ auto base = slice.makeApprox(lBound, rBound, P, A);
+
+ // Write n slice points into the output table.
+ double dp = (rBound - lBound) / static_cast<double>(n);
+ for (size_t i = 0; i < n; ++i) {
+ double par = lBound + dp * static_cast<double>(i);
+ double chi = base->eval(par, P);
+ col1->fromDouble(i, par0 + par);
+ col2->fromDouble(i, chi);
+ }
+
+ // Check if par0 is a minimum point of the chi squared
+ std::vector<double> AD;
+ // Calculate the derivative polynomial.
+ // AD are the Chebyshev expansion of the derivative.
+ base->derivative(A, AD);
+ // Find the roots of the derivative polynomial
+ std::vector<double> minima = base->roots(AD);
+ if (minima.empty()) {
+ minima.emplace_back(par0);
+ }
+
+ // If only 1 extremum is found assume (without checking) that it's a
+ // minimum.
+ // If there are more than 1, find the one with the smallest chi^2.
+ double chiMin = std::numeric_limits<double>::max();
+ double parMin = par0;
+ for (double minimum : minima) {
+ double value = base->eval(minimum, P);
+ if (value < chiMin) {
+ chiMin = value;
+ parMin = minimum;
+ }
+ }
+ // Get intersection of curve and line of constant q value to get confidence
+ // interval on parameter ip
+ valueColumn->fromDouble(row, par0);
+ minValueColumn->fromDouble(row, par0 + parMin);
+ for (size_t i = 0; i < qvalues.size(); i++) {
+ auto [rootsMin, rootsMax] =
+ getChiSquaredRoots(base, A, qvalues[i], rBound, lBound);
+ errorsTable->getColumn(3 + 2 * i)->fromDouble(row, rootsMin - parMin);
+ errorsTable->getColumn(4 + 2 * i)->fromDouble(row, rootsMax - parMin);
+ }
+
+ // Output the PDF
+ for (size_t i = 0; i < n; ++i) {
+ double chi = col2->toDouble(i);
+ col3->fromDouble(i, exp(-chi + chiMin));
+ }
+ // reset parameter values back to original value
+ m_function->setParameter(ip, par0);
+ }
+
+ // Square roots of the diagonals of the covariance matrix give
+ // the standard deviations in the quadratic approximation of the chi^2.
+ GSLMatrix V = getCovarianceMatrix();
+ for (size_t i = 0; i < freeParameters.size(); ++i) {
+ int ip = freeParameters[i];
+ quadraticErrColumn->fromDouble(i, sqrt(V.get(ip, ip)));
+ }
+}
+
+GSLMatrix ProfileChiSquared1D::getCovarianceMatrix() {
+ API::FunctionDomain_sptr domain;
+ API::FunctionValues_sptr values;
+ auto nParams = m_function->nParams();
+ m_domainCreator->createDomain(domain, values);
+ unfixParameters();
+ GSLJacobian J(*m_function, values->size());
+ m_function->functionDeriv(*domain, J);
+ refixParameters();
+ // Calculate the hessian at the current point.
+ GSLMatrix H;
+ H.resize(nParams, nParams);
+ for (size_t i = 0; i < nParams; ++i) {
+ for (size_t j = i; j < nParams; ++j) {
+ double h = 0.0;
+ for (size_t k = 0; k < values->size(); ++k) {
+ double w = values->getFitWeight(k);
+ h += J.get(k, i) * J.get(k, j) * w * w;
+ }
+ H.set(i, j, h);
+ if (i != j) {
+ H.set(j, i, h);
+ }
+ }
+ }
+ // Covariance matrix is inverse of hessian
+ GSLMatrix V(H);
+ V.invert();
+ return V;
+}
+
+/// Temporary unfix any fixed parameters.
+void ProfileChiSquared1D::unfixParameters() {
+ for (size_t i = 0; i < m_function->nParams(); ++i) {
+ if (!m_function->isActive(i)) {
+ m_function->unfix(i);
+ m_fixedParameters.emplace_back(i);
+ }
+ }
+}
+
+/// Restore the "fixed" status of previously unfixed paramters.
+void ProfileChiSquared1D::refixParameters() {
+ for (auto &fixedParameter : m_fixedParameters) {
+ m_function->fix(fixedParameter);
+ }
+ m_fixedParameters.clear();
+}
+
+std::tuple<double, double> ProfileChiSquared1D::getChiSquaredRoots(
+ const Functions::ChebfunBase_sptr &approximation,
+ std::vector<double> &coeffs, double qvalue, double rBound,
+ double lBound) const {
+ // Points of intersections with line chi^2 = 1 give an estimate of
+ // the standard deviation of this parameter if it's uncorrelated with the
+ // others.
+ // Cache original value of A0
+ auto Aold = coeffs[0];
+ // Now find roots of curve when quantile is subtracted
+ coeffs[0] = Aold - qvalue;
+ std::vector<double> roots = approximation->roots(coeffs);
+ std::sort(roots.begin(), roots.end());
+ if (roots.empty()) {
+ // Something went wrong; use the whole interval.
+ roots.resize(2);
+ roots[0] = lBound;
+ roots[1] = rBound;
+ } else if (roots.size() == 1) {
+ // Only one root found; use a bound for the other root.
+ if (roots.front() < 0) {
+ roots.emplace_back(rBound);
+ } else {
+ roots.insert(roots.begin(), lBound);
+ }
+ } else if (roots.size() > 2) {
+ // More than 2 roots; use the smallest and the biggest
+ auto smallest = roots.front();
+ auto biggest = roots.back();
+ roots.resize(2);
+ roots[0] = smallest;
+ roots[1] = biggest;
+ }
+ coeffs[0] = Aold;
+ return {roots[0], roots[1]};
+}
+
+} // namespace Algorithms
+} // namespace CurveFitting
+} // namespace Mantid
\ No newline at end of file
diff --git a/Framework/CurveFitting/test/Algorithms/CalculateChiSquaredTest.h b/Framework/CurveFitting/test/Algorithms/CalculateChiSquaredTest.h
index 4396cf4c4b745ff48bc932ecfab3b4c6e9e1df15..4c99b5eca7ed75c127045eaf1e52b9a3b18e5feb 100644
--- a/Framework/CurveFitting/test/Algorithms/CalculateChiSquaredTest.h
+++ b/Framework/CurveFitting/test/Algorithms/CalculateChiSquaredTest.h
@@ -169,99 +169,6 @@ public:
TS_ASSERT_DELTA(tester.chiSquared, 7151.0, 1.0);
}
- void test_errors_Thurber() {
- double x[] = {-3.067, -2.981, -2.921, -2.912, -2.840, -2.797, -2.702,
- -2.699, -2.633, -2.481, -2.363, -2.322, -1.501, -1.460,
- -1.274, -1.212, -1.100, -1.046, -0.915, -0.714, -0.566,
- -0.545, -0.400, -0.309, -0.109, -0.103, 0.010, 0.119,
- 0.377, 0.790, 0.963, 1.006, 1.115, 1.572, 1.841,
- 2.047, 2.200};
- double y[] = {80.574, 84.248, 87.264, 87.195, 89.076, 89.608,
- 89.868, 90.101, 92.405, 95.854, 100.696, 101.060,
- 401.672, 390.724, 567.534, 635.316, 733.054, 759.087,
- 894.206, 990.785, 1090.109, 1080.914, 1122.643, 1178.351,
- 1260.531, 1273.514, 1288.339, 1327.543, 1353.863, 1414.509,
- 1425.208, 1421.384, 1442.962, 1464.350, 1468.705, 1447.894,
- 1457.628};
- std::string fun("name=UserFunction,Formula=(b1 + b2*x + b3*x^2 + b4*x^3) / "
- "(1 + b5*x + b6*x^2 + b7*x^3),"
- "b1=1.2881396800E+03, b2=1.4910792535E+03, "
- "b3=5.8323836877E+02, b4=7.5416644291E+01, "
- "b5=9.6629502864E-01, b6=3.9797285797E-01, "
- "b7=4.9727297349E-02");
- double sigma[] = {4.6647963344E+00, 3.9571156086E+01, 2.8698696102E+01,
- 5.5675370270E+00, 3.1333340687E-02, 1.4984928198E-02,
- 6.5842344623E-03};
- size_t ndata = sizeof(x) / sizeof(double);
- Tester tester;
- tester.setTestCaseForErrorCalculations(ndata, x, y, fun);
- tester.runAlgorithm();
- tester.check1DSpectrum();
- tester.checkErrors(sigma);
- }
-
- void test_errors_Chwirut1() {
- double x[] = {
- 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
- 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.625, 0.625, 0.625, 0.625,
- 0.625, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75,
- 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.875, 0.875, 0.875,
- 0.875, 0.875, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1.25, 1.25, 1.25, 1.25, 1.25, 1.5, 1.5, 1.5, 1.5,
- 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.75, 1.75, 1.75, 1.75,
- 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75,
- 1.75, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
- 2, 2, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25,
- 2.25, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.75,
- 2.75, 2.75, 2.75, 2.75, 2.75, 2.75, 3, 3, 3, 3, 3,
- 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
- 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
- 3, 3, 3, 3.25, 3.25, 3.25, 3.25, 3.25, 3.75, 3.75, 3.75,
- 3.75, 3.75, 3.75, 3.75, 4, 4, 4, 4, 4, 4, 4.25,
- 4.25, 4.25, 4.25, 4.25, 4.75, 4.75, 4.75, 4.75, 4.75, 5, 5,
- 5, 5, 5, 5, 5.25, 5.25, 5.25, 5.25, 5.25, 5.75, 5.75,
- 5.75, 5.75, 5.75, 6, 6, 6, 6, 6, 6, 6, 6,
- 6, 6, 6, 6, 6};
- double y[] = {
- 76.8, 70.5, 78, 70.8, 90.6, 78, 74.1, 81.5,
- 81.3, 92.9, 75.8, 66.7, 81, 81.7, 78, 80,
- 80.7, 76.8, 79, 66, 76.9, 78.7, 67.3, 59.2,
- 60.9, 62, 60, 63.5, 59.9, 61.6, 59.5, 61,
- 63.8, 54.7, 71.6, 62, 61.3, 64.2, 62.4, 60.8,
- 64.5, 55.5, 63.6, 57.2, 58, 64.9, 48.5, 47.8,
- 48.8, 47.7, 53.2, 57.1, 40.8, 54, 47.5, 48,
- 50.3, 39.2, 42.5, 43.3, 41, 37.8, 29, 35.5,
- 35.2, 35.8, 32.9, 32, 32.5, 33.8, 39.8, 33.2,
- 30.7, 29.3625, 29.4, 28.4, 29.175, 25.5, 31.1, 26.85,
- 25.95, 28.9, 26.8, 28.95, 29.81, 29.3, 28.69, 29.8,
- 31.05, 21.67, 21, 20.32, 20, 22.57, 22.2, 25.7,
- 24, 29.8, 29.62, 25.99, 24.56, 21.15, 20.4, 23.6,
- 22.125, 21.4, 20.4, 23.775, 21.07, 20.2, 21, 21,
- 19.31, 16.95, 18.82, 16.3, 18.67, 17.7, 23.7, 23.81,
- 16.7625, 17.5125, 17.17, 16.65, 16.4625, 17.7375, 13.87, 15.64,
- 12.75, 13.95, 13.69, 15.64, 14.62, 15.19, 12.94, 13.12,
- 12.75, 13.31, 12.56, 13.35, 13.31, 11.81, 16.24, 13.87,
- 14.62, 17.7, 12.75, 13.84, 16.12, 12.07, 11.81, 18.07,
- 25.76, 15.56, 24.19, 13.12, 12.41, 13.2, 14.25, 12.525,
- 13.8, 9.67, 10.875, 9.45, 10.05, 10.39, 11.5875, 10.5375,
- 7.76, 8.17, 10.42, 11.25, 8.62, 8.74, 11.55, 9.15,
- 9.4125, 8.5875, 5.44, 8.175, 7.9125, 7.725, 7.125, 4.87,
- 7.35, 7.31, 9, 7.2, 7.87, 12.07, 8.55, 7.35,
- 6.1125, 4.01, 5.9625, 8.475, 5.9625, 3.75, 5.625, 6.1125,
- 8.025, 7.42, 6.67, 5.06, 4.87, 5.44, 6.64, 5.63,
- 5.44, 8.51, 8.74, 3.94, 5.63, 10.12};
- std::string fun("name=UserFunction,Formula=exp(-b1*x)/"
- "(b2+b3*x),b1=1.9027818370E-01,b2=6.1314004477E-03,b3=1."
- "0530908399E-02");
- double sigma[] = {2.1938557035E-02, 3.4500025051E-04, 7.9281847748E-04};
- size_t ndata = sizeof(x) / sizeof(double);
- Tester tester;
- tester.setTestCaseForErrorCalculations(ndata, x, y, fun);
- tester.runAlgorithm();
- tester.check1DSpectrum();
- tester.checkErrors(sigma);
- }
-
private:
class Tester {
// input parameters
@@ -281,7 +188,6 @@ private:
double EndX;
bool ignoreInvalidData;
std::string outputName;
- bool Weighted;
void makeXValues() {
size_t dlt = isHisto ? 1 : 0;
@@ -325,14 +231,12 @@ private:
double chiSquaredDividedByDOF;
double chiSquaredWeighted;
double chiSquaredWeightedDividedByDOF;
- ITableWorkspace_sptr errors;
- ITableWorkspace_sptr pdfs;
bool isExecuted;
Tester(size_t np = 3, size_t nd = 10, bool histo = true)
: nParams(np), nData(nd), isHisto(histo), xMin(-10), xMax(10),
workspaceIndex(0), StartX(EMPTY_DBL()), EndX(EMPTY_DBL()),
- ignoreInvalidData(false), Weighted(false),
+ ignoreInvalidData(false),
// output
chiSquared(-1), chiSquaredDividedByDOF(-1), chiSquaredWeighted(-1),
chiSquaredWeightedDividedByDOF(-1), isExecuted(false) {
@@ -353,10 +257,6 @@ private:
TS_ASSERT_THROWS_NOTHING(alg.setProperty("StartX", StartX));
TS_ASSERT_THROWS_NOTHING(alg.setProperty("EndX", EndX));
}
- if (!outputName.empty()) {
- TS_ASSERT_THROWS_NOTHING(alg.setProperty("Output", outputName));
- TS_ASSERT_THROWS_NOTHING(alg.setProperty("Weighted", Weighted));
- }
TS_ASSERT_THROWS_NOTHING(alg.execute());
isExecuted = alg.isExecuted();
if (isExecuted) {
@@ -365,12 +265,6 @@ private:
chiSquaredWeighted = alg.getProperty("ChiSquaredWeighted");
chiSquaredWeightedDividedByDOF =
alg.getProperty("ChiSquaredWeightedDividedByDOF");
- if (!outputName.empty()) {
- errors = AnalysisDataService::Instance().retrieveWS<ITableWorkspace>(
- "out_errors");
- pdfs = AnalysisDataService::Instance().retrieveWS<ITableWorkspace>(
- "out_pdf");
- }
}
}
@@ -482,8 +376,6 @@ private:
FunctionValues y(x);
function->function(x, y);
double tmp = yValues[i] - y[0];
- // std::cerr << "test " << xValue << ' ' << yValues[i] << ' ' << y[0]
- // << '\n';
sum2 += tmp * tmp;
tmp /= eValues[i];
sum2w += tmp * tmp;
@@ -500,20 +392,6 @@ private:
TS_ASSERT_DELTA(sum2, chiSquaredDividedByDOF, 1e-10);
TS_ASSERT_DELTA(sum2w, chiSquaredWeightedDividedByDOF, 1e-10);
}
-
void checkFailed() { TS_ASSERT(!isExecuted); }
-
- void checkErrors(double *sigma) {
- size_t np = function->nParams();
- TS_ASSERT_EQUALS(np, errors->rowCount());
- for (size_t i = 0; i < np; ++i) {
- TS_ASSERT_LESS_THAN(fabs(errors->Double(i, 6)), 1e-4);
- TS_ASSERT_DELTA(errors->Double(i, 1) / errors->Double(i, 2), 1.0, 1e-5);
- TS_ASSERT_DELTA(0.5 * (errors->Double(i, 4) - errors->Double(i, 3)) /
- errors->Double(i, 5),
- 1.0, 2.0);
- TS_ASSERT_DELTA(errors->Double(i, 5) / sigma[i], 1.0, 0.01);
- }
- }
};
};
diff --git a/Framework/CurveFitting/test/Algorithms/ProfileChiSquared1DTest.h b/Framework/CurveFitting/test/Algorithms/ProfileChiSquared1DTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..2580ca461128bf9d0afa9078c7c66ebb4dab5a80
--- /dev/null
+++ b/Framework/CurveFitting/test/Algorithms/ProfileChiSquared1DTest.h
@@ -0,0 +1,169 @@
+// Mantid Repository : https://github.com/mantidproject/mantid
+//
+// Copyright © 2021 ISIS Rutherford Appleton Laboratory UKRI,
+// NScD Oak Ridge National Laboratory, European Spallation Source,
+// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
+// SPDX - License - Identifier: GPL - 3.0 +
+#pragma once
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidAPI/AnalysisDataService.h"
+#include "MantidCurveFitting/Algorithms/ProfileChiSquared1D.h"
+#include "MantidDataObjects/TableWorkspace.h"
+
+using Mantid::CurveFitting::Algorithms::ProfileChiSquared1D;
+using namespace Mantid;
+using namespace Mantid::API;
+using namespace Mantid::DataObjects;
+
+namespace {
+constexpr char *linearFunctionString =
+ "name = LinearBackground, A0=0.8753627851076761, A1 = "
+ "2.026706319695708 ";
+}
+
+class ProfileChiSquared1DTest : public CxxTest::TestSuite {
+public:
+ static ProfileChiSquared1DTest *createSuite() {
+ return new ProfileChiSquared1DTest();
+ }
+ static void destroySuite(ProfileChiSquared1DTest *suite) {
+ AnalysisDataService::Instance().clear();
+ delete suite;
+ }
+
+ void loadLinearData(const std::string &workspaceName) {
+ auto algo = AlgorithmManager::Instance().create("Load");
+ algo->setPropertyValue("Filename", "ProfileChiSquared1DData_linear.nxs");
+ algo->setPropertyValue("OutputWorkspace", workspaceName);
+ algo->execute();
+ }
+
+ void executeAlgorithmOnLinearData(const std::string &outputName) {
+ std::string wsName = "ProfileChiSquared1DData_linear";
+ loadLinearData(wsName);
+ auto ws = AnalysisDataService::Instance().retrieveWS<Workspace>(wsName);
+ auto functionString = std::string(linearFunctionString);
+ auto profileAlg = ProfileChiSquared1D();
+ profileAlg.initialize();
+ profileAlg.setProperty("Function", functionString);
+ profileAlg.setProperty("InputWorkspace", ws);
+ profileAlg.setProperty("Output", outputName);
+ profileAlg.execute();
+ }
+
+ void test_Init() {
+ ProfileChiSquared1D alg;
+ TS_ASSERT_THROWS_NOTHING(alg.initialize())
+ TS_ASSERT(alg.isInitialized())
+ }
+
+ void test_Alg_produces_expected_outputs() {
+ executeAlgorithmOnLinearData("OutputName1");
+ TS_ASSERT(AnalysisDataService::Instance().doesExist("OutputName1_errors"));
+ TS_ASSERT(AnalysisDataService::Instance().doesExist("OutputName1_pdf"));
+
+ // if name is empty, workspaces will use default name of ProfileChiSquared1D
+ executeAlgorithmOnLinearData("");
+ TS_ASSERT(AnalysisDataService::Instance().doesExist(
+ "ProfileChiSquared1D_errors"));
+ TS_ASSERT(
+ AnalysisDataService::Instance().doesExist("ProfileChiSquared1D_pdf"));
+
+ AnalysisDataService::Instance().clear();
+ }
+ // the tests for these linear problems are comparing against analytical
+ // calculations which can be solved in closed form for a linear function
+ void test_errors_for_linear_function_are_correct() {
+ executeAlgorithmOnLinearData("OutputName1");
+ TableWorkspace_sptr errorsTable;
+ TS_ASSERT_THROWS_NOTHING(
+ errorsTable =
+ AnalysisDataService::Instance().retrieveWS<TableWorkspace>(
+ "OutputName1_errors"));
+
+ TS_ASSERT_EQUALS(errorsTable->String(0, 0), "A0");
+ TS_ASSERT_EQUALS(errorsTable->String(1, 0), "A1");
+ // check a0 parameter
+ TS_ASSERT_DELTA(errorsTable->Double(0, 1), 0.8753627851076761, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(0, 2), 0.8753627851076761, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(0, 3), -0.03447782006595415, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(0, 4), 0.03447782006595401, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(0, 5), -0.06895564013190685, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(0, 6), 0.06895564013190683, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(0, 7), -0.10343346019785989, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(0, 8), 0.1034334601978597, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(0, 9), 0.03447782006595295, 1e-6);
+ // check a1 parameter
+ TS_ASSERT_DELTA(errorsTable->Double(1, 1), 2.026706319695708, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(1, 2), 2.026706319695708, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(1, 3), -0.006137378377995283, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(1, 4), 0.006137378377995297, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(1, 5), -0.012274756755989097, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(1, 6), 0.012274756755989113, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(1, 7), -0.01841213513398322, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(1, 8), 0.01841213513398322, 1e-6);
+ TS_ASSERT_DELTA(errorsTable->Double(1, 9), 0.006137378377994362, 1e-6);
+ AnalysisDataService::Instance().clear();
+ }
+
+ void test_pdf_values_for_linear_function_are_correct() {
+ executeAlgorithmOnLinearData("OutputName2");
+ TableWorkspace_sptr pdfTable;
+ TS_ASSERT_THROWS_NOTHING(
+ pdfTable = AnalysisDataService::Instance().retrieveWS<TableWorkspace>(
+ "OutputName2_pdf"));
+ // check some values of pdf table
+ TS_ASSERT_DELTA(pdfTable->Double(0, 0), 0.696088486717624, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(0, 1), 27.03687327118116, 1e-2);
+ TS_ASSERT_DELTA(pdfTable->Double(0, 3), 2.0007644788036028, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(0, 4), 17.86634783053958, 1e-2);
+
+ TS_ASSERT_DELTA(pdfTable->Double(25, 0), 0.78572563591265, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(25, 1), 6.75921831779535, 1e-2);
+ TS_ASSERT_DELTA(pdfTable->Double(25, 3), 2.0137353992496556, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(25, 4), 4.466586957634646, 1e-2);
+
+ TS_ASSERT_DELTA(pdfTable->Double(50, 0), 0.8753627851076761, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(50, 1), -5.32907051820075e-14, 1e-2);
+ TS_ASSERT_DELTA(pdfTable->Double(50, 3), 2.026706319695708, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(50, 4), -3.197442310920451e-13, 1e-2);
+
+ TS_ASSERT_DELTA(pdfTable->Double(75, 0), 0.9649999343027021, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(75, 1), 6.759218317794955, 1e-2);
+ TS_ASSERT_DELTA(pdfTable->Double(75, 3), 2.0396772401417604, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(75, 4), 4.4665869576346875, 1e-2);
+
+ TS_ASSERT_DELTA(pdfTable->Double(99, 0), 1.051051597529927, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(99, 1), 25.96621308964161, 1e-2);
+ TS_ASSERT_DELTA(pdfTable->Double(99, 3), 2.052129323769971, 1e-6);
+ TS_ASSERT_DELTA(pdfTable->Double(99, 4), 17.15884045645028, 1e-2);
+
+ AnalysisDataService::Instance().clear();
+ }
+
+ void test_errors_table_has_correct_shape() {
+ executeAlgorithmOnLinearData("OutputName3");
+ TableWorkspace_sptr errorsTable;
+ TS_ASSERT_THROWS_NOTHING(
+ errorsTable =
+ AnalysisDataService::Instance().retrieveWS<TableWorkspace>(
+ "OutputName3_errors"));
+ TS_ASSERT_EQUALS(errorsTable->columnCount(), 10);
+ TS_ASSERT_EQUALS(errorsTable->rowCount(), 2);
+ AnalysisDataService::Instance().clear();
+ }
+
+ void test_pdf_table_has_correct_shape() {
+ executeAlgorithmOnLinearData("OutputName4");
+ TableWorkspace_sptr pdfTable;
+ TS_ASSERT_THROWS_NOTHING(
+ pdfTable = AnalysisDataService::Instance().retrieveWS<TableWorkspace>(
+ "OutputName4_pdf"));
+ TS_ASSERT_EQUALS(pdfTable->columnCount(), 6);
+ TS_ASSERT_EQUALS(pdfTable->rowCount(), 100);
+ AnalysisDataService::Instance().clear();
+ };
+};
\ No newline at end of file
diff --git a/Testing/Data/UnitTest/ProfileChiSquared1DData_linear.nxs.md5 b/Testing/Data/UnitTest/ProfileChiSquared1DData_linear.nxs.md5
new file mode 100644
index 0000000000000000000000000000000000000000..c64a34c81195c515f1767b423cb299fb5df184fe
--- /dev/null
+++ b/Testing/Data/UnitTest/ProfileChiSquared1DData_linear.nxs.md5
@@ -0,0 +1 @@
+d3526ef8e18fc293ac5c81a5208722bd
diff --git a/docs/source/algorithms/CalculateChiSquared-v1.rst b/docs/source/algorithms/CalculateChiSquared-v1.rst
index e96c3a6b10ef2af55cd50f433b4d741f5c7cd61a..4e038fab7ccf3008d7abbdb24cd1165f1d5da788 100644
--- a/docs/source/algorithms/CalculateChiSquared-v1.rst
+++ b/docs/source/algorithms/CalculateChiSquared-v1.rst
@@ -33,35 +33,6 @@ Finally, ChiSquaredWeightedDividedByDOF is
:math:`\chi_{4}^{2} = \chi_{3}^{2} / DOF`
-Parameter errors
-################
-
-Setting the Output property to a non-empty string makes the algorithm explore the surface of the :math:`\chi^{2}`
-around its minimum and estimate the standard deviations for the parameters. The value of the property is a base name
-for two output table workspaces: '<Output>_errors' and '<Output>_pdf'. The former workspace contains parameter error
-estimates and the latter shows :math:`\chi^{2}`'s 1d slices along each parameter (keeping all other fixed).
-
-The error table has the following columns:
-
-=============== ===========
-Column Description
-=============== ===========
-Parameter Parameter name
-Value Parameter value passed with the Function property
-Value at Min The minimum point of the 1d slice of the :math:`\chi^{2}`. If the Function is at the minimum then
- Value at Min should be equal to Value.
-Left Error The negative deviation from the minimum point equivalent to :math:`1\sigma`. Estimated from analisys
- of the surface.
-Right Error The positive deviation from the minimum point equivalent to :math:`1\sigma`. Estimated from analisys
- of the surface.
-Quadratic Error :math:`1\sigma` standard deviation in the quadratic approximation of the surface.
-Chi2 Min The value of :math:`\chi^{2}` at the minimum relative to the test point.
-=============== ===========
-
-The pdf table contains slices of the :math:`\chi^{2}` along each parameter. It has 3 column per parameter. The first column of the 3
-is the parameter values, the second has the :math:`\chi^{2}` and the third is the probability density function normalised to
-have 1 at the maximum.
-
Usage
-----
**Example 1**
@@ -119,38 +90,6 @@ Output:
Chi squared weighted / DOF is 0.0
Chi squared weighted / NDATA is 0.0
-**Example 2**
-
-.. testcode::
-
- import numpy as np
- # Create a workspace and fill it with some gaussian data and some noise
- n = 100
- x = np.linspace(-10,10,n)
- y = np.exp(-x*x/2) + np.random.normal(0.0, 0.01, n)
- e = [1] * n
- ws = CreateWorkspace(x,y,e)
-
- # Gefine a Gaussian with exactly the same parameters that were used to
- # generate the data
- fun_t = 'name=Gaussian,Height=%s,PeakCentre=%s,Sigma=%s'
- fun = fun_t % (1, 0, 1)
- # Test the chi squared.
- CalculateChiSquared(fun,ws,Output='Test0')
- # Check the Test0_errors table and see that the parameters are not at minimum
-
- # Fit the function
- res = Fit(fun,ws,Output='out')
- # res[3] is a table with the fitted parameters
- nParams = res[3].rowCount() - 1
- params = [res[3].cell(i,1) for i in range(nParams)]
- # Build a new function and populate it with the fitted parameters
- fun = fun_t % tuple(params)
- # Test the chi squared.
- CalculateChiSquared(fun,ws,Output='Test1')
- # Check the Test1_errors table and see that the parameters are at minimum now
-
-
.. categories::
.. sourcelink::
diff --git a/docs/source/algorithms/ProfileChiSquared1D-v1.rst b/docs/source/algorithms/ProfileChiSquared1D-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..f2bab9e2a4bb46adf07e93edde56e7c03a143ff1
--- /dev/null
+++ b/docs/source/algorithms/ProfileChiSquared1D-v1.rst
@@ -0,0 +1,184 @@
+
+.. algorithm::
+
+.. summary::
+
+.. relatedalgorithms::
+
+.. properties::
+
+Description
+-----------
+
+This algorithm explores the surface of the :math:`\chi^{2}` around its minimum and estimates the standard deviations for the parameters.
+The value of the output property is a base name for two output table workspaces: '<Output>_errors' and '<Output>_pdf'.
+The former workspace contains parameter error estimates and the latter shows :math:`\chi^{2}`'s 1d slices along each parameter.
+
+The procedure for calculating errors of parameters is described in Chapter 15 of Numerical recipes in C [1] and Chapter 9
+of Statistical Data Analysis [2]. Here, we summarise the main results.
+
+Consider the input dataset :math:`D_0`, with fit parameters :math:`\mathbf a_0`. Assuming Gaussian noise in the data, it is possible
+to create :math:`n` artificial datasets :math:`D_j` for :math:`j=1,2,..,n`. If we were to run the fit on each dataset,
+a set of parameters :math:`\mathbf a_j` would be obtained. The distribution of these parameters,
+about our original parameter set :math:`\delta \mathbf a = \mathbf a_j - \mathbf a_0` is described by the multivariate normal distribution,
+
+.. math::
+ P(\delta \mathbf a ) \propto exp (-\Delta \chi^2)
+
+where :math:`\Delta \chi^2=\chi^2(\mathbf a_J) - \chi^2(\mathbf a_0)`, is the difference between the chi squared statistics obtained for each parameter set,
+
+.. math::
+ \Delta \chi^2 = \sum_{i \in D_0} \left ( \frac{y_i -f(x_i;\mathbf a_J)}{\sigma_i}\right)^2 - \left ( \frac{y_i -f(x_i;\mathbf a_0)}{\sigma_i}\right)^2
+
+If we consider the variation of a single parameter :math:`a_k`, while the other parameters are the values that minimize :math:`\chi^2`,
+the quantity :math:`\Delta \chi^2` will be distributed as a chi squared distributed, :math:`f_{\chi^2}(x; \nu)` with 1 degree of freedom,
+:math:`\nu=1`. From this distribution, the probability of finding increases less than :math:`\Delta \chi^2` can be found from the integral,
+
+.. math::
+ p = P(X < \Delta \chi^2 ) = \int_0^{\Delta \chi^2} f_{\chi^2}(x; 1) dx
+
+where :math:`p` is the desired confidence level. For instance, if we consider a confidence level of :math:`p=68.3\%` (1-sigma),
+we find a critical value for the chi squared fluctuation of :math:`\Delta \chi^2 < 1`. From this desired confidence level,
+we can then calculate the parameter variation, (:math:`\sigma_{left}`, :math:`\sigma_{right}`), that increases chi squared by 1
+and obtain a :math:`68.3\%` confidence interval for the single parameter :math:`a_k`:
+
+.. math::
+ P( a_k - \sigma_{left} < a_k < a_k + \sigma_{right}) = 68.3\%
+
+This algorithm obtains the 1-sigma confidence level by varying a single parameter at a time and recording the extremums
+of the values that increase :math:`\chi^2` by :math:`1`. The results are outputted in the table '<Output>_errors',
+which reports the left and right 1-sigma errors. This procedure is repeated for 2-sigma and 3-sigma,
+with the results displayed in the columns (2-sigma) and (3-sigma). An additional value is also reported, termed the quadratic error.
+This is the 1-sigma error obtained from a Taylor expansion of the chi squared function about its minimum:
+
+.. math::
+ \chi^2(\mathbf a + \delta \mathbf a) = \chi^2(\mathbf a)_{min} + \delta \mathbf a^T \mathbf C \delta \mathbf a
+
+where :math:`\mathbf{C}` is the covariance matrix obtained from least squares fitting of the data, see the :ref:`Fitting <algm-Fit>` documentation for further details.
+For a linear model (or approximately linear model) the quadratic error should be equal to the left and right errors, i.e it will be symmetric.
+This can be seen in the following example
+
+.. code-block:: python
+
+ from mantid.simpleapi import *
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+ # create the data
+ x = np.linspace(1, 10, 250)
+ np.random.seed(0)
+ y = 1 + 2.0*x + 0.4*np.random.randn(x.size)
+ ws = CreateWorkspace(x, y)
+
+ # run the fit
+ func = "name=LinearBackground, A0=1, A1=2";
+ func_ = FunctionFactory.Instance().createInitialized(func)
+
+ fit_output = Fit(InputWorkspace=ws, Function=func_, Output="LinearFit")
+ a0_fit = fit_output.OutputParameters.column(1)[0]
+ a1_fit = fit_output.OutputParameters.column(1)[1]
+
+ # explore the chi squared profile for the fit parameters
+ fitted_func = "name=LinearBackground, A0={}, A1={}".format(a0_fit, a1_fit);
+ ProfileChiSquared1D(fitted_func, ws, Output="LinearProfile")
+
+ # print left and right errors of parameters
+ # you should note that they are approx equal to the quadratic error for this linear model
+ error_table = mtd["LinearProfile_errors"]
+ lerror_a0 = error_table.column(3)[0]
+ rerror_a0= error_table.column(4)[0]
+ qerror_a0 = error_table.column(9)[0]
+ print("1-sigma error bounds of A0 are {} and {}, with quadratic estimate {}".format(lerror_a0, rerror_a0, qerror_a0))
+
+ lerror_a1 = error_table.column(3)[1]
+ rerror_a1= error_table.column(4)[1]
+ qerror_a1 = error_table.column(9)[1]
+ print("1-sigma error bounds of A1 are {} and {}, with quadratic estimate {}".format(lerror_a1, rerror_a1, qerror_a1))
+
+
+For a non-linear model, it's possible that the left and right variances will not be equal, leading to an asymmetric error.
+This is shown in the example below:
+
+.. code-block:: python
+
+ # import mantid algorithms, numpy and matplotlib
+ from mantid.simpleapi import *
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+
+ # create decaying exponential data
+ x = np.linspace(1, 10, 250)
+ np.random.seed(0)
+ y = 3.0*np.exp(-x/2) + 0.1*np.random.randn(x.size)
+ ws = CreateWorkspace(x, y)
+
+ # run the fit
+ func = "name=ExpDecay,Height=3.0, Lifetime=0.5";
+ func_ = FunctionFactory.Instance().createInitialized(func)
+
+ fit_output = Fit(InputWorkspace=ws, Function=func_, Output="ExpFit")
+ height_fit = fit_output.OutputParameters.column(1)[0]
+ lifetime_fit = fit_output.OutputParameters.column(1)[1]
+
+ # explore the chi squared profile for the fit parameters
+ fitted_func = "name=ExpDecay, Height={}, Lifetime={}".format(height_fit, lifetime_fit);
+ ProfileChiSquared1D(fitted_func, ws, Output="ExpProfile")
+
+ # print left and right errors of parameters
+ # you should note that they differ from the quadratic errors
+ error_table = mtd["ExpProfile_errors"]
+ lerror_height = error_table.column(3)[0]
+ rerror_height= error_table.column(4)[0]
+ qerror_height = error_table.column(9)[0]
+ print("1-sigma error bounds of Height are {} and {}, with quadratic estimate {}".format(lerror_height, rerror_height, qerror_height))
+
+ lerror_lifetime = error_table.column(3)[1]
+ rerror_lifetime= error_table.column(4)[1]
+ qerror_lifetime = error_table.column(9)[1]
+ print("1-sigma error bounds of Lifetime are {} and {}, with quadratic estimate {}".format(lerror_lifetime, rerror_lifetime, qerror_lifetime))
+
+
+For each problem, the error table has the following columns:
+
+====================== ==============
+Column Description
+====================== ==============
+Parameter Parameter name
+Value Parameter value passed with the Function property
+Value at Min The minimum point of the 1d slice of the :math:`\chi^{2}`. If the Function is at the minimum then
+ Value at Min should be equal to Value.
+Left Error (1-sigma) The negative deviation from the minimum point equivalent to :math:`1\sigma`. Estimated from analisys
+ of the surface.
+Right Error (1-sigma) The positive deviation from the minimum point equivalent to :math:`1\sigma`. Estimated from analisys
+ of the surface.
+Left Error (2-sigma) The negative deviation from the minimum point equivalent to :math:`2\sigma`. Estimated from analisys
+ of the surface.
+Right Error (2-sigma) The positive deviation from the minimum point equivalent to :math:`2\sigma`. Estimated from analisys
+ of the surface.
+Left Error (3-sigma) The negative deviation from the minimum point equivalent to :math:`3\sigma`. Estimated from analisys
+ of the surface.
+Right Error (3-sigma) The positive deviation from the minimum point equivalent to :math:`3\sigma`. Estimated from analisys
+ of the surface.
+Quadratic Error :math:`1\sigma` standard deviation in the quadratic approximation of the surface.
+====================== ==============
+
+This algorithm also reports a probability density function (PDF) for each parameter, stored in the '<Output>_pdf' table.
+This PDF is found from Eq. (1), which relates the increase in chi squared, to the probability of the parameter variation.
+The pdf table also contains slices of the :math:`\chi^{2}` along each parameter. It has 3 column per parameter. The first column of the 3
+is the parameter values, the second has the :math:`\chi^{2}` and the third is the probability density function normalised to
+have 1 at the maximum. Plotting the second column of each parameter will show the change in :math:`\chi^{2}` with respect to
+the parameter value.
+
+References
+----------
+
+[1] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. 1992.
+Numerical recipes in C (2nd ed.): the art of scientific computing. Cambridge University Press, USA.
+
+[2] G. Cowan, Statistical Data Analysis, Clarendon, Oxford, 1998
+
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/concepts/Fitting.rst b/docs/source/concepts/Fitting.rst
index a7384ee07aabaa01350386bb9f391ae8078f0c43..0f918a7180eb7a20fbeebf7dcae7e8126d594afa 100644
--- a/docs/source/concepts/Fitting.rst
+++ b/docs/source/concepts/Fitting.rst
@@ -132,7 +132,7 @@ is :ref:`FitPeak <algm-FitPeak>` which aims to fit single peaks
with some checks to ensure results are physical.
If you remain unsure if a given fit was successful then e.g. try the
-tool :ref:`CalculateChiSquared <algm-CalculateChiSquared>`, which allows
+tool :ref:`ProfileChiSquared1D <algm-ProfileChiSquared1D>`, which allows
inspection of the cost function in the neighbourhood of a found minimum.
diff --git a/docs/source/release/v6.1.0/mantidworkbench.rst b/docs/source/release/v6.1.0/mantidworkbench.rst
index e33d38dfcbea5df9ba402e5b291087b365adb3b1..48012792269b866fb0825dd9a04cda8a50e440e8 100644
--- a/docs/source/release/v6.1.0/mantidworkbench.rst
+++ b/docs/source/release/v6.1.0/mantidworkbench.rst
@@ -18,6 +18,10 @@ Added Floating/On Top setting for all the windows that are opened by workbench (
- A new empty facility with empty instrument is the default facility now, and
user has to select their choice of facility (including ISIS) and instrument for the first time
+
+- Added an algorithm ProfileChiSquared1D to profile chi squared after a fit. This can be used
+ to find better estimates of parameter errors.
+
- Instrument view: when in tube selection mode, the sum of pixel counts is now output to the selection pane.
Bugfixes