Fitter::FitFCN and no FitResult::GetConfidenceIntervals

Hi,

I’ve defined my own chi2 function (myfcn.C) and model (mymodel.C) for some data and I would like to calculate confidence intervals like this (test.py):

[code]import ROOT
from array import array

ROOT.gROOT.ProcessLine(".L mymodel.C++g")
ROOT.gROOT.ProcessLine(".L myfcn.C++g")

x = array(‘d’,[1.])
p = array(‘d’,[1,2])

fitter = ROOT.Fit.Fitter()
fitdata = ROOT.Fit.BinData(1000,1,ROOT.Fit.BinData.kNoError)
datasize = 5
ci = array(‘d’,[0]*datasize)
for i in range(datasize):
fitdata.Add(i,i)

mymodel = ROOT.MyModel()
myfun = ROOT.MyFCN(fitdata,mymodel)
fitter.FitFCN(myfun,p,datasize)
fitter.Result().GetConfidenceIntervals(fitdata, ci,0.95)
print ci[/code]

output of this short script is (python test.py):

Info in <TUnixSystem::ACLiC>: creating shared library /home/jiri/workspace/physics/elscat/./mymodel_C.so Info in <TUnixSystem::ACLiC>: creating shared library /home/jiri/workspace/physics/elscat/./myfcn_C.so Error in <ROOT::Math::FitResult::GetConfidenceIntervals>: cannot compute Confidence Intervals without fitter function array('d', [0.0, 0.0, 0.0, 0.0, 0.0])

i.e. no confidence intervals as I would expect… :frowning:

The problem is that Fitter::FitFCN() function sets the “fitter function” to 0 before calling Fitter::DoMinimization<> and then FitResult gives error if one tries to call GetConfidenceIntervals() (see the code root.cern.ch/viewvc/trunk/math/m … iew=markup).

Maybe I missed something basic but the code confuse me quite a lot :confused: . Is there any possibility how to calculate the confidence intervals?

Thank you in advance,
Jiri
test.py (499 Bytes)
myfcn.C (313 Bytes)
mymodel.C (1.04 KB)

Hi,
the confidence intervals are computed on the model function used for fitting. The fit model function needs to be known in order to compute these intervals. The problem is, if you use FitFCN, that you are passing to the fitter class just the objective function (i.e. the least square or likelihood function ).
I would probably need to change the interface in order to be able to pass in addition to the objective function also the model function, so it can be then later. I could probably add optionally the possibility to pass a pointer to the model function in FitFCN.
For curiosity, why do you use FitFCN ? Do you need to use a customized least square or likelihood function ?

Best Regards

Lorenzo

Hi,

Ok, I looked to the code and I see the problem now. In my case MyFCN class is derived from ROOT::Fit::Chi2Function and contains both the objective and model function but if FitFCN is called than really just the objective function is passed to the fitter class.

Yes, I have some data, a model function and a customized least square function.The Fitter class looks like realy usefull for fitting because with FitFCN method I can minimize my customized least square function. Moreover, common interface for different minimizers, algorithms, possibility to set max function/iteration calls, etc. - interact more with the minimization process is exatly what I need. I’m just looking for a way how to calculate also the confidence intervals for a model function once all the free parameters were fitted to data using a customized least square function.

Cheers,
Jiri

Hi,

if I try to overload virtual double Chi2Function::DoEval (const double * x) function like (see updated myfcn.C):

#include "Fit/Chi2FCN.h"

class MyFCN: public ROOT::Fit::Chi2Function{
  public:
    MyFCN(const ROOT::Fit::BinData&  data, const ROOT::Fit::Chi2Function::IModelFunction&  func): ROOT::Fit::Chi2Function(data,func) {}
    ~MyFCN(){}

    private:
    virtual double  DoEval(const double* x) const{
        double ret = ROOT::Fit::Chi2Function::DoEval(x);
        std::cout << "MyFCN::DoEval" << ret << endl;
        assert(0);
        return ret;
    }
};

and try to perform the fit similarly as in the previous post (see updated test.py) then this overloaded MyFCN::DoEval function is not called at all and it seems that Chi2Function::DoEval is called istead.

Why is MyFCN::DoEval not called at all? Is it a proper way how to use the class Fitter and the method FitFCN to perform a fit with customized least square function?

Best Regards,
Jiri
test.py (541 Bytes)
mymodel.C (1.08 KB)
myfcn.C (455 Bytes)

Hi,

to re-implement the Chi2Function class you need to re-implement also the Clone() method, otherwise when the function is copied by the minimization program it will be sliced.
There is also another potential problem, you cannot call from the MyFCN class the method Chi2Function::DoEval since it is private. If you work using CINT (als also Python) it will work, but not in a C++ only compiled program. In that case if you need to use that method, the best is to have the Chi2Function class used as a member of your class. If you however want to customize the least square function, probably you don;t need to call Chi2Function::DoEval.

I attached below a corrected implementation of MyFCN (with the Clone method)

Best Regards

Lorenzo
myfcn.C (730 Bytes)

Hi,

thank you very much for correcting the MyFCN class. It works without problems after re-implementing also the Clone() method :smiley: .

Yes, I want to customize the least square function, so I don’t need to call Chi2Function::DoEval (I may call public Chi2Function::operator() which call the DoEval method instead…).

Best Regards,
Jiri

Hi,

I found this post since I was trying to do the same.

In may case the “combined” fit is done more or less exactly how is done in

root.cern.ch/root/html/tutorials … Fit.C.html

and I also obtain the error about the absence of the model function.

For me is not so easy to modify my code following the prescriptions said here since really I was starting differently (with the struct as in the example and not with a class, and by the way I don’t find the documentation for FitFCN(int npars, pointer to struct, …).

I tried also to pass to FitResult the model via ROOT::Fit::Chi2Function::ModelFunction() but I’m not able since any object I try to pass or is const while a non-const is needed or, If I try to build a non-const one via the copy constructor I cannot since everything I try is a pure abstract class…

Where I’m wrong? Can somebody help me?

Thanks,
Matteo

Hi,
The problem is that the FitResults when computing the confidence intervals needs to know the model function used and its parameters (i.e. their errors and their correlation matrix).
Now, in an example as in combinedFit.C you have actually two model functions, depending on the data set used.
A possible solution is to have an extended interface in FitResult::GetConfidenceInterval allowing the user to pass a model function and also a list of the parameter names (or indices) as they are defined in the FitResult class.

The method ROOT::Fit::Chi2Function::ModelFunction() is a method to return a pointer to the function and not to set a new object in the FitResult class.

If you need it, I could implement this extended interface in the next days

Best Regards

Lorenzo

Thanks a lot for the answer.

I have to admit that I didn’t understand completely the whole fitting machinery but I’m happy that at the end it revealed to be not my fault but it was really not possible!

I don’t want to ask you to implement this just for my stupid fits (is for an AMS paper but I’m even not sure if the collaboration would like to put this band at all or if the one I made in the while [just from error propagation and so a “variance” band and not a CL one…] cannot be enough).
Viceversa if you really plan to implement this since it could be useful in general it will be rally great!

Thanks a lot!
Matteo

Ciao,

reading again your answer I would like to clarify what I’m doing because I think there will be still another complication…
You mentioned about the possibility of a new interface allowing to pass a list of parameters (names or indices). I’m afraid this is not enough…

In your combined fit example you assume 6 parameters, shared between two distributions and one (the B) is using just two parameters, while the other (the S+B) is using five parameters, with one in common…

In my case instead I have 7 parameters and one component is using all of them while the other is using 5 parameters coming from the combination of the previous 7. And in facts my struct is this:

struct GlobalChi2 { 
  
  GlobalChi2(  ROOT::Math::IMultiGenFunction & f1,  
	       ROOT::Math::IMultiGenFunction & f2) : 
    fChi2_1(&f1), fChi2_2(&f2) {}
  
  double operator() (const double *par) const {    
    double p0[7];//flux
    p0[0] = par[0];//Cele
    p0[1] = par[1];//gele
    p0[2] = par[2];//Cs
    p0[3] = par[3];//gs
    p0[4] = par[4];//1.0/Es
    p0[5] = par[5];//Cpos
    p0[6] = par[6];//gpos

    double p1[5];//frac
    p1[0] = par[1]-par[6];//gele-gpos
    p1[1] = par[1]-par[3];//gele-gs
    p1[2] = par[5]/par[0];//Cpos/Cele
    p1[3] = par[2]/par[0];//Cs/Cele
    p1[4] = par[4];//1.0/Es

    return (*fChi2_1)(p0) + (*fChi2_2)(p1);
  } 
  
  const  ROOT::Math::IMultiGenFunction * fChi2_1;
  const  ROOT::Math::IMultiGenFunction * fChi2_2;
};

and while to “draw” the TF1 model correctly for the 7-parameters component was enough to do:

  int ppp[7] = {0, 1, 2, 3, 4, 5, 6};
  func->SetFitResult(result, ppp);

for the 5-parameters one I had to “replicate” by hand the SetFitResult method computing the errors by hand:

  if (result.IsEmpty()) { 
    Warning("SetFitResult","Empty Fit result - nathing is set in TF1");
    return NULL;      
  }
  if (func->GetNpar()!=5)) { 
    Error("SetFitResult","Invalid Fit result passed - number of parameter is %d , different than TF1::GetNpar() = %d", func->GetNpar(), result.NPar());
    return NULL;
  }
  if (result.Chi2() > 0) 
    func->SetChisquare(result.Chi2());
  else 
    func->SetChisquare(result.MinFcnValue());
  
  func->SetNDF(result.Ndf());
  //  printf("---------- %f\n", (double)(result.NFreeParameters()));
  func->SetNumberFitPoints(result.Ndf()+result.NFreeParameters());
  
  double* p1;
  double* errp1;
  int frac_params = 5;

  p1 = new double[frac_params];
  p1[0] = result.Parameter(1)-result.Parameter(6);//gele-gpos
  p1[1] = result.Parameter(1)-result.Parameter(3);//gele-gs
  p1[2] = result.Parameter(5)/result.Parameter(0);//Cpos/Cele
  p1[3] = result.Parameter(2)/result.Parameter(0);//Cs/Cele
  p1[4] = result.Parameter(4);//1.0/Es

  errp1 = new double[frac_params];
  errp1[0] = sqrt(pow(result.ParError(1), 2) + pow(result.ParError(6), 2));//gele-gpos
  errp1[1] = sqrt(pow(result.ParError(1), 2) + pow(result.ParError(3), 2));//gele-gs
  errp1[2] = p1[2]*sqrt(pow(result.ParError(5)/result.Parameter(5), 2) + pow(result.ParError(0)/result.Parameter(0), 2));//Cpos/Cele
  errp1[3] = p1[3]*sqrt(pow(result.ParError(2)/result.Parameter(2), 2) + pow(result.ParError(0)/result.Parameter(0), 2));//Cs/Cele
  errp1[4] = result.ParError(4);//1.0/Es
 
  for (int ii=0; ii<frac_params; ii++) {  
    func->SetParameter(ii, p1[ii]);
    func->SetParError(ii, errp1[ii]);
  }

since I didn’t find the way to pass to SetFitResults the ‘combinations’ of parameters I needed…

So in somehow I would have need the possibility to pass strings of parameters that should be interpreted not only as the names of the parameters but should be parsed as formulas… I imagine this is not so straightforward…

Thanks again,
Matteo

Matteo

Hi,

Thanks for your comment. Yes, in general you can have as model function what-ever subset and combinations of the parameters. Probably the best we could do is to generalise the code in FitResult::GetConfidenceLevel to work on any given function, parameter values and correlation matrix. Then you can compute the correlation matrix for your function parameters from the one obtained from FitResult using the given parameter relations and error propagation.

Lorenzo

Ciao,

mmm I didn’t understand completely…

I tried to think how to do the complete propagation error in my case and for the errors of the 5 combination parameters, from the all 7, is quite straightforward (in the code I posted before I did but I forgot, for laziness, the second order terms with the correlations…).

Instead for the propagation to the model I would need the derivatives of the model wrt to the parameters… Is this second part that you would like to implement inside GetConfidence Interval? I understood as the user should compute the correlation matrix of the particular combination under study (in my case the one with 5 parameters) from the complete one, by itself and then could pass it, with the parameters values and the model to a generic GetConfidenceLevel method (even a static one?) computing the confidence level. But if I see this possible for the variance (simple propagation error) but not for the confidence level, for which is needed to have the full chi-sq (or likelihood) profile. Isn’t it?

Matteo

Hi,

The derivatives of the model wrt to the parameters are currently computed numerically in FitResult::GetConfidenceIntervals.

In your case, I think, if your function depends on the parameter q, which are a function of the fitted parameter p, the ones stored in FitResult, you should pass to the CondidenceLevel routine not the function f(q) , but the function g§ = f ( q § ). In this way the code will compute automatically the derivatives for you.

One more thing, the GetConfidenceLevel routine does simple error propagation without using the full likelihood (chi square profile). The profile likelihood is used only when computing the confidence interval in the single parameter (using Minos) and it is returned in FitResult::LowerError(i) and FitResult::UpeerError(i).

So if you want a 95% band on the function values, it is computed rescaling the error assuming a parabolic chi square. In case your chi2 comes very different than one, the errors are re-normalised using the observed chi2 value.

Best Regards

Lorenzo

Hi,

thanks a lot for the precious explanations.
Now I understand better the various pieces and I know where to look to create my own version of the GetConfidenceLevel.

I didn’t know how the confidence level could have been computed and now it’s much more clear.

Thanks!
Matteo