fit: stable variance inversion and SVD dynamic range

lib/Statistics/XYSampleData.cpp

@@ -62,6 +62,11 @@ double SampleFitResult::getPValue(const Index s) const
     return Math::chi2PValue(getChi2(s), getNDof());
 }
 
+double SampleFitResult::getSvdRangeDb(void) const
+{
+    return svdRangeDb_;
+}
+
 double SampleFitResult::getCcdf(const Index s) const
 {
     return Math::chi2Ccdf(getChi2(s), getNDof());
@@ -107,9 +112,11 @@ void SampleFitResult::print(const bool printXsi, ostream &out) const
         getChi2(), static_cast<int>(getNDof()), getChi2PerDof(), getCcdf(), 
         getPValue());
     out << buf << endl;
+    sprintf(buf, "correlation dynamic range= %.1f dB", getSvdRangeDb());
+    out << buf << endl;
     for (Index p = 0; p < pMax; ++p)
     {
-        sprintf(buf, "%8s= % e +/- %e", parName_[p].c_str(),
+        sprintf(buf, "%12s= % e +/- %e", parName_[p].c_str(),
                 (*this)[central](p), err(p));
         out << buf << endl;
     }
@@ -249,6 +256,20 @@ const DMat & XYSampleData::getFitVarMatPInv(void)
     return data_.getFitVarMatPInv();
 }
 
+const DMat & XYSampleData::getFitCorrMat(void)
+{
+    computeVarMat();
+    
+    return data_.getFitCorrMat();
+}
+
+const DMat & XYSampleData::getFitCorrMatPInv(void)
+{
+    computeVarMat();
+    
+    return data_.getFitCorrMatPInv();
+}
+
 // set data to a particular sample /////////////////////////////////////////////
 void XYSampleData::setDataToSample(const Index s)
 {
@@ -319,9 +340,10 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
         }
     }
     minimizer.back()->setVerbosity(verbCopy);
-    result.nPar_    = sampleResult.getNPar();
+    result.nPar_       = sampleResult.getNPar();
-    result.nDof_    = sampleResult.nDof_;
+    result.nDof_       = sampleResult.nDof_;
-    result.parName_ = sampleResult.parName_;
+    result.parName_    = sampleResult.parName_;
+    result.svdRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
     
     return result;
 }

lib/Statistics/XYSampleData.hpp

@@ -49,6 +49,7 @@ public:
     double                       getNDof(void) const;
     Index                        getNPar(void) const;
     double                       getPValue(const Index s = central) const;
+    double                       getSvdRangeDb(void) const;
     double                       getCcdf(const Index s = central) const;
     const DoubleFunction &       getModel(const Index s = central,
                                           const Index j = 0) const;
@@ -60,6 +61,7 @@ public:
                std::ostream &out = std::cout) const;
 private:
     DSample                           chi2_;
+    double                            svdRangeDb_{0.};
     Index                             nDof_{0}, nPar_{0};
     std::vector<DoubleFunctionSample> model_;
     std::vector<std::string>          parName_;
@@ -91,9 +93,11 @@ public:
     const DMat &       getXYVar(const Index i, const Index j);
     DVec               getXError(const Index i);
     DVec               getYError(const Index j);
-    // get total fit variance matrix and its pseudo-inverse
+    // get total fit variance & correlation matrices and their pseudo-inverse
     const DMat & getFitVarMat(void);
     const DMat & getFitVarMatPInv(void);
+    const DMat & getFitCorrMat(void);
+    const DMat & getFitCorrMatPInv(void);
     // set data to a particular sample
     void setDataToSample(const Index s);
     // get internal XYStatData

lib/Statistics/XYStatData.cpp

@@ -60,6 +60,11 @@ double FitResult::getCcdf(void) const
     return Math::chi2Ccdf(getChi2(), getNDof());;
 }
 
+double FitResult::getSvdRangeDb(void) const
+{
+    return svdRangeDb_;
+}
+
 const DoubleFunction & FitResult::getModel(const Index j) const
 {
     return model_[j];
@@ -75,9 +80,11 @@ void FitResult::print(const bool printXsi, ostream &out) const
         getChi2(), static_cast<int>(getNDof()), getChi2PerDof(), getCcdf(), 
         getPValue());
     out << buf << endl;
+    sprintf(buf, "correlation dynamic range= %.1f dB", getSvdRangeDb());
+    out << buf << endl;
     for (Index p = 0; p < pMax; ++p)
     {
-        sprintf(buf, "%8s= %e", parName_[p].c_str(), (*this)(p));
+        sprintf(buf, "%12s= %e", parName_[p].c_str(), (*this)(p));
         out << buf << endl;
     }
 }
@@ -259,6 +266,20 @@ const DMat & XYStatData::getFitVarMatPInv(void)
     return fitVarInv_;
 }
 
+const DMat & XYStatData::getFitCorrMat(void)
+{
+    updateFitVarMat();
+    
+    return fitCorr_;
+}
+
+const DMat & XYStatData::getFitCorrMatPInv(void)
+{
+    updateFitVarMat();
+    
+    return fitCorrInv_;
+}
+
 // fit /////////////////////////////////////////////////////////////////////////
 FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
                           const vector<const DoubleModel *> &v)
@@ -337,9 +358,10 @@ FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
         result    = (*m)(chi2);
         totalInit = result;
     }
-    result.chi2_ = chi2(result);
+    result.svdRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
-    result.nPar_ = nPar;
+    result.chi2_       = chi2(result);
-    result.nDof_ = layout.totalYSize - nPar;
+    result.nPar_       = nPar;
+    result.nDof_       = layout.totalYSize - nPar;
     result.model_.resize(v.size());
     for (unsigned int j = 0; j < v.size(); ++j)
     {
@@ -551,8 +573,11 @@ void XYStatData::updateFitVarMat(void)
         chi2DataVec_.resize(layout.totalSize);
         chi2ModVec_.resize(layout.totalSize);
         chi2Vec_.resize(layout.totalSize);
-        fitVar_    = fitVar_.cwiseProduct(makeCorrFilter());
+        fitVar_     = fitVar_.cwiseProduct(makeCorrFilter());
-        fitVarInv_ = fitVar_.pInverse(getSvdTolerance());
+        fitCorr_    = Math::varToCorr(fitVar_);
+        fitCorrInv_ = fitCorr_.pInverse(getSvdTolerance());
+        fitVarInv_  = Math::corrToVar(fitCorrInv_, fitVar_.diagonal().cwiseInverse());
+
         scheduleFitVarMatInit(false);
     }
 }

lib/Statistics/XYStatData.hpp

@@ -48,12 +48,13 @@ public:
     Index                  getNPar(void) const;
     double                 getPValue(void) const;
     double                 getCcdf(void) const;
+    double                 getSvdRangeDb(void) const;
     const DoubleFunction & getModel(const Index j = 0) const;
     // IO
     void print(const bool printXsi = false,
                std::ostream &out = std::cout) const;
 private:
-    double                      chi2_{0.};
+    double                      chi2_{0.}, svdRangeDb_{0.};
     Index                       nDof_{0}, nPar_{0};
     std::vector<DoubleFunction> model_;
     std::vector<std::string>    parName_;
@@ -88,9 +89,11 @@ public:
     DVec           getXError(const Index i) const;
     DVec           getYError(const Index j) const;
     DMat           getTable(const Index i, const Index j) const;
-    // get total fit variance matrix and its pseudo-inverse
+    // get total fit variance & correlation matrices and their pseudo-inverse
     const DMat & getFitVarMat(void);
     const DMat & getFitVarMatPInv(void);
+    const DMat & getFitCorrMat(void);
+    const DMat & getFitCorrMatPInv(void);
     // fit
     FitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
                   const std::vector<const DoubleModel *> &v);
@@ -131,7 +134,7 @@ private:
     std::vector<DVec>                    xData_;
     std::vector<DVec>                    xMap_;
     Mat<DMat>                            xxVar_, yyVar_, xyVar_;
-    DMat                                 fitVar_, fitVarInv_;
+    DMat                                 fitVar_, fitVarInv_, fitCorr_, fitCorrInv_;
     DVec                                 chi2DataVec_, chi2ModVec_, chi2Vec_;
     DVec                                 xBuf_;
     bool                                 initXMap_{true};