Non Linear Least Square Fit with GSL MultiFit

Hello,

i am trying to build a non linear least square fit by using the GSL MultiFit Algorithm but i do not get along with it. Could the fitting experts of you please help me out?

I tried to make use of the “Plug-In Manager” for which an example can be found here:
root.cern.ch/drupal/content/nume … nimization

However, i write a little chisquare function similar as it suggest here for general fitting of datapoints:
root.cern.ch/drupal/content/how- … ata-points

When i try to run the program i get the following error message:

GSLNLSMinimizer: Invalid function set - only Chi2Func supported
Error in ROOT::Math::GSLNLSMinimizer::Minimize: Function has not been set

The function i have written is:

void Chi2Func(const double *xx){

    ostringstream FitData;
    FitData << "fit_mass_table_beta_1_qmass_0.01_tc.dat";

    const Int_t nbins = det_entries(FitData);
    Int_t i;

    //calculate chisquare
    Double_t chisq = 0;
    Double_t delta;
    for (i=0; i<nbins; i++) {
    delta  = (dMassFit.at(i)-(xx[0] + xx[1]/(dVolumeFit.at(i)*xx[2])+xx[1]*dChargeFit.at(i)/(dVolumeFit.at(i)*dVolumeFit.at(i)*xx[2]*xx[2])))/dErrorFit.at(i);
    chisq += delta*delta;
    }   
    f = chisq;

}

Thank you very much for any help

1 Like

The function that you feed to the ROOT::Math::Functor needs to return a double. The one you pasted seems to return void and instead sets the result in some global variable “f”.

Jean-François

I’ve rewritten the function into a function with a return value, but that did not solve my problem.

What i do now is trying to adapt the example from the fitting tutorials “combinedFit.C”
"http://root.cern.ch/root/html/tutorials/fit/combinedFit.C.html"
and use the GSL Multi Fit Algorithm but this is giving me a hard time.
I will explain you what i want to do and hopefully you can help me.

In the example two histograms are being fitted. In my case, for now, i just want to fit datapoints with a fit model which is the following:

double Extrapolation(double x*, double par*){

     return (par[0] + par[1]*par[2]/x[0] + par[1]*par[2]*par[2]*x[1]*x[1]/(x[0]*x[0]))

}

It is a two dimensional function with 3 parameters. My next step would be to fit several of those functions at the same time since all of them share the parameter par[2]. But for the moment i would be fully satisfied to learn how to fit one of those functions by adapting the example quoted above.

So i am filling a TGraph2DErrors with data points. Next i instantiate an object

TF2 *f2 = new TF2(“f2”,Extrapolation,xmin,xmax,ymin,ymax,npar);

and set its parameters. What next? In the example
ROOT::MATH::WrappedMultiTF1
is used but i am afraid there is nothing like
"ROOT::MATH::WrappedMultiTF2".

I heavily appreciate your help since getting this very soon is utterly important for me.

Regards,
Chris

Hi,
If you want to fit using the GSLMultiFit algorithm (the Levenberg-Marquardt algorithm) you can do in two ways:

  • either using the Fit method of the histogram of graph you want to fit.
  • or by providing yourself the chi2 function and using the ROOT::Fit::Fitter interface. In this second case your function must implement the ROOT::Math::FitMethodFunction interface, otherwise it will now work, because that algorithm needs to know all the chi2 contributions.
    See root.cern.ch/root/htmldoc/ROOT__ … iDim_.html

Best Regards

Lorenzo

Hi,

thank you for the answer but i am not sure if i understood since i am quite unexperienced.

I did the following, which does not work:


class MyParamFunc{

        public:

        double Extrapolation(double *x, double *par) {

                double value = (par[0] + par[1]*par[2]/x[0] + par[1]*par[2]*par[2]*x[1]*x[1]/(x[0]*x[0]));

                return value;
        }
};


void NonLinLSFit() { 

  
        int nEntriesFit = 8;

        ostringstream ossFit;
        ossFit << "fit_mass_table_beta_1_qmass_0.01_tc.dat";

        //Function for filling the global vectors
        read_fit_file(ossFit, nEntriesFit);
    
        TGraph2DErrors *gr = new TGraph2DErrors(nEntriesFit);

        //Setting the points for the graph
        for(int i = 0; i < nEntriesFit; i++){

                gr->SetPoint(i,dVolumeFit[i],dChargeFit[i],dMassFit[i]);
                gr->SetPointError(i,0,0,dErrorFit[i]);
        }   
    
        MyParamFunc* fobj = new MyParamFunc();

        int npar = 3;

        TF1 * MyFunc = new    TF1("MyFunc",fobj,&MyParamFunc::Extrapolation,0,1,npar,"MyParamFunction","Extrapolation");

        ROOT::Math::WrappedMultiTF1 EP(*MyFunc,2);    

        ROOT::Fit::DataOptions opt;
        ROOT::Fit::DataRange rangeB;
    
        rangeB.SetRange(0,0.007);

        ROOT::Fit::BinData dataEP(opt,rangeB);
        ROOT::Fit::FillData(dataEP, gr);
    
        ROOT::Fit::Chi2Function EPChi2(dataEP, EP);

        ROOT::Fit::Fitter EPfitter;

        const int EPparams = 3;
        double FitterParams[EPparams] = {0.6,-0.3,0.3}; 

        EPfitter.Config().SetParamsSettings(3,FitterParams);

        EPfitter.Config().ParSettings(0).SetLimits(0.5,1.0);
        EPfitter.Config().ParSettings(1).SetLimits(0.,-0.5);
        EPfitter.Config().ParSettings(2).SetLimits(0.5,-0.5);
    
        EPfitter.Config().ParSettings(0).SetStepSize(0.0001);
        EPfitter.Config().ParSettings(1).SetStepSize(0.0001);
        EPfitter.Config().ParSettings(2).SetStepSize(0.0001);

        EPfitter.Config().MinimizerOptions().SetPrintLevel(1);
        EPfitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");

        EPfitter.FitFCN(EPChi2,0,dataEP.Size(),true);
        ROOT::Fit::FitResult result = EPfitter.Result();

        result.Print(std::cout);

}

I can compile and run the code but the output will be:

Minimize using GSLNLSMinimizer undefined
Error in ROOT::Math::GSLNLSMinimizer::Minimize: Error setting the
residual functions ; iret = -1


       Invalid FitResult            

Minimizer is GSLMultiFit
Chi2 = 0
NDf = 0
Edm = -1
NCalls = 0
Par_0 = 0.6 (limited)
Par_1 = -0.3
Par_2 = 0.3

Could you please have a short look at the code and tell me whether it is totally wrong and maybe give some hints?

Hi,

Yes, I can have a look at your code, but please post it, if possible.

Lorenzo

Hi,

you mean that I should not post it in the Code environment as I did before?

I started the following approach this time.

I oriented myself at the fitting tutorial combinedFit.C which can be found here
http://root.cern.ch/root/html/tutorials/fit/combinedFit.C.html
where two histograms are fitted simultaneously.

My fit model function though is a parametric two dimensional function and if I understood it correctly than the
WrappedMultiTF1 class can handle such functions. What led me to this conclusion is what is written here
http://root.cern.ch/root/html/TF1.html
under
D - A general C++ function object (functor) with parameters

The function I put into the TF1 comes from here
http://root.cern.ch/drupal/content/how-implement-mathematical-function-inside-framework
and is a
IParametricFunctionMultiDim

Finally my code is the following:

[code]#include “Math/IFunction.h”
#include “Math/IParamFunction.h”
#include “Fit/Fitter.h”
#include “Fit/BinData.h”
#include “Fit/Chi2FCN.h”
#include “TH1.h”
#include “TList.h”
#include “Math/WrappedMultiTF1.h”
#include “Math/FitMethodFunction.h”
#include “HFitInterface.h”
#include “TCanvas.h”
#include “TStyle.h”
#include “TGraph2DErrors.h”
#include
#include
#include
#include

using namespace std;

vector dMassFit;
vector dVolumeFit;
vector dChargeFit;
vector dErrorFit;

void read_fit_file(ostringstream& ossFit, int nEntriesFit){

string Value; 
char GetLine[1000];

ifstream file;
file.open(ossFit.str().c_str());

if(file.good() == 0){
	cout << "File does not exist, terminating ..." << endl;
	exit(0);
}

int nCounter = 0;

for(int i = 0; i < 3; i++){
	file.getline(GetLine,1000);
}

while((nCounter < nEntriesFit)){
	file >> Value;
	dMassFit.push_back(atof(Value.c_str()));	
	file >> Value;
	dErrorFit.push_back(atof(Value.c_str()));	
	file >> Value;
	dVolumeFit.push_back(atof(Value.c_str()));	
	file >> Value;
	dChargeFit.push_back(atoi(Value.c_str()));					

	for(int i = 0; i < 5; i++){
	file >> Value;
	}


	nCounter++;
}

file.close();

}

void DisplayFitData(){

for(int i = 0; i < 8; i++){

	cout << dMassFit.at(i) << " " << dVolumeFit.at(i) << " " << dChargeFit.at(i) << " " << dErrorFit.at(i) << endl;
}

}

class MyParametricFunction: public ROOT::Math::IParametricFunctionMultiDim
{
private:
const double* pars;

public:
	double DoEvalPar(const double* x, const double* p) const
	{
		return (p[0] + p[1]*p[2]/x[0] + p[1]*p[2]*p[2]*x[1]*x[1]/(x[0]*x[0]));		
	}

	unsigned int NDim() const
	{
		return 2;
	}

	ROOT::Math::IParametricFunctionMultiDim* Clone() const
	{
		return new MyParametricFunction();
	}

	const double* Parameters() const 
	{
		return pars;
	}

	void SetParameters(const double* p)
	{
		pars = p;
	}

	unsigned int NPar() const
	{
		return 2;
	}

};

void NonLinLSFit() {

int nEntriesFit = 8;

ostringstream ossFit;
ossFit << "fit_mass_table_beta_1_qmass_0.01_tc.dat";

read_fit_file(ossFit, nEntriesFit);

TGraph2DErrors *gr = new TGraph2DErrors(nEntriesFit);

DisplayFitData();

//Setting the points for the graph
for(int i = 0; i < nEntriesFit; i++){

	gr->SetPoint(i,dVolumeFit[i],dChargeFit[i],dMassFit[i]);
	gr->SetPointError(i,0,0,dErrorFit[i]);
}

MyParametricFunction* fobj = new MyParametricFunction();

int npar = 3;
//TF1 * MyFunc = new TF1("MyFunc",fobj,&MyParamFunc::Extrapolation,0,1,npar,"MyParamFunction","Extrapolation");
TF1 * MyFunc = new TF1("MyFunc",fobj,0,1,npar,"MyParametricFunction");

ROOT::Math::WrappedMultiTF1 EP(*MyFunc,2);	

ROOT::Fit::DataOptions opt;
ROOT::Fit::DataRange rangeB;

rangeB.SetRange(0,0.007);

ROOT::Fit::BinData dataEP(opt,rangeB);
ROOT::Fit::FillData(dataEP, gr);

ROOT::Fit::Chi2Function EPChi2(dataEP, EP);

ROOT::Fit::Fitter EPfitter;

const int EPparams = 3;
double FitterParams[EPparams] = {0.6,-0.3,0.3}; 

EPfitter.Config().SetParamsSettings(3,FitterParams);

EPfitter.Config().ParSettings(0).SetLimits(0.5,1.0);
EPfitter.Config().ParSettings(1).SetLimits(0.,-0.5);
EPfitter.Config().ParSettings(2).SetLimits(0.5,-0.5);

EPfitter.Config().ParSettings(0).SetStepSize(0.0001);
EPfitter.Config().ParSettings(1).SetStepSize(0.0001);
EPfitter.Config().ParSettings(2).SetStepSize(0.0001);

EPfitter.Config().MinimizerOptions().SetPrintLevel(1);
EPfitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");

EPfitter.FitFCN(EPChi2,0,dataEP.Size(),true);
ROOT::Fit::FitResult result = EPfitter.Result();

result.Print(std::cout);

}[/code]

I can compile the code but the message I get after calling the function is

Minimize using GSLNLSMinimizer undefined
Error in ROOT::Math::GSLNLSMinimizer::Minimize: Error setting the residual functions ; iret = -1


        Invalid FitResult            

Minimizer is GSLMultiFit
Chi2 = 0
NDf = 0
Edm = -1
NCalls = 0
Par_0 = 0.6 (limited)
Par_1 = -0.3
Par_2 = 0.3

I have also uploaded the Code + the data file
fit_mass_table_beta_1_qmass_0.01_tc.dat (437 Bytes)
NonLinearLeastSquareFit2.C (3.63 KB)

Hi,

Your code seems fine. You have only one thing I do not understand. You are calling:

rangeB.SetRange(0,0.007);

This will remove all the graph points, since the range applies to the first graph coordinate. So you are fitting an empty set. Maybe the message is misleading and I need to fix it

Best Regards

Lorenzo

Hey,

this solved my problem! Thank you so much! Simply commenting it out did already do the job.

Please do not close the thread yet since I will shortly come up with a closely linked question about the simultaneous fitting!

Best Regards,
Chris

Hi,

now I am back with new problems.

I have expanded the my script to the case of simultaneously fitting two equal model functions which have 3 parameters each of which they have 1 in common.

I oriented myself on the example another user has posted here
[url]Simultaneous hist fit; ROOT::Math::IMultiGenFunction trouble

This was of great help since it is my plan to make the script applicable for an arbitrary number of equal model functions of which all share 1 common parameter which is what I tried to do.

I can compile the script and everything works so far just not the fitting.
The error message I get is:

GSLNLSMinimizer: Invalid function set - only Chi2Func supported
Error in ROOT::Math::GSLNLSMinimizer::Minimize: Function has not been set
Error: Empty FitResult !

For me it is not clear why this happens since everything is very similar to the tutorial
http://root.cern.ch/root/html/tutorials/fit/combinedFit.C.html

Could you please have a look on my code one more time? I think you could directly see where I am going wrong.

I attached my code and the two necessary files for fitting.

Best regards,
Chris
fit_mass_table_beta_6_qmass_0.03_tc_WL.dat (561 Bytes)
fit_mass_table_beta_6_qmass_0.03_tc_Pion.dat (565 Bytes)
TestCombFit.C (7.05 KB)

OK,

sorry for multiple posting.

I understood from the example written by Lorenzo which can be found here
[url]Locate like command
that I need a more complex object in order to do simultaneous fitting with GSLMultiFit.

The global chi2 function must be a class with a certain structure. I tried to rewrite the script for my purpose and succeeded partially.

What I did was formulating everything in such a way that in principle the user can feed an arbitrary number of fitting functions. The fitting functions are parametric, two dimensional and of equal structure. They share one parameter.

Everything seems fine and I can also compile the code. It just does not work with GSLMultiFit.
The error message I get is:

“Minimize using GSLNLSMinimizer undefined”

For some reason the script runs with Fumili but in the end the fitting results would be nonsense anyway :frowning:.

I would deeply appreciate it if Lorenzo or one of the other fitting experts around could have a look at my code. Maybe you can determine where I wrong and tell me.
I must get this script running within the next days.

Best regards and thanks a lot in advance,
Chris
fit_mass_table_beta_6_qmass_0.03_tc_WL.dat (561 Bytes)
fit_mass_table_beta_6_qmass_0.03_tc_Pion.dat (565 Bytes)
TestCombFit2.C (10 KB)

Hi,

Your code requires a simple change to run. Attached is the fixed version. However the fit does not work, either GSLMultiFit, or Minuit or Fumili. The problem I think is your function. It would need to be debug and understood why it gives a problem. I see you are fitting the different graphs all with the same fit function. Is it correct ?
In that case it is maybe better to build a combined data set instead of building a combined chi2.
To do this just call many times FillData passing the different Graph objects.

Best Regards

Lorenzo
TestCombFit2.C (10.2 KB)

Hi Lorenzo,

I deeply appreciate that you had a look at my code and fixed it. I need this in order to finish my Master Thesis very quickly.
I saw that you made a change about that g array and yes you are right…I am trying to fit all the graphs with only one function.

Before I come to your suggestion:
As far as I understand you the problem must be somewhere in putting my function into this WrappedMultiTF1 object.
Maybe I could reduce the two dimensional parametric function to a one dimensional function by passing the y-coordinate as a parameter since it takes only values and usually runs from 0 to 2 or 3 at the most.
Then things would be much easier I guess since I would not need a class for my function any more and could use a normal WrappedTF1 object.

Now to you suggestion:
To be honest I have no idea where to start. Is there some similar example in the tutorials where I could learn from?

Thank you and best regards,
Chris

Hi,

I some some further investigations on my problems and figured out where the error was. Now I got things working partially. At least I can do simultaneous fitting. Unfortunately some new strange problem occurred.

I have two version of the same script now which both only differ slightly.
The first script is the more general one for an arbitrary number of fitting models and the second one is the more specialized one for just two fitting models.

I went sure several times that I am
- feeding the same data to both of the scripts.
- setting the same starting parameters + limits.

Hi,

sorry for the post before. I sent it accidentally.
I made some some further investigations on my problems and figured out where the error was. Now I got things working partially. At least I can do simultaneous fitting.

The problem seems to be in my global chi2 function. I will do further investigations on that.

Unfortunately some new strange problem occurred which is by far more important right now since it affects the fitting results.

I have two version of the same script now which both only differ slightly.
The first script is the more general one for an arbitrary number of fitting models and the second one is the more specialized one for just two fitting models.

I went sure several times that I am
- feeding the same data to both of the scripts.
- setting the same starting parameters + limits.
- using the same fitter (GSL) for both scripts.
- using the same global chi2 function for both scripts.
- using the same fitting function for both scripts.

For some reason the fits will give different results. That is totally strange. I get no error message at any time or whatsoever. What is yet strange is that for both fits it says directly before fitting:
“Minimize using GSLNLSMinimizer undefined”
but afterwards it directly starts fitting using the GSLNLSMinimizer.
And I also should note that both scripts are running on the same machine.
Though the fit of the generalized script seems much better which means a lower chi2/dof and the parameters are slightly close to the “real” values.

It would be great if someone could have a look at the scripts and give a suggestion.

I will attach both of them with the two fitted data files. Both just must be compiled and run with no further adjustings.

Thanks and best regards,
Chris

EDIT: Could my be some of the moderators move this thread into the stat and math tool support? I think it rather belongs there!
fit_mass_table_beta_5_qmass_0.03_tc_WL_exp_1.dat (322 Bytes)
fit_mass_table_beta_5_qmass_0.03_tc_Pion.dat (459 Bytes)
GenCombFit.C (13.5 KB)
CombinedTopologyNLLSFit.C (9.36 KB)

Hello,

I just wanted to announce that all my problems with the code are solved. I located to problem in the part of the global chi2 function where the decision is met to which local chi2 function the i-th data point is passed.

In principle the code works for an arbitrary number of equal fit models.
In case some one should be interested i will attach the running code to this post.

I guess this thread can be considered solved. But the thread should still be moved to the stat and math tool support forum.

Lorenzo once again thank you for all the efforts!

Best regards,
Chris
fit_mass_table_beta_5_qmass_0.03_tc_WL_exp_1.dat (322 Bytes)
fit_mass_table_beta_5_qmass_0.03_tc_Pion.dat (459 Bytes)
GenCombFit.C (11.2 KB)

Hi,

here is a little update of the script with a few bug fixes and improvements.

The files posted are:

(1) The fitting script
(2) A parameter file where the initial parameters must be filled in in the specified format
(3) + (4) two sample data files for testing the script

In the beginning of the main function NonLinLSFit() the observable vector must be specified.
Furthermore the parameter file must be edited. That’s all.

The script allows for fitting an arbitrary number of observables with the same fitting model.
The fitting works equally well with GSLMultiFit, Minuit2 and Fumili.

This thread can considered solved.

Best regards,
Chris
NonLinFittingScript.zip (5.84 KB)

Dear ROOTers,

I have tried to run this example code on ROOT (5.34/01) on MS Windows (10.0.15063). For who has/wants to do the same, there are some tips I would like to spot out:

A normal loading and calling procedure will not work on a standar ROOT session called from cmd.exe, this is because some dynamic libraries have to be built before. Doing the altter will lead to a crash. For this reason first we build the libraries from the prompt of VS.
Before starting, it is CRUCIAL to have installed the EXACT SAME version of Visual Studio, with which ROOT has been built. I stress that this is CRUCIAL!
In the case of VS2010, the procedure to set up its prompt, can be found at this link.

When this is done, the rest is quite simple:

  1. open the Command Prompt of VS
  2. load the macro from your specific path adding “+” at the end, in order to have access to ACLiC (example .L GenCombFit.C+)
  3. Run the function (i.e. NonLinLSFit() ). NOTE: now it is possible to load the macro and call the fucntion from a standard comamdn prompt cmd.exe

Hope that works for everybody.

Cheers,

Giammarco