Hello. With 6.10/08, today I typed in the terminal:
f = new TF1("f","gaus(0)",-20,20)
f->SetParameters(0.13,5.4428,0.01)
f->Mean(-20,20)
(double) 5.632071 //!!
f->Mean(-10,10)
Error in <TF1::Moment>: Integral zero over range
(double) 0.000000
f->SetNpx(1e6) //Hoped this would increase accuracy of TF1::Integral
f->Mean(-20,20,0,1e-100) //Decreasing "epsilon" argument
(double) 5.632071 //The same
f->Variance(-20,20,0,1e-100)
(double) 0.000000
With f->SetParameter(2,0.05) (less peaked) I get the correct numbers. Values in between (e.g. f->SetParameter(2,0.029)) can give very strange results. Since the way TF1::Moment works, this probably has to do with TF1::Integral.
In the end, as a workaround I used h = f->GetHistogram().
So:
How can I instruct ROOT to calculate the Integral better in order to get an estimation of mean and variance of a TF1 comparable to the one obtained by degrading to an histogram? Or is the latter the best route?
Is there some sign/flag/rule of thumb I can check to understand if the result may not be accurate? The TF1 I was actually looking at was more complicated, and I didn’t realize immediately the problem.
I see. I suspected someone else already had to ask this, but I wasn’t able to find it by searching (we all start from different problems, it seems). Thank you for the references! (the big green grinning face as preview is lovely)
In particular I copypaste for future readers:
( Altough they all gave the same 5.632071 for my f->Mean(-20,20) )
And:
ROOT’s “integrators” are well known to misbehave in all cases when the function “changes rapidly” or if there are “narrow / sharp peaks” inside. I think that, in vast majority of cases like yours, people “blindly” use the returned “integral” without noticing that it is wrong. What’s even worse, “integrals” are quietly calculated and used internally by some methods in various classes (users may not even realize that this happens).
(If I may ask, would you guys mind write this in big letters in the TF1 documentation? By now, in TF1::Integral it says “IntegralOneDim or analytical integral”, so for a gaussian I would assume it always gives the exact result!)
So, I interpret your answer as: the method giving the most accurate results is changing in my code every TF1::Mean etc. with
f->SetNpx( (ULong_t) (1e3*(xmax-xmin)) ); //Big Number Here
h = f->GetHistogram();
h->GetMean();
h->GetStdDev();
h->Integral("width");
It should be ok to calculate moments of any order (by building a new TF1 = x^n f(x) ) this way.