Calculate the derivative with respect to each parameter

ROOT
, ,
1

Dear Experts,
I have a function y = f(x,a,b) with known a and b parameters. I would like to calculate the uncertainty of the function. To do that I have to calculate the derivative of f as a function of a and b. In this case, dy should be: dy = sqrt(((df/da)^2 * da^2 + (df/db)^2 * db^2 + (df/da)(df/db)dadb).

I have the covariance matrix of da^2, db^2 and dadb as:

conv= [a2, ab, b2 ] = [0.0009778, 0.009128, 1.016]

I couldn’t find any example showing how I can do that. I found some example suggest to use:

void GradientPar(const Double_t *x, Double_t *grad, Double_t eps)

But I don’t understand how I should use it and how I can provide the covariance matrix

Can you help me by showing some examples

Here is my function:

float * GEnFunction(float GEn_arr[], int length){
    Double_t a = 3.1320, b = 2.97626;
    Double_t mass=0.93957, lambda2=0.71;
    for (int i = 0; i < length; ++i){
        Double_t tau = q2_bin[i]/(4*pow(mass,2.0));
        Double_t GD  = 1.0/pow(1 + (q2_bin[i]/lambda2),2.0);
        GEn_arr[i] =(a*tau/(1 + (b*tau)))/GD;

        
    }
   return GEn_arr;
}

and the Covariance matrix:

conv= [a2, ab, b2 ] = [0.0009778, 0.009128, 1.016]

Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided

2

Define a TF1 object with your function and have a look at this ROOT Forum question.

3

Thanks @Eddy_Offermann for your reply by following the link you mentioned, I defined TF1 object and follow what @moneta explained. I just have issued of the way that I can provide my own covariance matrix, conv= [a2, ab, b2 ] = [0.0009778, 0.009128, 1.016], in the GradientPar ?
Can help me to do that please @moneta @Eddy_Offermann

Below is the code I use:

#include "TF1.h"
#include "TH1D.h"
#include "TFitResult.h"
#include "TMath.h"
#include <assert.h>
#include <iostream>
#include <cmath>
 
TF1 * fitFunc;  // fit function pointer

Double_t GEnFunction(double Q2) {
  
  Double_t mass=0.939565, lambda2=0.71;
  Double_t a = 3.1320, b = 2.97626;
  //  calculate GEp with parameters from Kelly fit.
  Double_t tau = Q2/(4*pow(mass,2.0));
  Double_t GD  = 1.0/pow(1 + (Q2/lambda2),2.0);
  double GEp = (a*tau/(1 + (b*tau)));

  return GEp;
}

void formfunc() {
   // example of script showing how to divide a canvas
   // into adjacent subpads + axis labels on the top and right side
   // of the pads.
   //Author; Rene Brun
   gROOT->Reset();
   auto myGEp = new TF1("myGEp","GEnFunction(x)",1,5);
   Double_t derivate[2];
   Double_t x[1];
   x[0] = 5.;
   myGEp->FixParameter(0, 3.1320); 
   myGEp->FixParameter(1, 2.97626); 
   
   myGEp->GradientPar(x, derivate);

   std::cout << derivate[0] << std::endl;

   // conv= [a2, ab, b2 ] = [0.0009778, 0.009128, 1.016]


   myGEp->Draw(); // draw it.
}
4

So you got now from the call to GradientPar the values for df/da and df/db. Insert those and your covariance values in the above formula.

5

It seems to me that, assuming covariance matrix cov(a, b) = {{va, cab}, {cab, vb}}, it should be:
dy = sqrt( (df/da)^2 * va + (df/db)^2 * vb + 2 * (df/da) * (df/db) * cab )

6

By following your way, the result I got is 0, while the result I should get is 1.9203e-3.
The df/da result I got is 2.30002e-314 and df/db is 6.95159e-310.
I am not quite sure where is the mistake I have done in my code
@Eddy_Offermann @Wile_E_Coyote
Here is my code

#include "TF1.h"
#include "TH1D.h"
#include "TFitResult.h"
#include "TMath.h"
#include <assert.h>
#include <iostream>
#include <cmath>
 
TF1 * fitFunc;  // fit function pointer

Double_t GEnFunction(double Q2) {
  
  Double_t mass=0.939565, lambda2=0.71;
  Double_t a = 3.1320, b = 2.97626;
  //  calculate GEp with parameters from Kelly fit.
  Double_t tau = Q2/(4*pow(mass,2.0));
  Double_t GD  = 1.0/pow(1 + (Q2/lambda2),2.0);
  double GEp = (a*tau/(1 + (b*tau)));

  return GEp;
}

void formfunc() {
   // example of script showing how to divide a canvas
   // into adjacent subpads + axis labels on the top and right side
   // of the pads.
   //Author; Rene Brun
   gROOT->Reset();
   auto myGEp = new TF1("myGEp","GEnFunction(x)",1,5);
   myGEp->FixParameter(0, 3.1320); 
   myGEp->FixParameter(1, 2.97626);

   Double_t derivate[2];
   Double_t x[1];
   x[0] = 2.;
   
   myGEp->GradientPar(x, derivate);

   
   Double_t dfa = derivate[0];
   Double_t dfb = derivate[1];
   // conv= [a2, ab, b2 ] = [0.0009778, 0.009128, 1.016]
   Double_t dy = sqrt( (dfa*dfa) * 0.0009778 + (dfb*dfb) * 1.016 + 2 * (dfa*dfb) * 0.009128 );
   
   std::cout << derivate[0] << std::endl;
   std::cout << derivate[1] << std::endl;
   std::cout << dy << std::endl;


   // myGEp->Draw(); // draw it.
}
7

myGEp->Print(); // see Npar

8

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