Commit 6ed67d3c authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Fixed 'chi2' calculation for M-estimators

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@134 281f6f2b-e318-4fd1-8bce-1a4ba7aab212
parent 438138d3
......@@ -7,7 +7,7 @@ PROJECT(GBL)
# project version
SET( ${PROJECT_NAME}_VERSION_MAJOR 2 )
SET( ${PROJECT_NAME}_VERSION_MINOR 1 )
SET( ${PROJECT_NAME}_VERSION_PATCH 0 )
SET( ${PROJECT_NAME}_VERSION_PATCH 1 )
# make life easier and simply use the ilcsoft default settings
# load default ilcsoft settings (install prefix, build type, rpath, etc.)
......
......@@ -128,6 +128,7 @@ private:
double thePrecision; ///< Precision (1/sigma**2)
unsigned int theTrajectory; ///< Trajectory number
unsigned int thePoint; ///< Point number (on trajectory)
unsigned int theDWMethod; ///< Down-weighting method (0: None, 1: Tukey, 2: Huber, 3: Cauchy)
double theDownWeight; ///< Down-weighting factor (0-1)
double thePrediction; ///< Prediction from fit
// standard local parameters (curvature, offsets), fixed size
......
......@@ -45,7 +45,7 @@ namespace gbl {
GblData::GblData(unsigned int aLabel, dataBlockType aType, double aValue,
double aPrec, unsigned int aTraj, unsigned int aPoint) :
theLabel(aLabel), theRow(0), theType(aType), theValue(aValue), thePrecision(
aPrec), theTrajectory(aTraj), thePoint(aPoint), theDownWeight(
aPrec), theTrajectory(aTraj), thePoint(aPoint), theDWMethod(0), theDownWeight(
1.), thePrediction(0.), theNumLocal(0), moreParameters(), moreDerivatives() {
}
......@@ -92,6 +92,7 @@ void GblData::setPrediction(const VVector &aVector) {
*/
double GblData::setDownWeighting(unsigned int aMethod) {
theDWMethod = aMethod;
double aWeight = 1.;
double scaledResidual = fabs(theValue - thePrediction) * sqrt(thePrecision);
if (aMethod == 1) // Tukey
......@@ -117,11 +118,32 @@ double GblData::setDownWeighting(unsigned int aMethod) {
/// Calculate Chi2 contribution.
/**
* For down-weighting with M-estimators the corresponding objective function is used.
*
* \return (down-weighted) Chi2
*/
double GblData::getChi2() const {
double aDiff = theValue - thePrediction;
return aDiff * aDiff * thePrecision * theDownWeight;
double scaledResidual = fabs(theValue - thePrediction) * sqrt(thePrecision);
double chi2 = scaledResidual * scaledResidual;
if (theDWMethod == 1) // Tukey
{
if (scaledResidual < 4.6851) {
chi2 = (1.0
- pow(1.0 - 0.045558 * scaledResidual * scaledResidual, 3))
/ (3. * 0.045558);
} else {
chi2 = 1.0 / (3. * 0.045558);
}
} else if (theDWMethod == 2) //Huber
{
if (scaledResidual >= 1.345) {
chi2 = 1.345 * (2. * scaledResidual - 1.345);
}
} else if (theDWMethod == 3) //Cauchy
{
chi2 = log(1.0 + (scaledResidual * scaledResidual / 5.6877)) * 5.6877;
}
return chi2;
}
/// Print data block.
......
......@@ -315,6 +315,8 @@ class GblData(object):
self.__value = aValue
## precision (diagonal element of inverse covariance matrix); float
self.__precision = aPrec
## down weighting method; int
self.__dwMethod = 0
## down weighting factor (M-estimators); float
self.__downWeight = 1.
## prediction (for value from fit); float
......@@ -402,10 +404,11 @@ class GblData(object):
# @return weight (0..1); float
#
def setDownWeighting(self, aMethod):
self.__dwMethod = aMethod
scaledResidual = abs(self.__value - self.__prediction) * math.sqrt(self.__precision)
if (aMethod == 1): # Tukey
if (scaledResidual < 4.6851):
aWeight = (1.0 - 0.045558 * scaledResidual ** 2) ** 2
aWeight = (1.0 - (scaledResidual / 4.6851) ** 2) ** 2
else:
aWeight = 0.
elif (aMethod == 2): # Huber
......@@ -419,11 +422,24 @@ class GblData(object):
return aWeight
## Calculate Chi2 (contribution) from data.
#
#
# For down-weighting with M-estimators the corresponding objective function is used.
#
# @return Chi2; float
#
def getChi2(self):
Chi2 = (self.__value - self.__prediction) ** 2 * self.__precision * self.__downWeight
scaledResidual = abs(self.__value - self.__prediction) * math.sqrt(self.__precision)
Chi2 = scaledResidual ** 2
if (self.__dwMethod == 1): # Tukey
if (scaledResidual < 4.6851):
Chi2 = 4.6851 ** 2 / 3. * (1. - (1. - (scaledResidual / 4.6851) ** 2) ** 3)
else:
Chi2 = 4.6851 ** 2 / 3.
elif (self.__dwMethod == 2): # Huber
if (scaledResidual > 1.345):
Chi2 = 1.345 * (2.*scaledResidual - 1.345)
elif (self.__dwMethod == 3): # Cauchy
Chi2 = math.log(1. + (scaledResidual / 2.3849) ** 2) * 2.3849 ** 2
return Chi2
## Get Label.
......
......@@ -73,7 +73,7 @@ def example1():
measErr = np.array([ 0.001, 0.001]) # 10 mu
measPrec = 1.0 / measErr ** 2
# scattering error
scatErr = 0.001 # 1 mread
scatErr = 0.001 # 1 mrad
# RMS of CurviLinear track parameters (Q/P, slopes, offsets)
clErr = np.array([0.001, -0.1, 0.2, -0.15, 0.25])
# precision matrix for external seed (in local system)
......@@ -158,8 +158,9 @@ def example1():
# point with scatterer
jac = np.dot(proC2l, np.dot(jacPointToPoint, proL2c))
point = GblPoint(jac)
scatP = local.getScatPrecision(scatErr)
point.addScatterer([scat, scatP])
if scatErr > 0:
scatP = local.getScatPrecision(scatErr)
point.addScatterer([scat, scatP])
iLabel = traj.addPoint(point)
if iLabel == abs(seedLabel):
clSeed = np.linalg.inv(clCov)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment