Parent Directory
| 
Revision Log
import changes from math development branches for subdirectory math. List of changes in detail:
mathcore:
---------
MinimizerOptions:
new class for storing Minimizer option, with static default values that can be
changed by the user
FitConfig:
- use default values from MinimizerOption class
- rename method to create parameter settings from a function
FitUtil.cxx:
improve the derivative calculations used in the effective chi2 and in Fumili and
fix a bug for evaluation of likelihood or chi2 terms.
In EvaluatePdf() work and return the log of the pdf.
FitResult:
- improve the class by adding extra information like, num. of free parameters,
minimizer status, global correlation coefficients, information about fixed
and bound parameters.
- add method for getting fit confidence intervals
- improve print method
DataRange:
add method SetRange to distinguish from AddRange. SetRange deletes the existing
ranges.
ParamsSettings: make few methods const
FCN functions (Chi2FCN, LogLikelihoodFCN, etc..)
move some common methods and data members in base class (FitMethodFunction)
RootFinder: add template Solve() for any callable function.
mathmore:
--------
minimizer classes: fill status information
GSLNLSMinimizer: return error and covariance matrix
minuit2:
-------
Minuit2Minimizer: fill status information
DavidonErrorUpdator: check that delgam or gvg are not zero ( can happen when dg = 0)
FumiliFCNAdapter: work on the log of pdf
minuit:
-------
TLinearMinimizer: add support for robust fitting
TMinuitMinimizer: fill status information and fix a bug in filling the correlation matrix.
fumili:
------
add TFumiliMinimizer:
wrapper class for TFumili using Minimizer interface
// @(#)root/mathcore:$Id$
// Author: L. Moneta Wed Aug 30 11:05:34 2006
/**********************************************************************
* *
* Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
* *
* *
**********************************************************************/
// Implementation file for class FitResult
#include "Fit/FitResult.h"
#include "Fit/FitConfig.h"
#include "Fit/BinData.h"
#include "Math/Minimizer.h"
#include "Math/IParamFunction.h"
#include "Math/OneDimFunctionAdapter.h"
#include "Math/DistFunc.h"
#include "TMath.h"
#include "Math/RichardsonDerivator.h"
#include "Math/Error.h"
#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
namespace ROOT {
namespace Fit {
FitResult::FitResult() :
fValid(false), fNormalized(false), fNFree(0), fNdf(0), fNCalls(0),
fStatus(-1), fVal(0), fEdm(0), fChi2(0), fFitFunc(0)
{
// Default constructor implementation.
}
FitResult::FitResult(ROOT::Math::Minimizer & min, const FitConfig & fconfig, IModelFunction & func, bool isValid, unsigned int sizeOfData, const ROOT::Math::IMultiGenFunction * chi2func, bool minosErr, unsigned int ncalls ) :
fValid(isValid),
fNormalized(false),
fNFree(min.NFree() ),
fNCalls(min.NCalls()),
fStatus(min.Status() ),
fVal (min.MinValue()),
fEdm (min.Edm()),
fFitFunc(&func),
fParams(std::vector<double>(min.X(), min.X() + min.NDim() ) )
{
// Constructor from a minimizer, fill the data. ModelFunction is passed as non const
// since it will be managed by the FitResult
if (sizeOfData > 0) fNdf = sizeOfData - min.NFree();
// set right parameters in function (in case minimizer did not do before)
// do also when fit is not valid
fFitFunc->SetParameters(&fParams.front());
if (min.Errors() != 0)
fErrors = std::vector<double>(min.Errors(), min.Errors() + min.NDim() ) ;
// check for fixed or limited parameters
for (unsigned int ipar = 0; ipar < fParams.size(); ++ipar) {
const ParameterSettings & par = fconfig.ParSettings(ipar);
if (par.IsFixed() ) fFixedParams.push_back(ipar);
if (par.IsBound() ) fBoundParams.push_back(ipar);
}
if (chi2func == 0)
fChi2 = fVal;
else {
// compute chi2 equivalent
fChi2 = (*chi2func)(&fParams[0]);
}
// replace ncalls if given (they are taken from the FitMethodFunction)
if (ncalls !=0) fNCalls = ncalls;
unsigned int n = min.NDim();
// // fill error matrix
// cov matrix rank
if (fValid) {
unsigned int r = n * ( n + 1 )/2;
fCovMatrix.reserve(r);
for (unsigned int i = 0; i < n; ++i)
for (unsigned int j = 0; j <= i; ++j)
fCovMatrix.push_back(min.CovMatrix(i,j) );
assert (fCovMatrix.size() == r );
// normalize errors if requested in configuration
if (fconfig.NormalizeErrors() ) NormalizeErrors();
// minos errors
if (minosErr) {
fMinosErrors.reserve(n);
for (unsigned int i = 0; i < n; ++i) {
double elow, eup;
bool ret = min.GetMinosError(0, elow, eup);
if (ret) fMinosErrors.push_back(std::make_pair(elow,eup) );
else fMinosErrors.push_back(std::make_pair(0.,0.) );
}
}
// globalCC
fGlobalCC.reserve(n);
for (unsigned int i = 0; i < n; ++i) {
double globcc = min.GlobalCC(i);
if (globcc < 0) break; // it is not supported by that minimizer
fGlobalCC.push_back(globcc);
}
}
fMinimType = fconfig.MinimizerType();
// append algorithm name for minimizer that support it
if ( (fMinimType.find("Fumili") == std::string::npos) &&
(fMinimType.find("GSLMultiFit") == std::string::npos)
) {
if (fconfig.MinimizerAlgoType() != "") fMinimType += " / " + fconfig.MinimizerAlgoType();
}
}
FitResult::FitResult(const FitResult &rhs) {
// Implementation of copy constructor
(*this) = rhs;
}
FitResult & FitResult::operator = (const FitResult &rhs) {
// Implementation of assignment operator.
if (this == &rhs) return *this; // time saving self-test
// Manages the fitted function
if (fFitFunc) delete fFitFunc;
fFitFunc = 0;
if (rhs.fFitFunc != 0 ) {
fFitFunc = dynamic_cast<IModelFunction *>( (rhs.fFitFunc)->Clone() );
assert(fFitFunc != 0);
}
// copy all other data members
fValid = rhs.fValid;
fNormalized = rhs.fNormalized;
fNFree = rhs.fNFree;
fNdf = rhs.fNdf;
fNCalls = rhs.fNCalls;
fStatus = rhs.fStatus;
fVal = rhs.fVal;
fEdm = rhs.fEdm;
fChi2 = rhs.fChi2;
fFixedParams = rhs.fFixedParams;
fBoundParams = rhs.fBoundParams;
fParams = rhs.fParams;
fErrors = rhs.fErrors;
fCovMatrix = rhs.fCovMatrix;
fGlobalCC = rhs.fGlobalCC;
fMinosErrors = rhs.fMinosErrors;
fMinimType = rhs.fMinimType;
return *this;
}
void FitResult::NormalizeErrors() {
// normalize errors and covariance matrix according to chi2 value
if (fNdf == 0 || fChi2 <= 0) return;
double s2 = fChi2/fNdf;
double s = std::sqrt(fChi2/fNdf);
for (unsigned int i = 0; i < fErrors.size() ; ++i)
fErrors[i] *= s;
for (unsigned int i = 0; i < fCovMatrix.size() ; ++i)
fCovMatrix[i] *= s2;
fNormalized = true;
}
double FitResult::Prob() const {
// fit probability
return ROOT::Math::chisquared_cdf_c(fChi2, static_cast<double>(fNdf) );
}
int FitResult::Index(const std::string & name) const {
// find index for given parameter name
unsigned int npar = fParams.size();
for (unsigned int i = 0; i < npar; ++i)
if ( fFitFunc->ParameterName(i) == name) return i;
return -1; // case name is not found
}
bool FitResult::IsParameterBound(unsigned int ipar) const {
for (unsigned int i = 0; i < fBoundParams.size() ; ++i)
if ( fBoundParams[i] == ipar) return true;
return false;
}
bool FitResult::IsParameterFixed(unsigned int ipar) const {
for (unsigned int i = 0; i < fFixedParams.size() ; ++i)
if ( fFixedParams[i] == ipar) return true;
return false;
}
void FitResult::Print(std::ostream & os, bool doCovMatrix) const {
// print the result in the given stream
// need to add minos errors , globalCC, etc..
if (!fValid) {
os << "\n****************************************\n";
os << " Invalid FitResult ";
os << "\n****************************************\n";
return;
}
os << "\n****************************************\n";
//os << " FitResult \n\n";
os << "Minimizer is " << fMinimType << std::endl;
unsigned int npar = fParams.size();
const unsigned int nw = 25;
if (fVal != fChi2)
os << std::setw(nw) << std::left << "Likelihood" << " =\t" << fVal << std::endl;
os << std::setw(nw) << std::left << "Chi2" << " =\t" << fChi2<< std::endl;
os << std::setw(nw) << std::left << "NDf" << " =\t" << fNdf << std::endl;
if (fMinimType.find("Linear") == std::string::npos) { // no need to print this for linear fits
os << std::setw(nw) << std::left << "Edm" << " =\t" << fEdm << std::endl;
os << std::setw(nw) << std::left << "NCalls" << " =\t" << fNCalls << std::endl;
}
assert(fFitFunc != 0);
for (unsigned int i = 0; i < npar; ++i) {
os << std::setw(nw) << std::left << fFitFunc->ParameterName(i) << " =\t" << fParams[i];
if (IsParameterFixed(i) )
os << " \t(fixed)";
else {
if (fErrors.size() != 0)
os << " \t+/-\t" << fErrors[i];
if (IsParameterBound(i) )
os << " \t (limited)";
}
os << std::endl;
}
if (doCovMatrix) PrintCovMatrix(os);
}
void FitResult::PrintCovMatrix(std::ostream &os) const {
// print the covariance and correlation matrix
if (!fValid) return;
if (fCovMatrix.size() == 0) return;
// os << "****************************************\n";
os << "\n Covariance Matrix \n\n";
unsigned int npar = fParams.size();
const int kPrec = 5;
const int kWidth = 10;
const int parw = 12;
const int matw = kWidth+4;
os << std::setw(parw) << " " << "\t";
for (unsigned int i = 0; i < npar; ++i)
if (!IsParameterFixed(i) ) {
os << std::setw(matw) << fFitFunc->ParameterName(i) ;
}
os << std::endl;
for (unsigned int i = 0; i < npar; ++i) {
if (!IsParameterFixed(i) ) {
os << std::setw(parw) << std::left << fFitFunc->ParameterName(i) << "\t";
for (unsigned int j = 0; j < npar; ++j) {
if (!IsParameterFixed(j) ) {
os.precision(kPrec); os.width(kWidth); os << std::setw(matw) << CovMatrix(i,j);
}
}
os << std::endl;
}
}
// os << "****************************************\n";
os << "\n Correlation Matrix \n\n";
os << std::setw(parw) << " " << "\t";
for (unsigned int i = 0; i < npar; ++i)
if (!IsParameterFixed(i) ) {
os << std::setw(matw) << fFitFunc->ParameterName(i) ;
}
os << std::endl;
for (unsigned int i = 0; i < npar; ++i) {
if (!IsParameterFixed(i) ) {
os << std::setw(parw) << std::left << fFitFunc->ParameterName(i) << "\t";
for (unsigned int j = 0; j < npar; ++j) {
if (!IsParameterFixed(j) ) {
os.precision(kPrec); os.width(kWidth); os << std::setw(matw) << Correlation(i,j);
}
}
os << std::endl;
}
}
}
void FitResult::GetConfidenceIntervals(unsigned int n, unsigned int stride1, unsigned int stride2, const double * x, double * ci, double cl ) const {
// stride1 stride in coordinate stride2 stride in dimension space
// i.e. i-th point in k-dimension is x[ stride1 * i + stride2 * k]
// compute the confidence interval of the fit on the given data points
// the dimension of the data points must match the dimension of the fit function
// confidence intervals are returned in array ci
// use student quantile
//double t = - TMath::StudentQuantile((1.-cl)/2, f->GetNDF());
double t = TMath::StudentQuantile(0.5 + cl/2, fNdf);
double chidf = TMath::Sqrt(fChi2/fNdf);
if (!fFitFunc) {
MATH_ERROR_MSG("FitResult::GetConfidenceIntervals","cannot compute Confidence Intervals without fitter function");
return;
}
unsigned int ndim = fFitFunc->NDim();
unsigned int npar = fFitFunc->NPar();
std::vector<double> xpoint(ndim);
std::vector<double> grad(npar);
std::vector<double> vsum(npar);
// loop on the points
for (unsigned int ipoint = 0; ipoint < n; ++ipoint) {
for (unsigned int kdim = 0; kdim < ndim; ++kdim) {
unsigned int i = ipoint * stride1 + kdim * stride2;
assert(i < ndim*n);
xpoint[kdim] = x[ipoint * stride1 + kdim * stride2];
}
// calculate gradient of fitted function w.r.t the parameters
// check first if fFitFunction provides parameter gradient or not
// does not provide gradient
// t.b.d : skip calculation for fixed parameters
ROOT::Math::RichardsonDerivator d;
for (unsigned int ipar = 0; ipar < npar; ++ipar) {
ROOT::Math::OneDimParamFunctionAdapter<const ROOT::Math::IParamMultiFunction &> fadapter(*fFitFunc,&xpoint.front(),&fParams.front(),ipar);
d.SetFunction(fadapter);
grad[ipar] = d(fParams[ipar] ); // evaluate df/dp
}
// multiply covariance matrix with gradient
vsum.assign(npar,0.0);
for (unsigned int ipar = 0; ipar < npar; ++ipar) {
for (unsigned int jpar = 0; jpar < npar; ++jpar) {
vsum[ipar] += CovMatrix(ipar,jpar) * grad[jpar];
}
}
// multiply gradient by vsum
double r2 = 0;
for (unsigned int ipar = 0; ipar < npar; ++ipar) {
r2 += grad[ipar] * vsum[ipar];
}
double r = std::sqrt(r2);
ci[ipoint] = r * t * chidf;
}
}
void FitResult::GetConfidenceIntervals(const BinData & data, double * ci, double cl ) const {
// implement confidence intervals from a given bin data sets
// currently copy the data from Bindata.
// could implement otherwise directly
unsigned int ndim = data.NDim();
unsigned int np = data.NPoints();
std::vector<double> xdata( ndim * np );
for (unsigned int i = 0; i < np ; ++i) {
const double * x = data.Coords(i);
std::vector<double>::iterator itr = xdata.begin()+ ndim * i;
std::copy(x,x+ndim,itr);
}
// points are arraned as x0,y0,z0, ....xN,yN,zN (stride1=ndim, stride2=1)
GetConfidenceIntervals(np,ndim,1,&xdata.front(),ci,cl);
}
} // end namespace Fit
} // end namespace ROOT
| Subversion Admin | ViewVC Help |
| Powered by ViewVC 1.0.9 |