Vavilov distribution

Roofit and RooStats
1

Dear ROOTers,

I would like to implement the Vavilov distribution in MathMore.
In CERNLIB there are two implementations: G115 wwwasdoc.web.cern.ch/wwwasdoc/sh … 6/top.html (VVIDEN etc), where G115 is supposedly faster, and G116 is more accurate.

Curently, G115 is implemented in TMath::Vavilov and TMath::VavilovI,
G116 is unavailable in ROOT (as far as I know).

Both methods need a relatively expensive setup routine to be called before the function can be evaluated. In TMath::Vavilov, the setup is repeated every time, which is probably not optimal in terms of speed. I think that a class is more suited to this problem, where the object holds the setup information for a specific combination of kappa and beta.

I propose the following:

  • Define a base class Vavilov, with two subclasses VavilovFast (G115) and VavilovAccurate (G116).
  • Vavilov has these methods:
    • set (kappa, beta2) // set or change kappa and beta2
    • pdf (x) //evaluates the pdf
    • pdf (x, kappa, beta2) // if necessary, calls set, then evaluates the pdf
    • cdf (x)
    • cdf (x, kappa, beta2)
    • quantile (z) // only for VavilovFast
    • quantile (z, kappa, beta2)

For many applications, free functions

  • vavilov_pdf
  • vavilov_cdf
  • vavilov_quantile
    are probably needed. I propose to add a static “instance” pointer to the Vavilov base class, plus something like “SetPolicy” (either kFast or kAccurate), and implement the free functions like this:

double vavilov_pdf (x, kappa, beta2) {
Vavilov *v = Vavilov::GetInstance();
return v->pdf (x, kappa, beta);
}

I have already written a first version of VavilovAccurate, by the way.

Any comments or suggestions?

Cheers, Benno

2

Hello Benno,

I agree with you , since the functions are expensive to compute due to a along set up time it would be good to have a class in order to evaluate quickly the function for the same parameter values. I agree with you also on the class method names. I would probably make it implement (for the pdf evaluation) the ROOT::Math::IParametricFunctionOneDim interface
(project-mathlibs.web.cern.ch/pro … neDim.html) to use easily for fitting and to create from it a TF1.

We could have also the free function where the set up is done in this case every time. The static instance is OK, but probably not needed, the time will be spent anyway in the Set routine and not in the class constructor. I would just do:

double vavilov_pdf (x, kappa, beta2) { 
VavilovAccurate v; 
return v.pdf (x, kappa, beta); 
}

I would add probably two different free function types,
vavilov_pdf (plus vavilob_cdf and vavilov_quantile) and
vavilov_approx_pdf (or vavilov_fast_pdf) instead of using SetPolicy.

Best Regards

Lorenzo

3

Hello Lorenzo,

thank you for the comments!
I will implement it along these lines.

Cheers, Benno

4

Hi Lorenzo,

the Vavilov classes are almost finished, including tests.
One issue remains:

[quote=“moneta”]Hello Benno,

I agree with you , since the functions are expensive to compute due to a along set up time it would be good to have a class in order to evaluate quickly the function for the same parameter values. I agree with you also on the class method names. I would probably make it implement (for the pdf evaluation) the ROOT::Math::IParametricFunctionOneDim interface
(project-mathlibs.web.cern.ch/pro … neDim.html) to use easily for fitting and to create from it a TF1.

Lorenzo[/quote]

I looked at the documentation of ROOT::Math::IParametricFunctionOneDim. The problem is that operator() is const, while changing the parameters kappa and beta2 triggers quite expensive calculations. Of course, I can make all internal variables mutable, although I don’t like that too much. Currently, I have two member functions

  double pdf (double x) const;
  double pdf (double x, double kappa, double beta2); // not const!

The alternative would be a separate class, something like

class ROOT::Math::VavilovPdf: public ROOT::Math::IParametricFunctionOneDim {
  public: 
    VavilovPdf (Vavilov& v_) : v(v_) {}
    VavilovPdf () : v(*new VavilovAccurate (1,0)) {}
    double operator() (double x, const double *p) const 
        { return v.pdf (x, p[0], p[1]); }
    double operator() (const double *x, const double *p) const 
        { return v.pdf (x[0], p[0], p[1]); }
  private:
    mutable Vavilov& v;
};

Somehow I like that better, and one could also add other adaptor classes like VavilovCdf.
Maybe operator() should also have free parameters for normalization, scaling and shifts:

double ROOT::Math::VavilovPdf::operator() (double x, const double *p) const { 
  // p[0]: norm, p[1]: x0, p[2]: width, p[3]: kappa, p[4]: beta2
  return p[0]/p[2]*v.pdf ((x-p[1])/p[2], p[3], p[4]); 
}

What do you think?

Cheers, Benno

5

Hi Benno,

operator()(x, p) is const in order to distinguish from a call to SetParameters() and then operator()(x).
The const is important in order to work in a multi-thread environment, for example when making a parallel fits. For this, I also don’t like making the variable mutable.
Yes, the best is probably to implement as you suggest a VavilovPdf (and maybe VavilovCdf classes). I would probably there implement both a double operator() (x,p) const and non-const. In the first case I would re-create the Vavilov object instead of re-using as a mutable, so the method is really const. I am not sure, but probably the majority of time will be spent anyway in the re-evaluation of the working arrays when changing the parameters and not so much in their allocation.

Cheers

Lorenzo