Extracting mean after fitting a distribution

Newbie
1

I am able a fit a distribution using a fitting function

"[1]*(exp(0.5*(log(2)/[0])*((log(2)/[0])*[2]**2-2*x))*TMath::Erfc(((log(2)/[0])*[2]**2-x)/(sqrt(2)*[2])))"

Now I am interested to get the mean of this distribution through fitting.

Here is the code :

   TF1 *f1 = new TF1("f1", "[1]*(exp(0.5*(log(2)/[0])*((log(2)/[0])*[2]**2-2*x))*TMath::Erfc(((log(2)/[0])*[2]**2-x)/(sqrt(2)*[2])))", 0.005, 100);
   f1->SetParameter(0, 6);
   f1->SetParameter(1, 25);
   f1->SetParameter(2, 10);

   f1->GetParameter(0);
   f1->GetParameter(1);
   f1->GetParameter(2);

How to obtain the mean in this fitting?

2

Hi @iamaanand,

thank you for the question. For the general, beginner level info on fitting please refer to: Fitting histograms - ROOT

Cheers,
Marta

3

Hi,
you could use the definition of mean.

Below you find a code to do it with a gaussian, you can change the gaussian with your function.

    double dx=1e-2;
    double x_min = 0.005;
    double x_max = 100;
    
    
    TF1 *f1= new TF1("f","TMath::Gaus(x,19.73, 3.5)",x_min, x_max);
    
    double sum=0;
    double norm=0;
    
    
    
    
    for(double x=x_min; x<=x_max; x+=dx)
    {
        sum+= x*f1->Eval(x)*dx;
        norm+= f1->Eval(x)*dx;
    }
    
    cout<<"mean is  "<< sum/norm<<endl;
1 Like
4

You can also use TF1::Mean(), it gives the same result as @Dilicus code:

void getmean(){
   double dx=1e-2;
   double x_min = 0.005;
   double x_max = 100;

   TF1 *f1= new TF1("f","TMath::Gaus(x,19.73, 3.5)",x_min, x_max);
   double sum=0;
   double norm=0;

   for (double x=x_min; x<=x_max; x+=dx) {
      sum+= x*f1->Eval(x)*dx;
      norm+= f1->Eval(x)*dx;
   }

   cout << "1) mean is  " << sum/norm<<endl;
   cout << "2) mean is  " << f1->Mean(x_min,x_max) << endl;
}

root [0] .x getmean.C
1) mean is  19.73
2) mean is  19.73
root [1]
2 Likes
5

I looked on TF1 documentation searching for GetMean() instead of Mean(), and since I did not found it, I wrote a a small code.
I just realised that I got confused with TH1.
Of course @couet solution is way better than mine in terms of readability and simplicity.

1 Like
6

I did the same and found Mean finally :grinning:

It is good to crosscheck anyway :+1:t2:

1 Like
7

As far as I know, in ROOT, once you’ve fit your data with a function (like the one you’ve described), the mean of the distribution can be obtained, but it requires some understanding of the fitted function. The mean of a distribution is not always directly given by the fit parameters, especially for complex functions. Your function looks quite specific and might not have a straightforward mean. But correct me if I’m wrong.

8

I know the topic’s been solved but I just wanted to suggest using numerical integration for your unique function to find the mean. This method is really handy for complex functions where the mean isn’t obvious.
I think it might look something like this:

double x_min = 0.005;
double x_max = 100;
TF1 *f1 = new TF1("f1", yourFunction, x_min, x_max);

double mean = f1->Integral(x_min, x_max, [](double *x, double *p) { return x[0] * f1->EvalPar(x, p); }) / f1->Integral(x_min, x_max);
cout << "Mean: " << mean << endl;

This code integrates your function multiplied by x (for the weighted average) over your range and normalizes it. It’s a great way to handle non-standard distributions.

What do you all think of this? Might it be a good solution?

9

Can you provide a running script ? The code you posted does not work if I put it my previous script.

10

In essence, the idea is to use numerical integration to calculate the weighted mean of a complex function, which is particularly useful when dealing with non-standard distributions. The code I provided demonstrates how to integrate the function multiplied by ‘x’ to find the weighted average, and then normalizes it over a specified range. This approach should be adaptable to a variety of functions, but the exact implementation might vary depending on the details of your function and script.

11

You provided a code snippet not a working script. Can you provide a working script please ?