Fitting with Bernstein polynomials

anikiten · March 10, 2024, 3:31pm

Dear experts,

I am trying to fit range with Bernstein polynomials:

 // background model
  RooRealVar a0("a0","coefficient a0",0.0,1.0);
  RooRealVar a1("a1","coefficient a1",0.0,1.0);
  RooRealVar a2("a2","coefficient a2",0.0,1.0);
  RooRealVar a3("a3","coefficient a3",0.0,1.0);
  RooRealVar a4("a4","coefficient a4",0.0,1.0);
  //
  RooBernstein bkg("bkg","background",dimuonmass2sideband,RooArgList(a0,a1,a2,a3,a4));
  //
  RooRealVar norm("norm","bkg normalization",1700.,0.0,3000.0);
  RooAddPdf model("model","model",RooArgList(bkg),RooArgList(norm));

  dimuonmass2sideband.setRange("R1",12.,25.);
  dimuonmass2sideband.setRange("R2",30.,70.);

  RooFitResult* r = model.fitTo(data,Range("R1,R2"),Save());

and I got a stupid result. All coefficients are set to 0.5 and curve does not describe my data at all

Number of data events = 1601
-log(L) at minimum = -3880.96
chi2 = 9.41812

RooRealVar:: a0 = 0.499798 +/- 0.000166767
RooRealVar:: a1 = 0.500182 +/- 9.49789e-09
RooRealVar:: a2 = 0.499923 +/- 3.1032e-06
RooRealVar:: a3 = 0.500068 +/- 7.51931e-07
RooRealVar:: a4 = 0.499655 +/- 9.76319e-06
RooRealVar:: norm = 1638.55 +/- 0.313795
norm = 1638.55 +/- 0.313795

Can you help, please.

Thanks, Sasha.

Danilo · March 10, 2024, 4:53pm

Hi Sasha,

I am sorry the fit does not seem to converge.
Before we involve the experts, can you verify the messages from the minimizer? Typically something about a fit gone wrong can be found there…

Best,
Danilo

anikiten · March 10, 2024, 6:26pm

Hi Danilo,

I got a lot of following messages:

[#0] WARNING:Minization – RooMinimizerFcn: Minimized function has error status.
Returning maximum FCN so far (16650) to force MIGRAD to back out of this region. Error log follows
Parameter values: a0=0.499798, a1=0.50019, a2=0.499923, a3=0.500068, a4=0.499655, norm=1638.55
RooBernstein::bkg[ x=dimuonmass2sideband coefList=(a0,a1,a2,a3,a4) ]
p.d.f value is less than zero (-0.000006), forcing value to zero @ x=dimuonmass2sideband=70, coefList=(a0 = 0.499798 +/- 3.07965e-07,a1 = 0.50019 +/- 8.79901e-09,a2 = 0.499923 +/- 2.74613e-08,a3 = 0.500068 +/- 2.95359e-08,a4 = 0.499655 +/- 8.64173e-08)

Cheers, Sasha.

anikiten · March 11, 2024, 8:02am

Hi again Danilo,

if you or they can try it I have my .C macro at

/afs/cern.ch/work/a/anikiten/public/CMSSW_11_0_4/src/HiggsAnalysis/CombinedLimit/test/13TeV_2017UL_RooFitTest_SR3_BERNSTEIN_FTEST.C

and root file with data used by this macro is at the same dir.

you do in interactive root

.L 13TeV_2017UL_RooFitTest_SR3_BERNSTEIN_FTEST.C
Draw()

Thanks, Sasha.

Danilo · March 11, 2024, 8:58am

Hi Sasha,

There is something going wrong with your fit. The fact that the pdf has a negative value is a sign. This does not seem to be a ROOT issue.
Are you sure this is the model you need to fit in that particular range? Is it maybe possible to facilitate convergence by more sensible parameters values?
I hope you manage to solve your issue soon!

Cheers,
Danilo

anikiten · March 11, 2024, 11:32am

Hi Danilo,

I wrote in ROOT forum that there was a problem to fit
with Bernstein the subrange of data…

Also, I just setup parameters to be varied in a range of 0-1.
I have no idea in advance what they should be. It gives
for all of them 0.5 and I see just straight line as a fit result
which have nothing to do with my data. I think it is something
wrong in implementation…

Cheers, Sasha.

jonas · March 11, 2024, 12:00pm

Hi @anikiten!

With ROOT 6.30/02, I don’t get any problems running your script:

   ------------------------------------------------------------------
  | Welcome to ROOT 6.30/02                        https://root.cern |
  | (c) 1995-2023, The ROOT Team; conception: R. Brun, F. Rademakers |
  | Built for linuxx8664gcc on Jan 25 2024, 20:04:04                 |
  | From heads/master@tags/v6-30-02                                  |
  | With c++ (GCC) 13.2.1 20230801                                   |
  | Try '.help'/'.?', '.demo', '.license', '.credits', '.quit'/'.q'  |
   ------------------------------------------------------------------

root [0] .L 13TeV_2017UL_RooFitTest_SR3_BERNSTEIN_FTEST.C
root [1] Draw()
[#1] INFO:DataHandling -- RooTreeDataStore::loadValues(data) Skipping event #3 because dimuonmass2sideband cannot accommodate the value 91.5787
[#1] INFO:DataHandling -- RooTreeDataStore::loadValues(data) Skipping event #8 because dimuonmass2sideband cannot accommodate the value 90.2896
[#1] INFO:DataHandling -- RooTreeDataStore::loadValues(data) Skipping event #12 because dimuonmass2sideband cannot accommodate the value 87.1738
[#1] INFO:DataHandling -- RooTreeDataStore::loadValues(data) Skipping event #14 because dimuonmass2sideband cannot accommodate the value 92.4627
[#1] INFO:DataHandling -- RooTreeDataStore::loadValues(data) Skipping ...
[#0] WARNING:DataHandling -- RooTreeDataStore::loadValues(data) Ignored 665 out-of-range events
[#1] INFO:Eval -- RooRealVar::setRange(dimuonmass2sideband) new range named 'R1' created with bounds [12,25]
[#1] INFO:Eval -- RooRealVar::setRange(dimuonmass2sideband) new range named 'R2' created with bounds [30,70]
[#1] INFO:Eval -- RooRealVar::setRange(dimuonmass2sideband) new range named 'Rtotal' created with bounds [12,70]
[#1] INFO:Eval -- RooRealVar::setRange(dimuonmass2sideband) new range named 'signalrange' created with bounds [25,30]
[#1] INFO:Minimization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Eval -- RooRealVar::setRange(dimuonmass2sideband) new range named 'fit_nll_model_data_R1' created with bounds [12,25]
[#1] INFO:Eval -- RooRealVar::setRange(dimuonmass2sideband) new range named 'fit_nll_model_data_R2' created with bounds [30,70]
[#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_model_data) constructing test statistic for sub-range named R1,R2
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (bkg)
Minuit2Minimizer: Minimize with max-calls 3000 convergence for edm < 1 strategy 1
Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator
Info in <Minuit2>: MnSeedGenerator Initial state: FCN =      -3854.634496 Edm =       4162.035796 NCalls =     25
Info in <Minuit2>: MnSeedGenerator Initial state
  Minimum value : -3854.634496
  Edm           : 4162.035796
  Internal parameters:	[                0                0                0                0                0     0.1337315894]	
  Internal gradient  :	[     -112.0551418     -95.80473469     -30.34783612      68.56709204      169.6180087     -41.58342261]	
  Internal covariance matrix:
[[    0.010923121              0              0              0              0              0]
 [              0    0.028469817              0              0              0              0]
 [              0              0    0.047980141              0              0              0]
 [              0              0              0     0.22256729              0              0]
 [              0              0              0              0     0.52680383              0]
 [              0              0              0              0              0   0.0016261633]]]
Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.001 with call limit = 3000
Info in <Minuit2>: VariableMetricBuilder    0 - FCN =      -3854.634496 Edm =       4162.035796 NCalls =     25
Info in <Minuit2>: VariableMetricBuilder    1 - FCN =       -4011.84414 Edm =       167.1999919 NCalls =     56
Info in <Minuit2>: VariableMetricBuilder    2 - FCN =      -4014.242847 Edm =        14.1727963 NCalls =     78
Info in <Minuit2>: VariableMetricBuilder    3 - FCN =      -4019.787684 Edm =      0.5519667999 NCalls =     99
Info in <Minuit2>: VariableMetricBuilder    4 - FCN =      -4020.241236 Edm =       0.267753583 NCalls =    113
Info in <Minuit2>: VariableMetricBuilder    5 - FCN =       -4020.77624 Edm =      0.1586104724 NCalls =    129
Info in <Minuit2>: VariableMetricBuilder    6 - FCN =       -4021.02251 Edm =     0.02813089728 NCalls =    145
Info in <Minuit2>: VariableMetricBuilder    7 - FCN =      -4021.090976 Edm =     0.02030282874 NCalls =    162
Info in <Minuit2>: VariableMetricBuilder    8 - FCN =       -4021.15916 Edm =     0.06146040058 NCalls =    178
Info in <Minuit2>: VariableMetricBuilder    9 - FCN =      -4021.248087 Edm =     0.04570298805 NCalls =    195
Info in <Minuit2>: VariableMetricBuilder   10 - FCN =      -4021.434175 Edm =      0.1403543632 NCalls =    219
Info in <Minuit2>: VariableMetricBuilder   11 - FCN =      -4021.491854 Edm =       0.154922566 NCalls =    236
Info in <Minuit2>: VariableMetricBuilder   12 - FCN =       -4021.58361 Edm =     0.08267535242 NCalls =    255
Info in <Minuit2>: VariableMetricBuilder   13 - FCN =      -4021.729821 Edm =      0.5523697551 NCalls =    274
Info in <Minuit2>: VariableMetricBuilder   14 - FCN =      -4021.940602 Edm =      0.3002726429 NCalls =    292
Info in <Minuit2>: VariableMetricBuilder   15 - FCN =      -4022.138749 Edm =      0.1547691695 NCalls =    310
Info in <Minuit2>: VariableMetricBuilder   16 - FCN =      -4022.258338 Edm =      0.0159670003 NCalls =    326
Info in <Minuit2>: VariableMetricBuilder   17 - FCN =      -4022.275158 Edm =   0.0001518387563 NCalls =    340
Warning in <Minuit2>: MnPosDef Matrix forced pos-def by adding to diagonal 0.0645936
Info in <Minuit2>: VariableMetricBuilder After Hessian
Info in <Minuit2>: VariableMetricBuilder   18 - FCN =      -4022.275158 Edm =   0.0001996152154 NCalls =    382
Warning in <Minuit2>: Minuit2Minimizer::Minimize Covar was made pos def
Minuit2Minimizer : Valid minimum - status = 1
FVAL  = -4022.27515826609215
Edm   = 0.000199615215363708223
Nfcn  = 382
a0	  = 0.431767	 +/-  0.39518	(limited)
a1	  = 0.990216	 +/-  0.725769	(limited)
a2	  = 0.481864	 +/-  0.324117	(limited)
a3	  = 0.136878	 +/-  0.268015	(limited)
a4	  = 0.135699	 +/-  0.146703	(limited)
norm	  = 1840.96	 +/-  47.9275	(limited)
Info in <Minuit2>: Minuit2Minimizer::Hesse Using max-calls 3000
Info in <Minuit2>: Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) p.d.f was fitted in a subrange and no explicit NormRange() was specified. Plotting / normalising in fit range. To override, do one of the following
	- Clear the automatic fit range attribute: <pdf>.removeStringAttribute("fitrange");
	- Explicitly specify the plotting range: Range("<rangeName>").
	- Explicitly specify where to compute the normalisation: NormRange("<rangeName>").
	The default (full) range can be denoted with Range("") / NormRange("").
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) only plotting range 'R1,R2'
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) p.d.f. curve is normalized using explicit choice of ranges 'fit_nll_model_data_R1,fit_nll_model_data_R2'
  Number of data events = 1601
-log(L) at minimum = -4022.28
chi2 = 1.03646
  1) RooRealVar::   a0 = 0.431767 +/- 0.244056
  2) RooRealVar::   a1 = 0.990216 +/- 0.649953
  3) RooRealVar::   a2 = 0.481864 +/- 0.322393
  4) RooRealVar::   a3 = 0.136878 +/- 0.11275
  5) RooRealVar::   a4 = 0.135699 +/- 0.0836204
  6) RooRealVar:: norm = 1840.96 +/- 46.664
  norm = 1840.96 +/- 46.664
Info in <TCanvas::Print>: pdf file 2017UL_SR3_CHEB_FTEST.pdf has been created

Is that the output that you expect? Maybe it is fixed by circumventing the multi-range fit problem I mentioned in the other thread, either by doing an extended fit or by using ROOT 6.26 or newer.

Cheers,
Jonas

anikiten · March 11, 2024, 12:14pm

ah, thank you very much, Jonas !

I suspected it might be the reason. I will move then to
ROOT 6.18/04.

Thanks again, cheers, Sasha.

jonas · March 11, 2024, 12:21pm

ROOT 6.18/04 is not new enough, it needs to be at least 6.26.00

anikiten · March 11, 2024, 12:36pm

Hi again Jonas,

sorry for stupid question. How to set ROOT 6.30/02. I am
working with CMSSW_11_3_4 setup which use ROOT 6.22/09
unfortunately.

Thanks, Sasha.

jonas · March 11, 2024, 1:02pm

Is it important for you to stick with the CMSSW environment? If yes, you can just sources a newer one, like CMSSW_14_0_1 from lxplus8 or lxplus9.

Or just log onto lxplus9 without any cmsenv, because ROOT 6.30 is the default there.

Cheers,
Jonas

anikiten · March 11, 2024, 2:13pm

Thanks, Jonas,

source /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.30.02/x86_64-almalinux8.8-gcc85-opt/bin/thisroot.csh

also works. I got my Bernstein running !

Cheers, Sasha.

system · March 25, 2024, 2:14pm

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.