ROOT   6.08/07 Reference Guide
VavilovAccurate.h
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: B. List 29.4.2010
3
4  /**********************************************************************
5  * *
6  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24
25 // Header file for class VavilovAccurate
26 //
27 // Created by: blist at Thu Apr 29 11:19:00 2010
28 //
29 // Last update: Thu Apr 29 11:19:00 2010
30 //
31 #ifndef ROOT_Math_VavilovAccurate
32 #define ROOT_Math_VavilovAccurate
33
34
35 #include "Math/Vavilov.h"
36
37 namespace ROOT {
38 namespace Math {
39
40 //____________________________________________________________________________
41 /**
42  Class describing a Vavilov distribution.
43
44  The probability density function of the Vavilov distribution
45  as function of Landau's parameter is given by:
46  \f[ p(\lambda_L; \kappa, \beta^2) =
47  \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} \phi(s) e^{\lambda_L s} ds\f]
48  where \f$\phi(s) = e^{C} e^{\psi(s)}\f$
49  with \f$C = \kappa (1+\beta^2 \gamma )\f$
50  and \f$\psi(s)= s \ln \kappa + (s+\beta^2 \kappa) 51 \cdot \left ( \int \limits_{0}^{1} 52 \frac{1 - e^{\frac{-st}{\kappa}}}{t} \,d t- \gamma \right ) 53 - \kappa \, e^{\frac{-s}{\kappa}}\f$.
54  \f$\gamma = 0.5772156649\dots\f$ is Euler's constant.
55
56  For the class VavilovAccurate,
57  Pdf returns the Vavilov distribution as function of Landau's parameter
58  \f$\lambda_L = \lambda_V/\kappa - \ln \kappa\f$,
59  which is the convention used in the CERNLIB routines, and in the tables
60  by S.M. Seltzer and M.J. Berger: Energy loss stragglin of protons and mesons:
61  Tabulation of the Vavilov distribution, pp 187-203
62  in: National Research Council (U.S.), Committee on Nuclear Science:
63  Studies in penetration of charged particles in matter,
64  Nat. Akad. Sci. Publication 1133,
65  Nucl. Sci. Series Report No. 39,
66  Washington (Nat. Akad. Sci.) 1964, 388 pp.
67  Available from
69
70  Therefore, for small values of \f$\kappa < 0.01\f$,
71  pdf approaches the Landau distribution.
72
73  For values \f$\kappa > 10\f$, the Gauss approximation should be used
74  with \f$\mu\f$ and \f$\sigma\f$ given by Vavilov::mean(kappa, beta2)
75  and sqrt(Vavilov::variance(kappa, beta2).
76
77  The original Vavilov pdf is obtained by
78  v.Pdf(lambdaV/kappa-log(kappa))/kappa.
79
80  For detailed description see
81  B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
82  <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>,
83  which has been implemented in
84  <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g116/top.html">
85  CERNLIB (G116)</A>.
86
87  The class stores coefficients needed to calculate \f$p(\lambda; \kappa, \beta^2)\f$
88  for fixed values of \f$\kappa\f$ and \f$\beta^2\f$.
89  Changing these values is computationally expensive.
90
91  The parameter \f$\kappa\f$ should be in the range \f$0.01 \le \kappa \le 10\f$.
92  In contrast to the CERNLIB implementation, all values of \f$\kappa \ge 0.001\f$ may be used,
93  but may result in slower running and/or inaccurate results.
94
95  The parameter \f$\beta^2\f$ must be in the range \f$0 \le \beta^2 \le 1\f$.
96
97  Two parameters which are fixed in the CERNLIB implementation may be set by the user:
98  - epsilonPM corresponds to \f$\epsilon^+ = \epsilon^-\f$ in Eqs. (2.1) and (2.2) of Schorr's paper.
99  epsilonPM gives an estimate on the integral of the cumulative distribution function
100  outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
101  where the approximation is valid.
102  Thus, it determines the support of the approximation used here (called $T_0 - T_1$ in the paper).
103  Schorr recommends \f$\epsilon^+ = \epsilon^- = 5\cdot 10^{-4}\f$.
104  The code from CERNLIB has been extended such that also smaller values are possible.
105
106  - epsilon corresponds to \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper.
107  It determines the accuracy of the series expansion.
108  Schorr recommends \f$\epsilon = 10^{-5}\f$.
109
110  For the quantile calculation, the algorithm given by Schorr is not used,
111  because it turns out to be very slow and still inaccurate.
112  Instead, an initial estimate is calculated based on a pre-calculated table,
113  which is subsequently improved by Newton iterations.
114
115  While the CERNLIB implementation calculates at most 156 terms in the series expansion
116  for the pdf and cdf calculation, this class calculates up to 500 terms, depending
117  on the values of epsilonPM and epsilon.
118
119  Average times on a Pentium Core2 Duo P8400 2.26GHz:
120  - 38us per call to SetKappaBeta2 or constructor
121  - 0.49us per call to Pdf, Cdf
122  - 8.2us per first call to Quantile after SetKappaBeta2 or constructor
123  - 0.83us per subsequent call to Quantile
124
125  Benno List, June 2010
126
127  @ingroup StatFunc
128  */
129
130
131 class VavilovAccurate: public Vavilov {
132
133 public:
134
135
136  /**
137  Initialize an object to calculate the Vavilov distribution
138
139  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
140  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
141  @param epsilonPM: \f$\epsilon^+ = \epsilon^-\f$ in Eqs. (2.1) and (2.2) of Schorr's paper; gives an estimate on the integral of the cumulative distribution function
142  outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
143  where the approximation is valid.
144  @param epsilon: \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper; determines the accuracy of the series expansion.
145  */
146
147  VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5);
148
149  /**
150  Destructor
151  */
152  virtual ~VavilovAccurate();
153
154
155 public:
156
157  /**
158  Evaluate the Vavilov probability density function
159
160  @param x The Landau parameter \f$x = \lambda_L\f$
161  */
162  double Pdf (double x) const;
163
164  /**
165  Evaluate the Vavilov probability density function,
166  and set kappa and beta2, if necessary
167
168  @param x The Landau parameter \f$x = \lambda_L\f$
169  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
170  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
171  */
172  double Pdf (double x, double kappa, double beta2);
173
174  /**
175  Evaluate the Vavilov cumulative probability density function
176
177  @param x The Landau parameter \f$x = \lambda_L\f$
178  */
179  double Cdf (double x) const;
180
181  /**
182  Evaluate the Vavilov cumulative probability density function,
183  and set kappa and beta2, if necessary
184
185  @param x The Landau parameter \f$x = \lambda_L\f$
186  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
187  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
188  */
189  double Cdf (double x, double kappa, double beta2);
190
191  /**
192  Evaluate the Vavilov complementary cumulative probability density function
193
194  @param x The Landau parameter \f$x = \lambda_L\f$
195  */
196  double Cdf_c (double x) const;
197
198  /**
199  Evaluate the Vavilov complementary cumulative probability density function,
200  and set kappa and beta2, if necessary
201
202  @param x The Landau parameter \f$x = \lambda_L\f$
203  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
204  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
205  */
206  double Cdf_c (double x, double kappa, double beta2);
207
208  /**
209  Evaluate the inverse of the Vavilov cumulative probability density function
210
211  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
212  */
213  double Quantile (double z) const;
214
215  /**
216  Evaluate the inverse of the Vavilov cumulative probability density function,
217  and set kappa and beta2, if necessary
218
219  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
220  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
221  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
222  */
223  double Quantile (double z, double kappa, double beta2);
224
225  /**
226  Evaluate the inverse of the complementary Vavilov cumulative probability density function
227
228  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
229  */
230  double Quantile_c (double z) const;
231
232  /**
233  Evaluate the inverse of the complementary Vavilov cumulative probability density function,
234  and set kappa and beta2, if necessary
235
236  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
237  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
238  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
239  */
240  double Quantile_c (double z, double kappa, double beta2);
241
242  /**
243  Change \f$\kappa\f$ and \f$\beta^2\f$ and recalculate coefficients if necessary
244
245  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
246  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
247  */
248  virtual void SetKappaBeta2 (double kappa, double beta2);
249
250
251  /**
252  (Re)Initialize the object
253
254  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
255  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
256  @param epsilonPM \f$\epsilon^+ = \epsilon^-\f$ in Eqs. (2.1) and (2.2) of Schorr's paper; gives an estimate on the integral of the cumulative distribution function
257  outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
258  where the approximation is valid.
259  @param epsilon \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper; determines the accuracy of the series expansion.
260  */
261  void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5);
262
263
264  /**
265  Return the minimum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
266  is nonzero in the current approximation
267  */
268  virtual double GetLambdaMin() const;
269
270  /**
271  Return the maximum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
272  is nonzero in the current approximation
273  */
274  virtual double GetLambdaMax() const;
275
276  /**
277  Return the current value of \f$\kappa\f$
278  */
279  virtual double GetKappa() const;
280
281  /**
282  Return the current value of \f$\beta^2\f$
283  */
284  virtual double GetBeta2() const;
285
286  /**
287  Return the value of \f$\lambda\f$ where the pdf is maximal
288  */
289  virtual double Mode() const;
290
291  /**
292  Return the value of \f$\lambda\f$ where the pdf is maximal function,
293  and set kappa and beta2, if necessary
294
295  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
296  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
297  */
298  virtual double Mode(double kappa, double beta2);
299
300  /**
301  Return the current value of \f$\epsilon^+ = \epsilon^-\f$
302  */
303
304  double GetEpsilonPM() const;
305
306  /**
307  Return the current value of \f$\epsilon\f$
308  */
309  double GetEpsilon() const;
310
311  /**
312  Return the number of terms used in the series expansion
313  */
314  double GetNTerms() const;
315
316  /**
317  Returns a static instance of class VavilovFast
318  */
319  static VavilovAccurate *GetInstance();
320
321  /**
322  Returns a static instance of class VavilovFast,
323  and sets the values of kappa and beta2
324
325  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
326  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
327  */
328  static VavilovAccurate *GetInstance(double kappa, double beta2);
329
330
331 private:
332  enum{MAXTERMS=500};
334  double fKappa, fBeta2;
336
337  mutable bool fQuantileInit;
338  mutable int fNQuant;
339  enum{kNquantMax=32};
340  mutable double fQuant[kNquantMax];
341  mutable double fLambda[kNquantMax];
342
343  void InitQuantile() const;
344
346
347  double G116f1 (double x) const;
348  double G116f2 (double x) const;
349
350  int Rzero (double a, double b, double& x0,
351  double eps, int mxf, double (VavilovAccurate::*f)(double)const) const;
352  static double E1plLog (double x); // Calculates log(|x|)+E_1(x)
353
354 };
355
356  /**
357  The Vavilov probability density function
358
359  @param x The Landau parameter \f$x = \lambda_L\f$
360  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
361  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
362
363  @ingroup PdfFunc
364  */
365 double vavilov_accurate_pdf (double x, double kappa, double beta2);
366
367  /**
368  The Vavilov cumulative probability density function
369
370  @param x The Landau parameter \f$x = \lambda_L\f$
371  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
372  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
373
374  @ingroup ProbFunc
375  */
376 double vavilov_accurate_cdf (double x, double kappa, double beta2);
377
378  /**
379  The Vavilov complementary cumulative probability density function
380
381  @param x The Landau parameter \f$x = \lambda_L\f$
382  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
383  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
384
385  @ingroup ProbFunc
386  */
387 double vavilov_accurate_cdf_c (double x, double kappa, double beta2);
388
389  /**
390  The inverse of the Vavilov cumulative probability density function
391
392  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
393  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
394  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
395
396  @ingroup QuantFunc
397  */
398 double vavilov_accurate_quantile (double z, double kappa, double beta2);
399
400  /**
401  The inverse of the complementary Vavilov cumulative probability density function
402
403  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
404  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
405  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
406
407  @ingroup QuantFunc
408  */
409 double vavilov_accurate_quantile_c (double z, double kappa, double beta2);
410
411 } // namespace Math
412 } // namespace ROOT
413
414 #endif /* ROOT_Math_VavilovAccurate */
Base class describing a Vavilov distribution.
Definition: Vavilov.h:121
virtual double GetKappa() const
Return the current value of .
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
TArc * a
Definition: textangle.C:12
int Rzero(double a, double b, double &x0, double eps, int mxf, double(VavilovAccurate::*f)(double) const) const
double G116f2(double x) const
double Cdf(double x) const
Evaluate the Vavilov cumulative probability density function.
double GetNTerms() const
Return the number of terms used in the series expansion.
double vavilov_accurate_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
double GetEpsilonPM() const
Return the current value of .
Double_t x[n]
Definition: legend1.C:17
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
double GetEpsilon() const
Return the current value of .
void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5)
(Re)Initialize the object
static VavilovAccurate * GetInstance()
Returns a static instance of class VavilovFast.
double Quantile(double z) const
Evaluate the inverse of the Vavilov cumulative probability density function.
double Cdf_c(double x) const
Evaluate the Vavilov complementary cumulative probability density function.
virtual double Mode() const
Return the value of where the pdf is maximal.
static VavilovAccurate * fgInstance
Double_t E()
Definition: TMath.h:54
VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5)
Initialize an object to calculate the Vavilov distribution.
REAL epsilon
Definition: triangle.c:617
double f(double x)
double vavilov_accurate_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
Class describing a Vavilov distribution.
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
double G116f1(double x) const
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
virtual double GetBeta2() const
Return the current value of .
static double E1plLog(double x)
double Pdf(double x) const
Evaluate the Vavilov probability density function.
double fLambda[kNquantMax]
double Quantile_c(double z) const
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
virtual ~VavilovAccurate()
Destructor.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
double vavilov_accurate_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function. ...
double vavilov_accurate_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cumulative probability density function.
double vavilov_accurate_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.