Dear experts,
I want to fit some datasets to the output of a numerical model, where 5 parameters are left free.
The problem is that this numerical model is very slow, even running it in parallel on a GPU or on multiple CPUs (>= 5 minutes), so it’s not feasible to call it from within a TF1.
In addition, I have multiple datasets to fit, so it’s not efficient to run the model with the same combination of parameters more than once, like it would happen with a TF1.
So I created a 5D grid spanning reasonable ranges in the parameters of interest, run the model for all the grid points (order of 1e6) and for each datasets I computed the chi2 for every grid points and saved them into separate TTrees.
Then I created a class to provide a FCN for ROOT::Fit::Fitter, by linearly interpolating the chi2 values from the 5D grid:
class GridChi2
{
public:
GridChi2(TTree *t_chi2, ModelsCollection *mc) :
t_chi2(t_chi2), mc(mc), npars(5)
{
t_chi2->SetBranchAddress("chi2", &chi2);
}
// 5D linear interpolation of chi2 from the grid saved in t_chi2
Double_t operator()(const Double_t *par)
{
UShort_t ipars[npars];
Double_t s[npars];
mc->GetLinearInterpolatingFactors(par[0], par[1], par[2], par[3], par[4], ipars, s);
Double_t chi2 = 0.;
for (UShort_t n = 0; n < 1<<npars; ++n)
{
UShort_t shift[npars] = { 0 };
Double_t k = 1.;
for (UShort_t ipar = 0; ipar < npars; ++ipar)
{
shift[ipar] = (n >> ipar) & 1;
if (shift[ipar]) k*= s[ipar];
else k *= 1 - s[ipar];
ipars[ipar] += shift[ipar];
}
Long64_t imodel = mc->iModel(ipars);
t_chi2->GetEntry(imodel);
chi2 += k*model_chi2;
for (UShort_t ipar = 0; ipar < npars; ++ipar)
{
ipars[ipar] -= shift[ipar];
}
}
return chi2;
}
private:
TTree *t_chi2;
ModelsCollection *mc;
Float_t model_chi2;
const UShort_t npars;
};
ModelsCollection is a class which contains information about the grid and the models output: its code is not important here.
I tested the multidimensional linear interpolating code and it does what it’s supposed to do.
Then I fit the data in this way:
ROOT::Fit::Fitter fitter;
fitter.Config().MinimizerOptions().SetMaxIterations(5e4);
fitter.Config().MinimizerOptions().SetMaxFunctionCalls(1e5);
fitter.Config().MinimizerOptions().SetTolerance(1);
fitter.Config().MinimizerOptions().SetPrintLevel(3);
fitter.Config().MinimizerOptions().SetStrategy(2);
fitter.Config().ParamsSettings().resize(npars);
for (UShort_t ipar = 0; ipar < npars; ++ipar)
{
fitter.Config().ParSettings(ipar).Set(par_name[ipar], par_best[ipar], grid_step[ipar]/2., max(0., par_min[ipar] - grid_step[ipar]), par_max[ipar] + grid_step[ipar]);
}
GridChi2 grid_chi2(t_chi2, &mc);
fitter.FitFCN(npars, grid_chi2, 0, dof, true);
where par_best is an array with the values of the parameters for which the chi2 in the grid is minimum.
I tried all combinations of minimizer and algorithm, but none of them is working.
Migrad, Minimizer and MigradImproved fail with status = 4 (in reality because the Hessian matrix is not positively defined), and so does Minos.
Simplex and Scan converge, but they do not move much far from the starting point (even if the starting point is very far from the minimum of the grid…) and the errors are too small to make sense.
All others algorithms (Seek, GSL*, Genetic) either do not converge or seems to run indefinitely.
For now, I’m doing a sort of manual fitting: for each parameter, I marginalize the chi2 vs the parameter, then compute the parabola passing through the lowest chi2 and its neighbor points (one left, one right of the minimum). This way I can extract for each parameter a “best fit” value and the classic 1 sigma error by imposing chi2_err = chi2_min + 1.
The problem is that sometimes:
- for some parameter the marginalized chi2 curve is not parabolic at all, so the minimum is ill-defined.
- each parameter has a different minimum chi2, because the parabola is different for each parameter: so what is the actual chi2 of the global fit?
- the minimum of the parabola can become negative if the marginalized chi2 curve is very narrow: this of course does not make any sense and will give a smaller error than the actual one.
Essentially, I’m reimplementing a very simple version of Minos, with no covariance matrix, nor 2D contours.
That’s why I wanted to use ROOT::Fit::Fitter…
Any idea if my GridChi2 class and the way I use it make sense?
Or how to solve the problems in my manual fitting method?
Thanks.
Claudio