Logo ROOT   6.16/01
Reference Guide
QuantFuncMathCore.cxx
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Authors: L. Moneta, A. Zsenei 08/2005
3
4
5#include "Math/Math.h"
7#include "SpecFuncCephes.h"
8#include <limits>
9
10
11namespace ROOT {
12namespace Math {
13
14
15
16 double beta_quantile_c(double z, double a, double b) {
17 // use Cephes and proprety of icomplete beta function
18 if ( z < 0.5)
19 return 1. - ROOT::Math::Cephes::incbi(b,a,z);
20 else
21 return ROOT::Math::Cephes::incbi(a,b,1.0-z);
22
23 }
24
25
26 double beta_quantile(double z, double a, double b ) {
27 // use Cephes function
29
30 }
31
32
33 double cauchy_quantile_c(double z, double b) {
34 // inverse of Caucy is simply the tan(PI(z-0.5))
35 if (z == 0) return std::numeric_limits<double>::infinity();
36 if (z == 1) return - std::numeric_limits<double>::infinity();
37 if (z < 0.5)
38 // use fact that tan(PI(0.5-z)) = 1/tan(PI*z)
39 return b / std::tan( M_PI * z );
40 else
41 return b * std::tan( M_PI * (0.5 - z ) );
42 }
43
44
45
46 double cauchy_quantile(double z, double b) {
47 // inverse of Caucy is simply the tan(PI(z-0.5))
48 if (z == 0) return - std::numeric_limits<double>::infinity();
49 if (z == 1) return + std::numeric_limits<double>::infinity();
50 if (z < 0.5)
51 // use fact that tan(PI(0.5-z)) = 1/tan(PI*z)
52 return - b / std::tan( M_PI * z );
53 else
54 return b * std::tan( M_PI * ( z - 0.5 ) );
55
56 }
57
58
59
60 double chisquared_quantile_c(double z, double r) {
61 // use Cephes igami which return inverse of complemented regularized gamma
62 return 2.* ROOT::Math::Cephes::igami( 0.5 *r, z);
63
64 }
65
66
67 double chisquared_quantile(double z, double r) {
68 // if possible, should use MathMore function ROOT::Math::chisquared_quantile for z close to zero
69 // otherwise will always return zero for z value smaller than eps
70 return 2.* ROOT::Math::Cephes::igami( 0.5 *r, 1. - z);
71 }
72
73
74 double exponential_quantile_c(double z, double lambda) {
75
76 return - std::log(z)/ lambda;
77
78 }
79
80
81
82 double exponential_quantile(double z, double lambda) {
83 // use log1p for avoid errors at small z
84 return - ROOT::Math::log1p(-z)/lambda;
85
86 }
87
88
89 double fdistribution_quantile_c(double z, double n, double m) {
90 // use cephes incbi function and use propreties of incomplete beta for case <> 0.5
91 if (n == 0) return 0; // is value of cdf for n = 0
92 if (z < 0.5) {
93 double y = ROOT::Math::Cephes::incbi( .5*m, .5*n, z);
94 return m/(n * y) - m/n;
95 }
96 else {
97 double y = ROOT::Math::Cephes::incbi( .5*n, .5*m, 1.0 - z);
98 // will lose precision for y approx to 1
99 return m * y /(n * ( 1. - y) );
100 }
101 }
102
103 double fdistribution_quantile(double z, double n, double m) {
104 // use cephes incbi function
105 if (n == 0) return 0; // is value of cdf for n = 0
106 double y = ROOT::Math::Cephes::incbi( .5*n, .5*m, z);
107 // will lose precision for y approx to 1
108 return m * y /(n * ( 1. - y) );
109 }
110
111
112 double gamma_quantile_c(double z, double alpha, double theta) {
113
114 return theta * ROOT::Math::Cephes::igami( alpha, z);
115
116 }
117
118 double gamma_quantile(double z, double alpha, double theta) {
119 // if possible, should use MathMore function ROOT::Math::gamma_quantile for z close to zero
120 // otherwise will always return zero for z value smaller than eps
121 return theta * ROOT::Math::Cephes::igami( alpha, 1.- z);
122 }
123
124
125
126 double normal_quantile_c(double z, double sigma) {
127 // use cephes and fact that ntri(1.-x) = - ndtri(x)
128 return - sigma * ROOT::Math::Cephes::ndtri(z);
129
130 }
131
132
133
134 double normal_quantile(double z, double sigma) {
135 // use cephes ndtri function
137
138 }
139
140
141
142
143 double lognormal_quantile_c(double z, double m, double s) {
144 // if y is log normal, u = exp(y) is log-normal distributed
145 double y = - s * ROOT::Math::Cephes::ndtri(z) + m;
146 return std::exp(y);
147 }
148
149
150
151 double lognormal_quantile(double z, double m, double s) {
152 // if y is log normal, u = exp(y) is log-normal distributed
153 double y = s * ROOT::Math::Cephes::ndtri(z) + m;
154 return std::exp(y);
155
156 }
157
158
159// double tdistribution_quantile_c(double z, double r) {
160
161// return gsl_cdf_tdist_Qinv(z, r);
162
163// }
164
165
166
167// double tdistribution_quantile(double z, double r) {
168
169// return gsl_cdf_tdist_Pinv(z, r);
170
171// }
172
173
174
175 double uniform_quantile_c(double z, double a, double b) {
176
177 return a * z + b * (1.0 - z);
178
179 }
180
181
182
183 double uniform_quantile(double z, double a, double b) {
184
185 return b * z + a * (1.0 - z);
186
187 }
188
189 double landau_quantile(double z, double xi) {
190 // LANDAU quantile : algorithm from CERNLIB G110 ranlan
191 // with scale parameter xi
192 // Converted by Rene Brun from CERNLIB routine ranlan(G110),
193 // Moved and adapted to QuantFuncMathCore by B. List 29.4.2010
194
195 static double f[982] = {
196 0 , 0 , 0 ,0 ,0 ,-2.244733,
197 -2.204365,-2.168163,-2.135219,-2.104898,-2.076740,-2.050397,
198 -2.025605,-2.002150,-1.979866,-1.958612,-1.938275,-1.918760,
199 -1.899984,-1.881879,-1.864385,-1.847451,-1.831030,-1.815083,
200 -1.799574,-1.784473,-1.769751,-1.755383,-1.741346,-1.727620,
201 -1.714187,-1.701029,-1.688130,-1.675477,-1.663057,-1.650858,
202 -1.638868,-1.627078,-1.615477,-1.604058,-1.592811,-1.581729,
203 -1.570806,-1.560034,-1.549407,-1.538919,-1.528565,-1.518339,
204 -1.508237,-1.498254,-1.488386,-1.478628,-1.468976,-1.459428,
205 -1.449979,-1.440626,-1.431365,-1.422195,-1.413111,-1.404112,
206 -1.395194,-1.386356,-1.377594,-1.368906,-1.360291,-1.351746,
207 -1.343269,-1.334859,-1.326512,-1.318229,-1.310006,-1.301843,
208 -1.293737,-1.285688,-1.277693,-1.269752,-1.261863,-1.254024,
209 -1.246235,-1.238494,-1.230800,-1.223153,-1.215550,-1.207990,
210 -1.200474,-1.192999,-1.185566,-1.178172,-1.170817,-1.163500,
211 -1.156220,-1.148977,-1.141770,-1.134598,-1.127459,-1.120354,
212 -1.113282,-1.106242,-1.099233,-1.092255,
213 -1.085306,-1.078388,-1.071498,-1.064636,-1.057802,-1.050996,
214 -1.044215,-1.037461,-1.030733,-1.024029,-1.017350,-1.010695,
215 -1.004064, -.997456, -.990871, -.984308, -.977767, -.971247,
216 -.964749, -.958271, -.951813, -.945375, -.938957, -.932558,
217 -.926178, -.919816, -.913472, -.907146, -.900838, -.894547,
218 -.888272, -.882014, -.875773, -.869547, -.863337, -.857142,
219 -.850963, -.844798, -.838648, -.832512, -.826390, -.820282,
220 -.814187, -.808106, -.802038, -.795982, -.789940, -.783909,
221 -.777891, -.771884, -.765889, -.759906, -.753934, -.747973,
222 -.742023, -.736084, -.730155, -.724237, -.718328, -.712429,
223 -.706541, -.700661, -.694791, -.688931, -.683079, -.677236,
224 -.671402, -.665576, -.659759, -.653950, -.648149, -.642356,
225 -.636570, -.630793, -.625022, -.619259, -.613503, -.607754,
226 -.602012, -.596276, -.590548, -.584825, -.579109, -.573399,
227 -.567695, -.561997, -.556305, -.550618, -.544937, -.539262,
228 -.533592, -.527926, -.522266, -.516611, -.510961, -.505315,
229 -.499674, -.494037, -.488405, -.482777,
230 -.477153, -.471533, -.465917, -.460305, -.454697, -.449092,
231 -.443491, -.437893, -.432299, -.426707, -.421119, -.415534,
232 -.409951, -.404372, -.398795, -.393221, -.387649, -.382080,
233 -.376513, -.370949, -.365387, -.359826, -.354268, -.348712,
234 -.343157, -.337604, -.332053, -.326503, -.320955, -.315408,
235 -.309863, -.304318, -.298775, -.293233, -.287692, -.282152,
236 -.276613, -.271074, -.265536, -.259999, -.254462, -.248926,
237 -.243389, -.237854, -.232318, -.226783, -.221247, -.215712,
238 -.210176, -.204641, -.199105, -.193568, -.188032, -.182495,
239 -.176957, -.171419, -.165880, -.160341, -.154800, -.149259,
240 -.143717, -.138173, -.132629, -.127083, -.121537, -.115989,
241 -.110439, -.104889, -.099336, -.093782, -.088227, -.082670,
242 -.077111, -.071550, -.065987, -.060423, -.054856, -.049288,
243 -.043717, -.038144, -.032569, -.026991, -.021411, -.015828,
244 -.010243, -.004656, .000934, .006527, .012123, .017722,
245 .023323, .028928, .034535, .040146, .045759, .051376,
246 .056997, .062620, .068247, .073877,
247 .079511, .085149, .090790, .096435, .102083, .107736,
248 .113392, .119052, .124716, .130385, .136057, .141734,
249 .147414, .153100, .158789, .164483, .170181, .175884,
250 .181592, .187304, .193021, .198743, .204469, .210201,
251 .215937, .221678, .227425, .233177, .238933, .244696,
252 .250463, .256236, .262014, .267798, .273587, .279382,
253 .285183, .290989, .296801, .302619, .308443, .314273,
254 .320109, .325951, .331799, .337654, .343515, .349382,
255 .355255, .361135, .367022, .372915, .378815, .384721,
256 .390634, .396554, .402481, .408415, .414356, .420304,
257 .426260, .432222, .438192, .444169, .450153, .456145,
258 .462144, .468151, .474166, .480188, .486218, .492256,
259 .498302, .504356, .510418, .516488, .522566, .528653,
260 .534747, .540850, .546962, .553082, .559210, .565347,
261 .571493, .577648, .583811, .589983, .596164, .602355,
262 .608554, .614762, .620980, .627207, .633444, .639689,
263 .645945, .652210, .658484, .664768,
264 .671062, .677366, .683680, .690004, .696338, .702682,
265 .709036, .715400, .721775, .728160, .734556, .740963,
266 .747379, .753807, .760246, .766695, .773155, .779627,
267 .786109, .792603, .799107, .805624, .812151, .818690,
268 .825241, .831803, .838377, .844962, .851560, .858170,
269 .864791, .871425, .878071, .884729, .891399, .898082,
270 .904778, .911486, .918206, .924940, .931686, .938446,
271 .945218, .952003, .958802, .965614, .972439, .979278,
272 .986130, .992996, .999875, 1.006769, 1.013676, 1.020597,
273 1.027533, 1.034482, 1.041446, 1.048424, 1.055417, 1.062424,
274 1.069446, 1.076482, 1.083534, 1.090600, 1.097681, 1.104778,
275 1.111889, 1.119016, 1.126159, 1.133316, 1.140490, 1.147679,
276 1.154884, 1.162105, 1.169342, 1.176595, 1.183864, 1.191149,
277 1.198451, 1.205770, 1.213105, 1.220457, 1.227826, 1.235211,
278 1.242614, 1.250034, 1.257471, 1.264926, 1.272398, 1.279888,
279 1.287395, 1.294921, 1.302464, 1.310026, 1.317605, 1.325203,
280 1.332819, 1.340454, 1.348108, 1.355780,
281 1.363472, 1.371182, 1.378912, 1.386660, 1.394429, 1.402216,
282 1.410024, 1.417851, 1.425698, 1.433565, 1.441453, 1.449360,
283 1.457288, 1.465237, 1.473206, 1.481196, 1.489208, 1.497240,
284 1.505293, 1.513368, 1.521465, 1.529583, 1.537723, 1.545885,
285 1.554068, 1.562275, 1.570503, 1.578754, 1.587028, 1.595325,
286 1.603644, 1.611987, 1.620353, 1.628743, 1.637156, 1.645593,
287 1.654053, 1.662538, 1.671047, 1.679581, 1.688139, 1.696721,
288 1.705329, 1.713961, 1.722619, 1.731303, 1.740011, 1.748746,
289 1.757506, 1.766293, 1.775106, 1.783945, 1.792810, 1.801703,
290 1.810623, 1.819569, 1.828543, 1.837545, 1.846574, 1.855631,
291 1.864717, 1.873830, 1.882972, 1.892143, 1.901343, 1.910572,
292 1.919830, 1.929117, 1.938434, 1.947781, 1.957158, 1.966566,
293 1.976004, 1.985473, 1.994972, 2.004503, 2.014065, 2.023659,
294 2.033285, 2.042943, 2.052633, 2.062355, 2.072110, 2.081899,
295 2.091720, 2.101575, 2.111464, 2.121386, 2.131343, 2.141334,
296 2.151360, 2.161421, 2.171517, 2.181648, 2.191815, 2.202018,
297 2.212257, 2.222533, 2.232845, 2.243195,
298 2.253582, 2.264006, 2.274468, 2.284968, 2.295507, 2.306084,
299 2.316701, 2.327356, 2.338051, 2.348786, 2.359562, 2.370377,
300 2.381234, 2.392131, 2.403070, 2.414051, 2.425073, 2.436138,
301 2.447246, 2.458397, 2.469591, 2.480828, 2.492110, 2.503436,
302 2.514807, 2.526222, 2.537684, 2.549190, 2.560743, 2.572343,
303 2.583989, 2.595682, 2.607423, 2.619212, 2.631050, 2.642936,
304 2.654871, 2.666855, 2.678890, 2.690975, 2.703110, 2.715297,
305 2.727535, 2.739825, 2.752168, 2.764563, 2.777012, 2.789514,
306 2.802070, 2.814681, 2.827347, 2.840069, 2.852846, 2.865680,
307 2.878570, 2.891518, 2.904524, 2.917588, 2.930712, 2.943894,
308 2.957136, 2.970439, 2.983802, 2.997227, 3.010714, 3.024263,
309 3.037875, 3.051551, 3.065290, 3.079095, 3.092965, 3.106900,
310 3.120902, 3.134971, 3.149107, 3.163312, 3.177585, 3.191928,
311 3.206340, 3.220824, 3.235378, 3.250005, 3.264704, 3.279477,
312 3.294323, 3.309244, 3.324240, 3.339312, 3.354461, 3.369687,
313 3.384992, 3.400375, 3.415838, 3.431381, 3.447005, 3.462711,
314 3.478500, 3.494372, 3.510328, 3.526370,
315 3.542497, 3.558711, 3.575012, 3.591402, 3.607881, 3.624450,
316 3.641111, 3.657863, 3.674708, 3.691646, 3.708680, 3.725809,
317 3.743034, 3.760357, 3.777779, 3.795300, 3.812921, 3.830645,
318 3.848470, 3.866400, 3.884434, 3.902574, 3.920821, 3.939176,
319 3.957640, 3.976215, 3.994901, 4.013699, 4.032612, 4.051639,
320 4.070783, 4.090045, 4.109425, 4.128925, 4.148547, 4.168292,
321 4.188160, 4.208154, 4.228275, 4.248524, 4.268903, 4.289413,
322 4.310056, 4.330832, 4.351745, 4.372794, 4.393982, 4.415310,
323 4.436781, 4.458395, 4.480154, 4.502060, 4.524114, 4.546319,
324 4.568676, 4.591187, 4.613854, 4.636678, 4.659662, 4.682807,
325 4.706116, 4.729590, 4.753231, 4.777041, 4.801024, 4.825179,
326 4.849511, 4.874020, 4.898710, 4.923582, 4.948639, 4.973883,
327 4.999316, 5.024942, 5.050761, 5.076778, 5.102993, 5.129411,
328 5.156034, 5.182864, 5.209903, 5.237156, 5.264625, 5.292312,
329 5.320220, 5.348354, 5.376714, 5.405306, 5.434131, 5.463193,
330 5.492496, 5.522042, 5.551836, 5.581880, 5.612178, 5.642734,
331 5.673552, 5.704634, 5.735986, 5.767610,
332 5.799512, 5.831694, 5.864161, 5.896918, 5.929968, 5.963316,
333 5.996967, 6.030925, 6.065194, 6.099780, 6.134687, 6.169921,
334 6.205486, 6.241387, 6.277630, 6.314220, 6.351163, 6.388465,
335 6.426130, 6.464166, 6.502578, 6.541371, 6.580553, 6.620130,
336 6.660109, 6.700495, 6.741297, 6.782520, 6.824173, 6.866262,
337 6.908795, 6.951780, 6.995225, 7.039137, 7.083525, 7.128398,
338 7.173764, 7.219632, 7.266011, 7.312910, 7.360339, 7.408308,
339 7.456827, 7.505905, 7.555554, 7.605785, 7.656608, 7.708035,
340 7.760077, 7.812747, 7.866057, 7.920019, 7.974647, 8.029953,
341 8.085952, 8.142657, 8.200083, 8.258245, 8.317158, 8.376837,
342 8.437300, 8.498562, 8.560641, 8.623554, 8.687319, 8.751955,
343 8.817481, 8.883916, 8.951282, 9.019600, 9.088889, 9.159174,
344 9.230477, 9.302822, 9.376233, 9.450735, 9.526355, 9.603118,
345 9.681054, 9.760191, 9.840558, 9.922186,10.005107,10.089353,
346 10.174959,10.261958,10.350389,10.440287,10.531693,10.624646,
347 10.719188,10.815362,10.913214,11.012789,11.114137,11.217307,
348 11.322352,11.429325,11.538283,11.649285,
349 11.762390,11.877664,11.995170,12.114979,12.237161,12.361791,
350 12.488946,12.618708,12.751161,12.886394,13.024498,13.165570,
351 13.309711,13.457026,13.607625,13.761625,13.919145,14.080314,
352 14.245263,14.414134,14.587072,14.764233,14.945778,15.131877,
353 15.322712,15.518470,15.719353,15.925570,16.137345,16.354912,
354 16.578520,16.808433,17.044929,17.288305,17.538873,17.796967,
355 18.062943,18.337176,18.620068,18.912049,19.213574,19.525133,
356 19.847249,20.180480,20.525429,20.882738,21.253102,21.637266,
357 22.036036,22.450278,22.880933,23.329017,23.795634,24.281981,
358 24.789364,25.319207,25.873062,26.452634,27.059789,27.696581,
359 28.365274,29.068370,29.808638,30.589157,31.413354,32.285060,
360 33.208568,34.188705,35.230920,36.341388,37.527131,38.796172,
361 40.157721,41.622399,43.202525,44.912465,46.769077,48.792279,
362 51.005773,53.437996,56.123356,59.103894 };
363
364 if (xi <= 0) return 0;
365 if (z <= 0) return -std::numeric_limits<double>::infinity();
366 if (z >= 1) return std::numeric_limits<double>::infinity();
367
368 double ranlan, u, v;
369 u = 1000*z;
370 int i = int(u);
371 u -= i;
372 if (i >= 70 && i < 800) {
373 ranlan = f[i-1] + u*(f[i] - f[i-1]);
374 } else if (i >= 7 && i <= 980) {
375 ranlan = f[i-1] + u*(f[i]-f[i-1]-0.25*(1-u)*(f[i+1]-f[i]-f[i-1]+f[i-2]));
376 } else if (i < 7) {
377 v = std::log(z);
378 u = 1/v;
379 ranlan = ((0.99858950+(3.45213058E1+1.70854528E1*u)*u)/
380 (1 +(3.41760202E1+4.01244582 *u)*u))*
381 (-std::log(-0.91893853-v)-1);
382 } else {
383 u = 1-z;
384 v = u*u;
385 if (z <= 0.999) {
386 ranlan = (1.00060006+2.63991156E2*u+4.37320068E3*v)/
387 ((1 +2.57368075E2*u+3.41448018E3*v)*u);
388 } else {
389 ranlan = (1.00001538+6.07514119E3*u+7.34266409E5*v)/
390 ((1 +6.06511919E3*u+6.94021044E5*v)*u);
391 }
392 }
393 return xi*ranlan;
394 }
395
396 double landau_quantile_c(double z, double xi) {
397 return landau_quantile(1.-z,xi);
398 }
399
400} // namespace Math
401} // namespace ROOT
SVector< double, 2 > v
Definition: Dict.h:5
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define M_PI
Definition: Rotated.cxx:105
double tan(double)
double exp(double)
double log(double)
double fdistribution_quantile_c(double z, double n, double m)
Inverse ( ) of the cumulative distribution function of the upper tail of the f distribution (fdistrib...
double uniform_quantile_c(double z, double a, double b)
Inverse ( ) of the cumulative distribution function of the upper tail of the uniform (flat) distribut...
double landau_quantile(double z, double xi=1)
Inverse ( ) of the cumulative distribution function of the lower tail of the Landau distribution (lan...
double chisquared_quantile_c(double z, double r)
Inverse ( ) of the cumulative distribution function of the upper tail of the distribution with degr...
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
double exponential_quantile_c(double z, double lambda)
Inverse ( ) of the cumulative distribution function of the upper tail of the exponential distribution...
double gamma_quantile_c(double z, double alpha, double theta)
Inverse ( ) of the cumulative distribution function of the upper tail of the gamma distribution (gamm...
double landau_quantile_c(double z, double xi=1)
Inverse ( ) of the cumulative distribution function of the upper tail of the landau distribution (lan...
double fdistribution_quantile(double z, double n, double m)
Inverse ( ) of the cumulative distribution function of the lower tail of the f distribution (fdistrib...
double cauchy_quantile_c(double z, double b)
Inverse ( ) of the cumulative distribution function of the upper tail of the Cauchy distribution (cau...
double uniform_quantile(double z, double a, double b)
Inverse ( ) of the cumulative distribution function of the lower tail of the uniform (flat) distribut...
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
double chisquared_quantile(double z, double r)
Inverse ( ) of the cumulative distribution function of the lower tail of the distribution with degr...
double beta_quantile_c(double x, double a, double b)
Inverse ( ) of the cumulative distribution function of the lower tail of the beta distribution (beta_...
double gamma_quantile(double z, double alpha, double theta)
Inverse ( ) of the cumulative distribution function of the lower tail of the gamma distribution (gamm...
double exponential_quantile(double z, double lambda)
Inverse ( ) of the cumulative distribution function of the lower tail of the exponential distribution...
double lognormal_quantile(double x, double m, double s)
Inverse ( ) of the cumulative distribution function of the lower tail of the lognormal distribution (...
double cauchy_quantile(double z, double b)
Inverse ( ) of the cumulative distribution function of the lower tail of the Cauchy distribution (cau...
double beta_quantile(double x, double a, double b)
Inverse ( ) of the cumulative distribution function of the upper tail of the beta distribution (beta_...
double lognormal_quantile_c(double x, double m, double s)
Inverse ( ) of the cumulative distribution function of the upper tail of the lognormal distribution (...
const Double_t sigma
Double_t y[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
double igami(double a, double y)
double incbi(double a, double b, double y)
double ndtri(double y)
double log1p(double x)
declarations for functions which are not implemented by some compilers
Definition: Math.h:100
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static constexpr double s
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12