Skip to content

Commit

Permalink
Re #11617. Added weighted sums.
Browse files Browse the repository at this point in the history
  • Loading branch information
@mantid-roman
mantid-roman committed May 22, 2015
1 parent 25c34a2 commit 28cb000
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 25 deletions.
39 changes: 27 additions & 12 deletions Code/Mantid/Framework/CurveFitting/src/CalculateChiSquared.cpp
View file
Expand Up @@ -42,6 +42,11 @@ void CalculateChiSquared::initConcrete() {
"Output value of chi squared divided by the "
"number of degrees of freedom (NofData "
"- nOfParams).", Direction::Output);
declareProperty("ChiSquaredWeighted", 0.0, "Output value of weighted chi squared.", Direction::Output);
declareProperty("ChiSquaredWeightedDividedByDOF", 0.0,
"Output value of weighted chi squared divided by the "
"number of degrees of freedom (NofData "
"- nOfParams).", Direction::Output);
}

//----------------------------------------------------------------------------------------------
Expand All @@ -61,39 +66,49 @@ void CalculateChiSquared::execConcrete() {
// Calculate function values.
m_function->function(*domain, *values);

// Get the number of free fitting parameters
size_t nParams = 0;
for (size_t i = 0; i < m_function->nParams(); ++i) {
if (!m_function->isFixed(i))
nParams += 1;
}
double dof = - static_cast<double>(nParams);

// Calculate the chi squared.
double chiSquared = 0.0;
double chiSquaredWeighted = 0.0;
for(size_t i = 0; i < values->size(); ++i) {
if (values->getFitWeight(i) > 0.0) {
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;
}
}
g_log.debug() << "Chi squared " << chiSquared << std::endl;
g_log.notice() << "Chi squared " << chiSquared << std::endl;
g_log.notice() << "Chi squared weighted " << chiSquaredWeighted << std::endl;

// Store the result.
setProperty("ChiSquared", chiSquared);
setProperty("chiSquaredWeighted", chiSquaredWeighted);

// Divide by the DOF
// Get the number of free fitting parameters
size_t nParams = 0;
for (size_t i = 0; i < m_function->nParams(); ++i) {
if (!m_function->isFixed(i))
nParams += 1;
}
double dof =
static_cast<double>(domain->size()) - static_cast<double>(nParams);
g_log.debug() << "DOF " << dof << std::endl;
if (dof <= 0.0) {
dof = 1.0;
g_log.warning() << "DOF has a non-positive value, changing to 1,0."
g_log.warning() << "DOF has a non-positive value, changing to 1.0."
<< std::endl;
}
chiSquared /= dof;
g_log.debug() << "Chi squared / DOF " << chiSquared << std::endl;
chiSquaredWeighted /= dof;
g_log.notice() << "Chi squared / DOF " << chiSquared << std::endl;
g_log.notice() << "Chi squared weighed / DOF " << chiSquaredWeighted << std::endl;

// Store the result.
setProperty("ChiSquaredDividedByDOF", chiSquared);
setProperty("ChiSquaredWeightedDividedByDOF", chiSquaredWeighted);
}

} // namespace CurveFitting
Expand Down
31 changes: 18 additions & 13 deletions Code/Mantid/Framework/CurveFitting/test/CalculateChiSquaredTest.h
View file
Expand Up @@ -129,22 +129,13 @@ class CalculateChiSquaredTest : public CxxTest::TestSuite {
TS_ASSERT_DELTA(tester.chiSquared, 1450.39, 1.0);
}

void test_1D_values_divide_by_dof() {
Tester tester;
tester.set1DFunction();
tester.set1DSpectrumValues();
tester.runAlgorithm();
tester.check1DSpectrum();
TS_ASSERT_DELTA(tester.chiSquared, 236.5, 1.0);
}

void test_1D_values_divide_by_dof_fixed_params() {
Tester tester(1);
tester.set1DFunction();
tester.set1DSpectrumValues();
tester.runAlgorithm();
tester.check1DSpectrum();
TS_ASSERT_DELTA(tester.chiSquared, 184.0, 1.0);
TS_ASSERT_DELTA(tester.chiSquared, 1655.0, 1.0);
}

void test_1D_values_divide_by_dof_zero() {
Expand Down Expand Up @@ -223,12 +214,16 @@ class CalculateChiSquaredTest : public CxxTest::TestSuite {
// algorithm output
double chiSquared;
double chiSquaredDividedByDOF;
double chiSquaredWeighted;
double chiSquaredWeightedDividedByDOF;
bool isExecuted;

Tester(size_t np = 3, size_t nd = 10, bool histo = true)
: nParams(np), nData(nd), isHisto(histo), workspaceIndex(0), xMin(-10),
xMax(10), StartX(EMPTY_DBL()), EndX(EMPTY_DBL()),
ignoreInvalidData(false), chiSquared(-1), chiSquaredDividedByDOF(-1), isExecuted(false) {
ignoreInvalidData(false), chiSquared(-1), chiSquaredDividedByDOF(-1),
chiSquaredWeighted(-1), chiSquaredWeightedDividedByDOF(-1),
isExecuted(false) {
makeXValues();
}

Expand All @@ -249,6 +244,8 @@ class CalculateChiSquaredTest : public CxxTest::TestSuite {
if (isExecuted) {
chiSquared = alg.getProperty("ChiSquared");
chiSquaredDividedByDOF = alg.getProperty("ChiSquaredDividedByDOF");
chiSquaredWeighted = alg.getProperty("ChiSquaredWeighted");
chiSquaredWeightedDividedByDOF = alg.getProperty("ChiSquaredWeightedDividedByDOF");
}
}

Expand Down Expand Up @@ -296,6 +293,7 @@ class CalculateChiSquaredTest : public CxxTest::TestSuite {
auto space = WorkspaceFactory::Instance().create("Workspace2D", nSpec,
nData + dn, nData);
space->dataX(0).assign(xBins.begin(), xBins.end());
space->dataE(0).assign(nData, 10.0);
workspace = space;
}

Expand Down Expand Up @@ -327,25 +325,32 @@ class CalculateChiSquaredTest : public CxxTest::TestSuite {
TS_ASSERT(isExecuted);
setDefaultXRange();
double sum2 = 0.0;
double sum2w = 0.0;
auto &yValues = dynamic_cast<MatrixWorkspace&>(*workspace).readY(workspaceIndex);
auto &eValues = dynamic_cast<MatrixWorkspace&>(*workspace).readE(workspaceIndex);
double dof = - double(nParams);
for (size_t i = 0; i < xValues.size(); ++i) {
const double xValue = xValues[i];
if (xValue >= StartX && xValue <= EndX && isGoodValue(yValues[i],eValues[i])) {
FunctionDomain1DVector x(xValue);
FunctionValues y(x);
function->function(x, y);
const double tmp = yValues[i] - y[0];
double tmp = yValues[i] - y[0];
//std::cerr << "test " << xValue << ' ' << yValues[i] << ' ' << y[0] << std::endl;
sum2 += tmp * tmp;
tmp /= eValues[i];
sum2w += tmp * tmp;
dof += 1.0;
}
}
TS_ASSERT_DIFFERS(sum2, 0);
TS_ASSERT_DELTA(sum2, chiSquared, 1e-10);
double dof = double(nData) - double(nParams);
TS_ASSERT_DELTA(sum2w, chiSquaredWeighted, 1e-10);
if (dof <= 0.0) dof = 1.0;
sum2 /= dof;
sum2w /= dof;
TS_ASSERT_DELTA(sum2, chiSquaredDividedByDOF, 1e-10);
TS_ASSERT_DELTA(sum2w, chiSquaredWeightedDividedByDOF, 1e-10);
}

void checkFailed() {
Expand Down

0 comments on commit 28cb000

Please sign in to comment.