portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2025-06-18 23:37:05 +01:00
Code Releases Activity

Merge from AP's develop.

Browse Source
This commit is contained in:
AndrewYongZhenNing
2022-02-23 13:52:27 +00:00
parent db99951dd4 47d0b3f040
commit cd1aeac669
30 changed files with 689 additions and 192 deletions
Download Patch File Download Diff File Expand all files Collapse all files

76
lib/Statistics/MatSample.hpp
View File

@ -103,6 +103,10 @@ public:
const Index nCol);
// resize all matrices
void resizeMat(const Index nRow, const Index nCol);
// covariance matrix
Mat<T> covarianceMatrix(const MatSample<T> &sample) const;
Mat<T> varianceMatrix(void) const;
Mat<T> correlationMatrix(void) const;
};
// non-member operators
@ -379,6 +383,78 @@ void MatSample<T>::resizeMat(const Index nRow, const Index nCol)
}
}
// covariance matrix ///////////////////////////////////////////////////////////
template <typename T>
Mat<T> MatSample<T>::covarianceMatrix(const MatSample<T> &sample) const
{
if (((*this)[central].cols() != 1) or (sample[central].cols() != 1))
{
LATAN_ERROR(Size, "samples have more than one column");
}
Index n1 = (*this)[central].rows(), n2 = sample[central].rows();
Index nSample = this->size();
Mat<T> tmp1(n1, nSample), tmp2(n2, nSample), res(n1, n2);
Mat<T> s1(n1, 1), s2(n2, 1), one(nSample, 1);
one.fill(1.);
s1.fill(0.);
s2.fill(0.);
for (unsigned int s = 0; s < nSample; ++s)
{
s1 += (*this)[s];
tmp1.col(s) = (*this)[s];
}
tmp1 -= s1*one.transpose()/static_cast<double>(nSample);
for (unsigned int s = 0; s < nSample; ++s)
{
s2 += sample[s];
tmp2.col(s) = sample[s];
}
tmp2 -= s2*one.transpose()/static_cast<double>(nSample);
res = tmp1*tmp2.transpose()/static_cast<double>(nSample - 1);
return res;
}
template <typename T>
Mat<T> MatSample<T>::varianceMatrix(void) const
{
if ((*this)[central].cols() != 1)
{
LATAN_ERROR(Size, "samples have more than one column");
}
Index n1 = (*this)[central].rows();
Index nSample = this->size();
Mat<T> tmp1(n1, nSample), res(n1, n1);
Mat<T> s1(n1, 1), one(nSample, 1);
one.fill(1.);
s1.fill(0.);
for (unsigned int s = 0; s < nSample; ++s)
{
s1 += (*this)[s];
tmp1.col(s) = (*this)[s];
}
tmp1 -= s1*one.transpose()/static_cast<double>(nSample);
res = tmp1*tmp1.transpose()/static_cast<double>(nSample - 1);
return res;
}
template <typename T>
Mat<T> MatSample<T>::correlationMatrix(void) const
{
Mat<T> res = varianceMatrix();
Mat<T> invDiag(res.rows(), 1);
invDiag = res.diagonal();
invDiag = invDiag.cwiseInverse().cwiseSqrt();
res = (invDiag*invDiag.transpose()).cwiseProduct(res);
return res;
}
END_LATAN_NAMESPACE

148
lib/Statistics/StatArray.hpp
View File

@ -51,14 +51,11 @@ public:
const T & operator[](const Index s) const;
// statistics
void bin(Index binSize);
T sum(const Index pos = 0, const Index n = -1) const;
T meanOld(const Index pos = 0, const Index n = -1) const;
T mean(const Index pos = 0, const Index n = -1) const;
T covariance(const StatArray<T, os> &array, const Index pos = 0,
const Index n = -1) const;
T covarianceMatrix(const StatArray<T, os> &array, const Index pos = 0,
const Index n = -1) const;
T variance(const Index pos = 0, const Index n = -1) const;
T varianceMatrix(const Index pos = 0, const Index n = -1) const;
T correlationMatrix(const Index pos = 0, const Index n = -1) const;
T covariance(const StatArray<T, os> &array) const;
T variance(void) const;
// IO type
virtual IoType getType(void) const;
public:
@ -66,7 +63,7 @@ public:
};
// reduction operations
namespace ReducOp
namespace StatOp
{
// general templates
template <typename T>
@ -148,128 +145,67 @@ void StatArray<T, os>::bin(Index binSize)
}
}
template <typename T, Index os>
T StatArray<T, os>::sum(const Index pos, const Index n) const
{
T result;
const Index m = (n >= 0) ? n : size();
result = (*this)[pos];
for (Index i = pos + 1; i < pos + m; ++i)
{
result += (*this)[i];
}
return result;
}
template <typename T, Index os>
T StatArray<T, os>::mean(const Index pos, const Index n) const
{
T result = T();
const Index m = (n >= 0) ? n : size();
if (m)
return sum(pos, n)/static_cast<double>(m);
}
template <typename T, Index os>
T StatArray<T, os>::covariance(const StatArray<T, os> &array) const
{
T s1, s2, res;
s1 = array.sum();
s2 = this->sum();
res = StatOp::prod<T>(array[0], (*this)[0]);
for (Index i = 1; i < size(); ++i)
{
result = this->segment(pos+os, m).redux(&ReducOp::sum<T>);
res += StatOp::prod<T>(array[i], (*this)[i]);
}
return result/static_cast<double>(m);
}
template <typename T, Index os>
T StatArray<T, os>::covariance(const StatArray<T, os> &array, const Index pos,
const Index n) const
{
T s1, s2, prs, res = T();
const Index m = (n >= 0) ? n : size();
if (m)
{
auto arraySeg = array.segment(pos+os, m);
auto thisSeg = this->segment(pos+os, m);
s1 = thisSeg.redux(&ReducOp::sum<T>);
s2 = arraySeg.redux(&ReducOp::sum<T>);
prs = thisSeg.binaryExpr(arraySeg, &ReducOp::prod<T>)
.redux(&ReducOp::sum<T>);
res = prs - ReducOp::prod(s1, s2)/static_cast<double>(m);
}
return res/static_cast<double>(m - 1);
}
template <typename T, Index os>
T StatArray<T, os>::covarianceMatrix(const StatArray<T, os> &array,
const Index pos, const Index n) const
{
T s1, s2, prs, res = T();
const Index m = (n >= 0) ? n : size();
if (m)
{
auto arraySeg = array.segment(pos+os, m);
auto thisSeg = this->segment(pos+os, m);
s1 = thisSeg.redux(&ReducOp::sum<T>);
s2 = arraySeg.redux(&ReducOp::sum<T>);
prs = thisSeg.binaryExpr(arraySeg, &ReducOp::tensProd<T>)
.redux(&ReducOp::sum<T>);
res = prs - ReducOp::tensProd(s1, s2)/static_cast<double>(m);
}
return res/static_cast<double>(m - 1);
}
template <typename T, Index os>
T StatArray<T, os>::variance(const Index pos, const Index n) const
{
return covariance(*this, pos, n);
}
template <typename T, Index os>
T StatArray<T, os>::varianceMatrix(const Index pos, const Index n) const
{
return covarianceMatrix(*this, pos, n);
}
template <typename T, Index os>
T StatArray<T, os>::correlationMatrix(const Index pos, const Index n) const
{
T res = varianceMatrix(pos, n);
T invDiag(res.rows(), 1);
invDiag = res.diagonal();
invDiag = invDiag.cwiseInverse().cwiseSqrt();
res = (invDiag*invDiag.transpose()).cwiseProduct(res);
res -= StatOp::prod<T>(s1, s2)/static_cast<double>(size());
res /= static_cast<double>(size() - 1);
return res;
}
// reduction operations ////////////////////////////////////////////////////////
namespace ReducOp
template <typename T, Index os>
T StatArray<T, os>::variance(void) const
{
template <typename T>
inline T sum(const T &a, const T &b)
{
return a + b;
}
return covariance(*this);
}
// reduction operations ////////////////////////////////////////////////////////
namespace StatOp
{
template <typename T>
inline T prod(const T &a, const T &b)
{
return a*b;
}
template <typename T>
inline T tensProd(const T &v1 __dumb, const T &v2 __dumb)
{
LATAN_ERROR(Implementation,
"tensorial product not implemented for this type");
}
template <>
inline Mat<double> prod(const Mat<double> &a, const Mat<double> &b)
{
return a.cwiseProduct(b);
}
template <>
inline Mat<double> tensProd(const Mat<double> &v1,
const Mat<double> &v2)
{
if ((v1.cols() != 1) or (v2.cols() != 1))
{
LATAN_ERROR(Size,
"tensorial product is only valid with column vectors");
}
return v1*v2.transpose();
}
}
// IO type /////////////////////////////////////////////////////////////////////

61
lib/Statistics/XYSampleData.cpp
View File

@ -62,6 +62,11 @@ double SampleFitResult::getPValue(const Index s) const
return Math::chi2PValue(getChi2(s), getNDof());
}
double SampleFitResult::getCorrRangeDb(void) const
{
return corrRangeDb_;
}
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", getCorrRangeDb());
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)
{
@ -330,9 +351,10 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
{
computeVarMat();
SampleFitResult result;
FitResult sampleResult;
DVec initCopy = init;
SampleFitResult result;
FitResult sampleResult;
DVec initCopy = init;
Minimizer::Verbosity verbCopy = minimizer.back()->getVerbosity();
result.resize(nSample_);
result.chi2_.resize(nSample_);
@ -348,9 +370,11 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
result.model_[j][s] = sampleResult.getModel(j);
}
}
result.nPar_ = sampleResult.getNPar();
result.nDof_ = sampleResult.nDof_;
result.parName_ = sampleResult.parName_;
minimizer.back()->setVerbosity(verbCopy);
result.nPar_ = sampleResult.getNPar();
result.nDof_ = sampleResult.nDof_;
result.parName_ = sampleResult.parName_;
result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
return result;
}
@ -382,6 +406,29 @@ XYSampleData XYSampleData::getResiduals(const SampleFitResult &fit)
return res;
}
XYSampleData XYSampleData::getNormalisedResiduals(const SampleFitResult &fit)
{
XYSampleData res(*this);
for (Index j = 0; j < getNYDim(); ++j)
{
const DoubleFunctionSample &f = fit.getModel(_, j);
for (auto &p: yData_[j])
{
res.y(p.first, j) -= f(x(p.first));
}
const DMat &var = res.getYYVar(j, j);
for (auto &p: yData_[j])
{
res.y(p.first, j) /= sqrt(var(p.first, p.first));
}
}
return res;
}
XYSampleData XYSampleData::getPartialResiduals(const SampleFitResult &fit,
const DVec &ref, const Index i)
{

7
lib/Statistics/XYSampleData.hpp
View File

@ -49,6 +49,7 @@ public:
double getNDof(void) const;
Index getNPar(void) const;
double getPValue(const Index s = central) const;
double getCorrRangeDb(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 corrRangeDb_{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
@ -118,6 +122,7 @@ public:
const DoubleModel &model, const Ts... models);
// residuals
XYSampleData getResiduals(const SampleFitResult &fit);
XYSampleData getNormalisedResiduals(const SampleFitResult &fit);
XYSampleData getPartialResiduals(const SampleFitResult &fit, const DVec &x,
const Index i);
private:

60
lib/Statistics/XYStatData.cpp
View File

@ -60,6 +60,11 @@ double FitResult::getCcdf(void) const
return Math::chi2Ccdf(getChi2(), getNDof());;
}
double FitResult::getCorrRangeDb(void) const
{
return corrRangeDb_;
}
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", getCorrRangeDb());
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;
}
}
@ -216,7 +223,7 @@ DVec XYStatData::getXError(const Index i) const
DVec XYStatData::getYError(const Index j) const
{
checkXDim(j);
checkYDim(j);
return yyVar_(j, j).diagonal().cwiseSqrt();
}
@ -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.nPar_ = nPar;
result.nDof_ = layout.totalYSize - nPar;
result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
result.chi2_ = chi2(result);
result.nPar_ = nPar;
result.nDof_ = layout.totalYSize - nPar;
result.model_.resize(v.size());
for (unsigned int j = 0; j < v.size(); ++j)
{
@ -379,6 +401,27 @@ XYStatData XYStatData::getResiduals(const FitResult &fit)
return res;
}
XYStatData XYStatData::getNormalisedResiduals(const FitResult &fit)
{
XYStatData res(*this);
for (Index j = 0; j < getNYDim(); ++j)
{
const DoubleFunction &f = fit.getModel(j);
const DVec err = getYError(j);
Index row = 0;
for (auto &p: yData_[j])
{
res.y(p.first, j) -= f(x(p.first));
res.y(p.first, j) /= err(row);
row++;
}
}
return res;
}
XYStatData XYStatData::getPartialResiduals(const FitResult &fit,
const DVec &ref, const Index i)
{
@ -530,8 +573,11 @@ void XYStatData::updateFitVarMat(void)
chi2DataVec_.resize(layout.totalSize);
chi2ModVec_.resize(layout.totalSize);
chi2Vec_.resize(layout.totalSize);
fitVar_ = fitVar_.cwiseProduct(makeCorrFilter());
fitVarInv_ = fitVar_.pInverse(getSvdTolerance());
fitVar_ = fitVar_.cwiseProduct(makeCorrFilter());
fitCorr_ = Math::varToCorr(fitVar_);
fitCorrInv_ = fitCorr_.pInverse(getSvdTolerance());
fitVarInv_ = Math::corrToVar(fitCorrInv_, fitVar_.diagonal().cwiseInverse());
scheduleFitVarMatInit(false);
}
}

10
lib/Statistics/XYStatData.hpp
View File

@ -48,12 +48,13 @@ public:
Index getNPar(void) const;
double getPValue(void) const;
double getCcdf(void) const;
double getCorrRangeDb(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.}, corrRangeDb_{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);
@ -104,6 +107,7 @@ public:
const DoubleModel &model, const Ts... models);
// residuals
XYStatData getResiduals(const FitResult &fit);
XYStatData getNormalisedResiduals(const FitResult &fit);
XYStatData getPartialResiduals(const FitResult &fit, const DVec &ref,
const Index i);
protected:
@ -130,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};