EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-01-30 10:22:10
0001 // @(#)root/mathmore:$Id$ 0002 // Authors: B. List 29.4.2010 0003 0004 /********************************************************************** 0005 * * 0006 * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * 0007 * * 0008 * This library is free software; you can redistribute it and/or * 0009 * modify it under the terms of the GNU General Public License * 0010 * as published by the Free Software Foundation; either version 2 * 0011 * of the License, or (at your option) any later version. * 0012 * * 0013 * This library is distributed in the hope that it will be useful, * 0014 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 0015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * 0016 * General Public License for more details. * 0017 * * 0018 * You should have received a copy of the GNU General Public License * 0019 * along with this library (see file COPYING); if not, write * 0020 * to the Free Software Foundation, Inc., 59 Temple Place, Suite * 0021 * 330, Boston, MA 02111-1307 USA, or contact the author. * 0022 * * 0023 **********************************************************************/ 0024 0025 // Header file for class VavilovAccurate 0026 // 0027 // Created by: blist at Thu Apr 29 11:19:00 2010 0028 // 0029 // Last update: Thu Apr 29 11:19:00 2010 0030 // 0031 #ifndef ROOT_Math_VavilovAccurate 0032 #define ROOT_Math_VavilovAccurate 0033 0034 0035 #include "Math/Vavilov.h" 0036 0037 namespace ROOT { 0038 namespace Math { 0039 0040 //____________________________________________________________________________ 0041 /** 0042 Class describing a Vavilov distribution. 0043 0044 The probability density function of the Vavilov distribution 0045 as function of Landau's parameter is given by: 0046 \f[ p(\lambda_L; \kappa, \beta^2) = 0047 \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} \phi(s) e^{\lambda_L s} ds\f] 0048 where \f$\phi(s) = e^{C} e^{\psi(s)}\f$ 0049 with \f$ C = \kappa (1+\beta^2 \gamma )\f$ 0050 and \f$\psi(s)= s \ln \kappa + (s+\beta^2 \kappa) 0051 \cdot \left ( \int \limits_{0}^{1} 0052 \frac{1 - e^{\frac{-st}{\kappa}}}{t} \,d t- \gamma \right ) 0053 - \kappa \, e^{\frac{-s}{\kappa}}\f$. 0054 \f$ \gamma = 0.5772156649\dots\f$ is Euler's constant. 0055 0056 For the class VavilovAccurate, 0057 Pdf returns the Vavilov distribution as function of Landau's parameter 0058 \f$\lambda_L = \lambda_V/\kappa - \ln \kappa\f$, 0059 which is the convention used in the CERNLIB routines, and in the tables 0060 by S.M. Seltzer and M.J. Berger: Energy loss stragglin of protons and mesons: 0061 Tabulation of the Vavilov distribution, pp 187-203 0062 in: National Research Council (U.S.), Committee on Nuclear Science: 0063 Studies in penetration of charged particles in matter, 0064 Nat. Akad. Sci. Publication 1133, 0065 Nucl. Sci. Series Report No. 39, 0066 Washington (Nat. Akad. Sci.) 1964, 388 pp. 0067 Available from 0068 <A HREF="http://books.google.de/books?id=kmMrAAAAYAAJ&lpg=PP9&pg=PA187#v=onepage&q&f=false">Google books</A> 0069 0070 Therefore, for small values of \f$\kappa < 0.01\f$, 0071 pdf approaches the Landau distribution. 0072 0073 For values \f$\kappa > 10\f$, the Gauss approximation should be used 0074 with \f$\mu\f$ and \f$\sigma\f$ given by Vavilov::mean(kappa, beta2) 0075 and sqrt(Vavilov::variance(kappa, beta2). 0076 0077 The original Vavilov pdf is obtained by 0078 v.Pdf(lambdaV/kappa-log(kappa))/kappa. 0079 0080 For detailed description see 0081 B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers, 0082 <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>, 0083 which has been implemented in 0084 <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g116/top.html"> 0085 CERNLIB (G116)</A>. 0086 0087 The class stores coefficients needed to calculate \f$p(\lambda; \kappa, \beta^2)\f$ 0088 for fixed values of \f$\kappa\f$ and \f$\beta^2\f$. 0089 Changing these values is computationally expensive. 0090 0091 The parameter \f$\kappa\f$ should be in the range \f$0.01 \le \kappa \le 10\f$. 0092 In contrast to the CERNLIB implementation, all values of \f$\kappa \ge 0.001\f$ may be used, 0093 but may result in slower running and/or inaccurate results. 0094 0095 The parameter \f$\beta^2\f$ must be in the range \f$0 \le \beta^2 \le 1\f$. 0096 0097 Two parameters which are fixed in the CERNLIB implementation may be set by the user: 0098 - epsilonPM corresponds to \f$\epsilon^+ = \epsilon^-\f$ in Eqs. (2.1) and (2.2) of Schorr's paper. 0099 epsilonPM gives an estimate on the integral of the cumulative distribution function 0100 outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$ 0101 where the approximation is valid. 0102 Thus, it determines the support of the approximation used here (called $T_0 - T_1$ in the paper). 0103 Schorr recommends \f$\epsilon^+ = \epsilon^- = 5\cdot 10^{-4}\f$. 0104 The code from CERNLIB has been extended such that also smaller values are possible. 0105 0106 - epsilon corresponds to \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper. 0107 It determines the accuracy of the series expansion. 0108 Schorr recommends \f$\epsilon = 10^{-5}\f$. 0109 0110 For the quantile calculation, the algorithm given by Schorr is not used, 0111 because it turns out to be very slow and still inaccurate. 0112 Instead, an initial estimate is calculated based on a pre-calculated table, 0113 which is subsequently improved by Newton iterations. 0114 0115 While the CERNLIB implementation calculates at most 156 terms in the series expansion 0116 for the pdf and cdf calculation, this class calculates up to 500 terms, depending 0117 on the values of epsilonPM and epsilon. 0118 0119 Average times on a Pentium Core2 Duo P8400 2.26GHz: 0120 - 38us per call to SetKappaBeta2 or constructor 0121 - 0.49us per call to Pdf, Cdf 0122 - 8.2us per first call to Quantile after SetKappaBeta2 or constructor 0123 - 0.83us per subsequent call to Quantile 0124 0125 Benno List, June 2010 0126 0127 @ingroup StatFunc 0128 */ 0129 0130 0131 class VavilovAccurate: public Vavilov { 0132 0133 public: 0134 0135 0136 /** 0137 Initialize an object to calculate the Vavilov distribution 0138 0139 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0140 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0141 @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 0142 outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$ 0143 where the approximation is valid. 0144 @param epsilon: \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper; determines the accuracy of the series expansion. 0145 */ 0146 0147 VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5); 0148 0149 /** 0150 Destructor 0151 */ 0152 ~VavilovAccurate() override; 0153 0154 0155 public: 0156 0157 /** 0158 Evaluate the Vavilov probability density function 0159 0160 @param x The Landau parameter \f$x = \lambda_L\f$ 0161 */ 0162 double Pdf (double x) const override; 0163 0164 /** 0165 Evaluate the Vavilov probability density function, 0166 and set kappa and beta2, if necessary 0167 0168 @param x The Landau parameter \f$x = \lambda_L\f$ 0169 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0170 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0171 */ 0172 double Pdf (double x, double kappa, double beta2) override; 0173 0174 /** 0175 Evaluate the Vavilov cumulative probability density function 0176 0177 @param x The Landau parameter \f$x = \lambda_L\f$ 0178 */ 0179 double Cdf (double x) const override; 0180 0181 /** 0182 Evaluate the Vavilov cumulative probability density function, 0183 and set kappa and beta2, if necessary 0184 0185 @param x The Landau parameter \f$x = \lambda_L\f$ 0186 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0187 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0188 */ 0189 double Cdf (double x, double kappa, double beta2) override; 0190 0191 /** 0192 Evaluate the Vavilov complementary cumulative probability density function 0193 0194 @param x The Landau parameter \f$x = \lambda_L\f$ 0195 */ 0196 double Cdf_c (double x) const override; 0197 0198 /** 0199 Evaluate the Vavilov complementary cumulative probability density function, 0200 and set kappa and beta2, if necessary 0201 0202 @param x The Landau parameter \f$x = \lambda_L\f$ 0203 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0204 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0205 */ 0206 double Cdf_c (double x, double kappa, double beta2) override; 0207 0208 /** 0209 Evaluate the inverse of the Vavilov cumulative probability density function 0210 0211 @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$ 0212 */ 0213 double Quantile (double z) const override; 0214 0215 /** 0216 Evaluate the inverse of the Vavilov cumulative probability density function, 0217 and set kappa and beta2, if necessary 0218 0219 @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$ 0220 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0221 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0222 */ 0223 double Quantile (double z, double kappa, double beta2) override; 0224 0225 /** 0226 Evaluate the inverse of the complementary Vavilov cumulative probability density function 0227 0228 @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$ 0229 */ 0230 double Quantile_c (double z) const override; 0231 0232 /** 0233 Evaluate the inverse of the complementary Vavilov cumulative probability density function, 0234 and set kappa and beta2, if necessary 0235 0236 @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$ 0237 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0238 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0239 */ 0240 double Quantile_c (double z, double kappa, double beta2) override; 0241 0242 /** 0243 Change \f$\kappa\f$ and \f$\beta^2\f$ and recalculate coefficients if necessary 0244 0245 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0246 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0247 */ 0248 void SetKappaBeta2 (double kappa, double beta2) override; 0249 0250 0251 /** 0252 (Re)Initialize the object 0253 0254 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0255 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0256 @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 0257 outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$ 0258 where the approximation is valid. 0259 @param epsilon \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper; determines the accuracy of the series expansion. 0260 */ 0261 void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5); 0262 0263 0264 /** 0265 Return the minimum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$ 0266 is nonzero in the current approximation 0267 */ 0268 double GetLambdaMin() const override; 0269 0270 /** 0271 Return the maximum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$ 0272 is nonzero in the current approximation 0273 */ 0274 double GetLambdaMax() const override; 0275 0276 /** 0277 Return the current value of \f$\kappa\f$ 0278 */ 0279 double GetKappa() const override; 0280 0281 /** 0282 Return the current value of \f$\beta^2\f$ 0283 */ 0284 double GetBeta2() const override; 0285 0286 /** 0287 Return the value of \f$\lambda\f$ where the pdf is maximal 0288 */ 0289 double Mode() const override; 0290 0291 /** 0292 Return the value of \f$\lambda\f$ where the pdf is maximal function, 0293 and set kappa and beta2, if necessary 0294 0295 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0296 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0297 */ 0298 double Mode(double kappa, double beta2) override; 0299 0300 /** 0301 Return the current value of \f$\epsilon^+ = \epsilon^-\f$ 0302 */ 0303 0304 double GetEpsilonPM() const; 0305 0306 /** 0307 Return the current value of \f$\epsilon\f$ 0308 */ 0309 double GetEpsilon() const; 0310 0311 /** 0312 Return the number of terms used in the series expansion 0313 */ 0314 double GetNTerms() const; 0315 0316 /** 0317 Returns a static instance of class VavilovFast 0318 */ 0319 static VavilovAccurate *GetInstance(); 0320 0321 /** 0322 Returns a static instance of class VavilovFast, 0323 and sets the values of kappa and beta2 0324 0325 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0326 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0327 */ 0328 static VavilovAccurate *GetInstance(double kappa, double beta2); 0329 0330 0331 private: 0332 constexpr static int MAXTERMS{500}; 0333 double fH[8], fT0, fT1, fT, fOmega, fA_pdf[MAXTERMS+1], fB_pdf[MAXTERMS+1], fA_cdf[MAXTERMS+1], fB_cdf[MAXTERMS+1], fX0; 0334 double fKappa, fBeta2; 0335 double fEpsilonPM, fEpsilon; 0336 0337 mutable bool fQuantileInit; 0338 mutable int fNQuant; 0339 constexpr static int kNquantMax{32}; 0340 mutable double fQuant[kNquantMax]; 0341 mutable double fLambda[kNquantMax]; 0342 0343 void InitQuantile() const; 0344 0345 static VavilovAccurate *fgInstance; 0346 0347 double G116f1 (double x) const; 0348 double G116f2 (double x) const; 0349 0350 int Rzero (double a, double b, double& x0, 0351 double eps, int mxf, double (VavilovAccurate::*f)(double)const) const; 0352 static double E1plLog (double x); // Calculates log(|x|)+E_1(x) 0353 0354 }; 0355 0356 /** 0357 The Vavilov probability density function 0358 0359 @param x The Landau parameter \f$x = \lambda_L\f$ 0360 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0361 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0362 0363 @ingroup PdfFunc 0364 */ 0365 double vavilov_accurate_pdf (double x, double kappa, double beta2); 0366 0367 /** 0368 The Vavilov cumulative probability density function 0369 0370 @param x The Landau parameter \f$x = \lambda_L\f$ 0371 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0372 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0373 0374 @ingroup ProbFunc 0375 */ 0376 double vavilov_accurate_cdf (double x, double kappa, double beta2); 0377 0378 /** 0379 The Vavilov complementary cumulative probability density function 0380 0381 @param x The Landau parameter \f$x = \lambda_L\f$ 0382 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0383 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0384 0385 @ingroup ProbFunc 0386 */ 0387 double vavilov_accurate_cdf_c (double x, double kappa, double beta2); 0388 0389 /** 0390 The inverse of the Vavilov cumulative probability density function 0391 0392 @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$ 0393 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0394 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0395 0396 @ingroup QuantFunc 0397 */ 0398 double vavilov_accurate_quantile (double z, double kappa, double beta2); 0399 0400 /** 0401 The inverse of the complementary Vavilov cumulative probability density function 0402 0403 @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$ 0404 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$ 0405 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$ 0406 0407 @ingroup QuantFunc 0408 */ 0409 double vavilov_accurate_quantile_c (double z, double kappa, double beta2); 0410 0411 } // namespace Math 0412 } // namespace ROOT 0413 0414 #endif /* ROOT_Math_VavilovAccurate */
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|