[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/mathutil.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2005 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.5.0, Dec 07 2006 ) */ 00008 /* The VIGRA Website is */ 00009 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00010 /* Please direct questions, bug reports, and contributions to */ 00011 /* koethe@informatik.uni-hamburg.de or */ 00012 /* vigra@kogs1.informatik.uni-hamburg.de */ 00013 /* */ 00014 /* Permission is hereby granted, free of charge, to any person */ 00015 /* obtaining a copy of this software and associated documentation */ 00016 /* files (the "Software"), to deal in the Software without */ 00017 /* restriction, including without limitation the rights to use, */ 00018 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00019 /* sell copies of the Software, and to permit persons to whom the */ 00020 /* Software is furnished to do so, subject to the following */ 00021 /* conditions: */ 00022 /* */ 00023 /* The above copyright notice and this permission notice shall be */ 00024 /* included in all copies or substantial portions of the */ 00025 /* Software. */ 00026 /* */ 00027 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00028 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00029 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00030 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00031 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00032 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00033 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00034 /* OTHER DEALINGS IN THE SOFTWARE. */ 00035 /* */ 00036 /************************************************************************/ 00037 00038 #ifndef VIGRA_MATHUTIL_HXX 00039 #define VIGRA_MATHUTIL_HXX 00040 00041 #include <cmath> 00042 #include <cstdlib> 00043 #include "config.hxx" 00044 #include "error.hxx" 00045 #include "tuple.hxx" 00046 #include "sized_int.hxx" 00047 #include "numerictraits.hxx" 00048 00049 /*! \page MathConstants Mathematical Constants 00050 00051 <TT>M_PI, M_SQRT2</TT> 00052 00053 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>" 00054 00055 Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized, 00056 we provide definitions here for those compilers that don't support them. 00057 00058 \code 00059 #ifndef M_PI 00060 # define M_PI 3.14159265358979323846 00061 #endif 00062 00063 #ifndef M_SQRT2 00064 # define M_SQRT2 1.41421356237309504880 00065 #endif 00066 \endcode 00067 */ 00068 #ifndef M_PI 00069 # define M_PI 3.14159265358979323846 00070 #endif 00071 00072 #ifndef M_SQRT2 00073 # define M_SQRT2 1.41421356237309504880 00074 #endif 00075 00076 namespace vigra { 00077 00078 /** \addtogroup MathFunctions Mathematical Functions 00079 00080 Useful mathematical functions and functors. 00081 */ 00082 //@{ 00083 // import functions into namespace vigra which VIGRA is going to overload 00084 00085 using VIGRA_CSTD::pow; 00086 using VIGRA_CSTD::floor; 00087 using VIGRA_CSTD::ceil; 00088 00089 // import abs(float), abs(double), abs(long double) from <cmath> 00090 // and abs(int), abs(long), abs(long long) from <cstdlib> 00091 using std::abs; 00092 00093 // define the missing variants of abs() to avoid 'ambigous overload' 00094 // errors in template functions 00095 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \ 00096 inline T abs(T t) { return t; } 00097 00098 VIGRA_DEFINE_UNSIGNED_ABS(bool) 00099 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char) 00100 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short) 00101 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int) 00102 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long) 00103 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long long) 00104 00105 #undef VIGRA_DEFINE_UNSIGNED_ABS 00106 00107 #define VIGRA_DEFINE_MISSING_ABS(T) \ 00108 inline T abs(T t) { return t < 0 ? -t : t; } 00109 00110 VIGRA_DEFINE_MISSING_ABS(signed char) 00111 VIGRA_DEFINE_MISSING_ABS(signed short) 00112 00113 #undef VIGRA_DEFINE_MISSING_ABS 00114 00115 /*! The rounding function. 00116 00117 Defined for all floating point types. Rounds towards the nearest integer 00118 such that <tt>abs(round(t)) == round(abs(t))</tt> for all <tt>t</tt>. 00119 00120 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00121 Namespace: vigra 00122 */ 00123 inline float round(float t) 00124 { 00125 return t >= 0.0 00126 ? floor(t + 0.5) 00127 : ceil(t - 0.5); 00128 } 00129 00130 inline double round(double t) 00131 { 00132 return t >= 0.0 00133 ? floor(t + 0.5) 00134 : ceil(t - 0.5); 00135 } 00136 00137 inline long double round(long double t) 00138 { 00139 return t >= 0.0 00140 ? floor(t + 0.5) 00141 : ceil(t - 0.5); 00142 } 00143 00144 /*! The square function. 00145 00146 sq(x) is needed so often that it makes sense to define it as a function. 00147 00148 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00149 Namespace: vigra 00150 */ 00151 template <class T> 00152 inline 00153 typename NumericTraits<T>::Promote sq(T t) 00154 { 00155 return t*t; 00156 } 00157 00158 namespace detail { 00159 00160 template <class T> 00161 class IntSquareRoot 00162 { 00163 public: 00164 static int sqq_table[]; 00165 static UInt32 exec(UInt32 v); 00166 }; 00167 00168 template <class T> 00169 int IntSquareRoot<T>::sqq_table[] = { 00170 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 00171 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 00172 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 00173 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118, 00174 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132, 00175 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145, 00176 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157, 00177 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168, 00178 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 00179 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 00180 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 00181 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206, 00182 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215, 00183 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223, 00184 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231, 00185 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, 00186 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 00187 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 00188 253, 254, 254, 255 00189 }; 00190 00191 template <class T> 00192 UInt32 IntSquareRoot<T>::exec(UInt32 x) 00193 { 00194 UInt32 xn; 00195 if (x >= 0x10000) 00196 if (x >= 0x1000000) 00197 if (x >= 0x10000000) 00198 if (x >= 0x40000000) { 00199 if (x >= (UInt32)65535*(UInt32)65535) 00200 return 65535; 00201 xn = sqq_table[x>>24] << 8; 00202 } else 00203 xn = sqq_table[x>>22] << 7; 00204 else 00205 if (x >= 0x4000000) 00206 xn = sqq_table[x>>20] << 6; 00207 else 00208 xn = sqq_table[x>>18] << 5; 00209 else { 00210 if (x >= 0x100000) 00211 if (x >= 0x400000) 00212 xn = sqq_table[x>>16] << 4; 00213 else 00214 xn = sqq_table[x>>14] << 3; 00215 else 00216 if (x >= 0x40000) 00217 xn = sqq_table[x>>12] << 2; 00218 else 00219 xn = sqq_table[x>>10] << 1; 00220 00221 goto nr1; 00222 } 00223 else 00224 if (x >= 0x100) { 00225 if (x >= 0x1000) 00226 if (x >= 0x4000) 00227 xn = (sqq_table[x>>8] >> 0) + 1; 00228 else 00229 xn = (sqq_table[x>>6] >> 1) + 1; 00230 else 00231 if (x >= 0x400) 00232 xn = (sqq_table[x>>4] >> 2) + 1; 00233 else 00234 xn = (sqq_table[x>>2] >> 3) + 1; 00235 00236 goto adj; 00237 } else 00238 return sqq_table[x] >> 4; 00239 00240 /* Run two iterations of the standard convergence formula */ 00241 00242 xn = (xn + 1 + x / xn) / 2; 00243 nr1: 00244 xn = (xn + 1 + x / xn) / 2; 00245 adj: 00246 00247 if (xn * xn > x) /* Correct rounding if necessary */ 00248 xn--; 00249 00250 return xn; 00251 } 00252 00253 } // namespace detail 00254 00255 using VIGRA_CSTD::sqrt; 00256 00257 /*! Signed integer square root. 00258 00259 Useful for fast fixed-point computations. 00260 00261 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00262 Namespace: vigra 00263 */ 00264 inline Int32 sqrti(Int32 v) 00265 { 00266 if(v < 0) 00267 throw std::domain_error("sqrti(Int32): negative argument."); 00268 return (Int32)detail::IntSquareRoot<UInt32>::exec((UInt32)v); 00269 } 00270 00271 /*! Unsigned integer square root. 00272 00273 Useful for fast fixed-point computations. 00274 00275 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00276 Namespace: vigra 00277 */ 00278 inline UInt32 sqrti(UInt32 v) 00279 { 00280 return detail::IntSquareRoot<UInt32>::exec(v); 00281 } 00282 00283 #ifndef VIGRA_HAS_HYPOT 00284 /*! Compute the Euclidean distance (length of the hypothenuse of a right-angled triangle). 00285 00286 The hypot() function returns the sqrt(a*a + b*b). 00287 It is implemented in a way that minimizes round-off error. 00288 00289 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00290 Namespace: vigra 00291 */ 00292 inline double hypot(double a, double b) 00293 { 00294 double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b); 00295 if (absa > absb) 00296 return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa)); 00297 else 00298 return absb == 0.0 00299 ? 0.0 00300 : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb)); 00301 } 00302 00303 #else 00304 00305 using ::hypot; 00306 00307 #endif 00308 00309 /*! The sign function. 00310 00311 Returns 1, 0, or -1 depending on the sign of \a t. 00312 00313 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00314 Namespace: vigra 00315 */ 00316 template <class T> 00317 T sign(T t) 00318 { 00319 return t > NumericTraits<T>::zero() 00320 ? NumericTraits<T>::one() 00321 : t < NumericTraits<T>::zero() 00322 ? -NumericTraits<T>::one() 00323 : NumericTraits<T>::zero(); 00324 } 00325 00326 /*! The binary sign function. 00327 00328 Transfers the sign of \a t2 to \a t1. 00329 00330 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00331 Namespace: vigra 00332 */ 00333 template <class T1, class T2> 00334 T1 sign(T1 t1, T2 t2) 00335 { 00336 return t2 >= NumericTraits<T2>::zero() 00337 ? abs(t1) 00338 : -abs(t1); 00339 } 00340 00341 #define VIGRA_DEFINE_NORM(T) \ 00342 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \ 00343 inline NormTraits<T>::NormType norm(T t) { return abs(t); } 00344 00345 VIGRA_DEFINE_NORM(bool) 00346 VIGRA_DEFINE_NORM(signed char) 00347 VIGRA_DEFINE_NORM(unsigned char) 00348 VIGRA_DEFINE_NORM(short) 00349 VIGRA_DEFINE_NORM(unsigned short) 00350 VIGRA_DEFINE_NORM(int) 00351 VIGRA_DEFINE_NORM(unsigned int) 00352 VIGRA_DEFINE_NORM(long) 00353 VIGRA_DEFINE_NORM(unsigned long) 00354 VIGRA_DEFINE_NORM(float) 00355 VIGRA_DEFINE_NORM(double) 00356 VIGRA_DEFINE_NORM(long double) 00357 00358 #undef VIGRA_DEFINE_NORM 00359 00360 template <class T> 00361 inline typename NormTraits<std::complex<T> >::SquaredNormType 00362 squaredNorm(std::complex<T> const & t) 00363 { 00364 return sq(t.real()) + sq(t.imag()); 00365 } 00366 00367 #ifdef DOXYGEN // only for documentation 00368 /*! The squared norm of a numerical object. 00369 00370 For scalar types: equals <tt>vigra::sq(t)</tt><br>. 00371 For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>. 00372 For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>. 00373 For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements). 00374 */ 00375 NormTraits<T>::SquaredNormType squaredNorm(T const & t); 00376 00377 #endif 00378 00379 /*! The norm of a numerical object. 00380 00381 For scalar types: implemented as <tt>abs(t)</tt><br> 00382 otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>. 00383 00384 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00385 Namespace: vigra 00386 */ 00387 template <class T> 00388 inline typename NormTraits<T>::NormType 00389 norm(T const & t) 00390 { 00391 typedef typename NormTraits<T>::SquaredNormType SNT; 00392 return sqrt(static_cast<typename SquareRootTraits<SNT>::SquareRootArgument>(squaredNorm(t))); 00393 } 00394 00395 namespace detail { 00396 00397 template <class T> 00398 T ellipticRD(T x, T y, T z) 00399 { 00400 double f = 1.0, s = 0.0, X, Y, Z, m; 00401 while(true) 00402 { 00403 m = (x + y + 3.0*z) / 5.0; 00404 X = 1.0 - x/m; 00405 Y = 1.0 - y/m; 00406 Z = 1.0 - z/m; 00407 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01) 00408 break; 00409 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z); 00410 s += f / (VIGRA_CSTD::sqrt(z)*(z + l)); 00411 f /= 4.0; 00412 x = (x + l)/4.0; 00413 y = (y + l)/4.0; 00414 z = (z + l)/4.0; 00415 } 00416 double a = X*Y; 00417 double b = sq(Z); 00418 double c = a - b; 00419 double d = a - 6.0*b; 00420 double e = d + 2.0*c; 00421 return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0) 00422 +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5); 00423 } 00424 00425 template <class T> 00426 T ellipticRF(T x, T y, T z) 00427 { 00428 double X, Y, Z, m; 00429 while(true) 00430 { 00431 m = (x + y + z) / 3.0; 00432 X = 1.0 - x/m; 00433 Y = 1.0 - y/m; 00434 Z = 1.0 - z/m; 00435 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01) 00436 break; 00437 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z); 00438 x = (x + l)/4.0; 00439 y = (y + l)/4.0; 00440 z = (z + l)/4.0; 00441 } 00442 double d = X*Y - sq(Z); 00443 double p = X*Y*Z; 00444 return (1.0 - d/10.0 + p/14.0 + sq(d)/24.0 - d*p*3.0/44.0) / VIGRA_CSTD::sqrt(m); 00445 } 00446 00447 } // namespace detail 00448 00449 /*! The incomplete elliptic integral of the first kind. 00450 00451 Computes 00452 00453 \f[ 00454 \mbox{F}(x, k) = \int_0^x \frac{1}{\sqrt{1 - k^2 \sin(t)^2}} dt 00455 \f] 00456 00457 according to the algorithm given in Press et al. "Numerical Recipes". 00458 00459 Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral 00460 functions must be k^2 rather than k. Check the documentation when results disagree! 00461 00462 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00463 Namespace: vigra 00464 */ 00465 inline double ellipticIntegralF(double x, double k) 00466 { 00467 double c2 = sq(VIGRA_CSTD::cos(x)); 00468 double s = VIGRA_CSTD::sin(x); 00469 return s*detail::ellipticRF(c2, 1.0 - sq(k*s), 1.0); 00470 } 00471 00472 /*! The incomplete elliptic integral of the second kind. 00473 00474 Computes 00475 00476 \f[ 00477 \mbox{E}(x, k) = \int_0^x \sqrt{1 - k^2 \sin(t)^2} dt 00478 \f] 00479 00480 according to the algorithm given in Press et al. "Numerical Recipes". The 00481 complete elliptic integral of the second kind is simply <tt>ellipticIntegralE(M_PI/2, k)</TT>. 00482 00483 Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral 00484 functions must be k^2 rather than k. Check the documentation when results disagree! 00485 00486 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00487 Namespace: vigra 00488 */ 00489 inline double ellipticIntegralE(double x, double k) 00490 { 00491 double c2 = sq(VIGRA_CSTD::cos(x)); 00492 double s = VIGRA_CSTD::sin(x); 00493 k = sq(k*s); 00494 return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0)); 00495 } 00496 00497 #ifndef VIGRA_HAS_ERF 00498 00499 namespace detail { 00500 00501 template <class T> 00502 double erfImpl(T x) 00503 { 00504 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x)); 00505 double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+ 00506 t*(0.09678418+t*(-0.18628806+t*(0.27886807+ 00507 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+ 00508 t*0.17087277))))))))); 00509 if (x >= 0.0) 00510 return 1.0 - ans; 00511 else 00512 return ans - 1.0; 00513 } 00514 00515 } // namespace detail 00516 00517 /*! The error function. 00518 00519 If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the 00520 new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error 00521 function 00522 00523 \f[ 00524 \mbox{erf}(x) = \int_0^x e^{-t^2} dt 00525 \f] 00526 00527 according to the formula given in Press et al. "Numerical Recipes". 00528 00529 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00530 Namespace: vigra 00531 */ 00532 inline double erf(double x) 00533 { 00534 return detail::erfImpl(x); 00535 } 00536 00537 #else 00538 00539 using VIGRA_CSTD::erf; 00540 00541 #endif 00542 00543 namespace detail { 00544 00545 template <class T> 00546 double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, T noncentrality, T arg) 00547 { 00548 double a = degreesOfFreedom + noncentrality; 00549 double b = (a + noncentrality) / sq(a); 00550 double t = (VIGRA_CSTD::pow((double)arg / a, 1.0/3.0) - (1.0 - 2.0 / 9.0 * b)) / VIGRA_CSTD::sqrt(2.0 / 9.0 * b); 00551 return 0.5*(1.0 + erf(t/VIGRA_CSTD::sqrt(2.0))); 00552 } 00553 00554 template <class T> 00555 void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans, unsigned int & j) 00556 { 00557 double tol = -50.0; 00558 if(lans < tol) 00559 { 00560 lans = lans + VIGRA_CSTD::log(arg / j); 00561 dans = VIGRA_CSTD::exp(lans); 00562 } 00563 else 00564 { 00565 dans = dans * arg / j; 00566 } 00567 pans = pans - dans; 00568 j += 2; 00569 } 00570 00571 template <class T> 00572 std::pair<double, double> noncentralChi2CDF(unsigned int degreesOfFreedom, T noncentrality, T arg, T eps) 00573 { 00574 vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0, 00575 "noncentralChi2P(): parameters must be positive."); 00576 if (arg == 0.0 && degreesOfFreedom > 0) 00577 return std::make_pair(0.0, 0.0); 00578 00579 // Determine initial values 00580 double b1 = 0.5 * noncentrality, 00581 ao = VIGRA_CSTD::exp(-b1), 00582 eps2 = eps / ao, 00583 lnrtpi2 = 0.22579135264473, 00584 probability, density, lans, dans, pans, sum, am, hold; 00585 unsigned int maxit = 500, 00586 i, m; 00587 if(degreesOfFreedom % 2) 00588 { 00589 i = 1; 00590 lans = -0.5 * (arg + VIGRA_CSTD::log(arg)) - lnrtpi2; 00591 dans = VIGRA_CSTD::exp(lans); 00592 pans = erf(VIGRA_CSTD::sqrt(arg/2.0)); 00593 } 00594 else 00595 { 00596 i = 2; 00597 lans = -0.5 * arg; 00598 dans = VIGRA_CSTD::exp(lans); 00599 pans = 1.0 - dans; 00600 } 00601 00602 // Evaluate first term 00603 if(degreesOfFreedom == 0) 00604 { 00605 m = 1; 00606 degreesOfFreedom = 2; 00607 am = b1; 00608 sum = 1.0 / ao - 1.0 - am; 00609 density = am * dans; 00610 probability = 1.0 + am * pans; 00611 } 00612 else 00613 { 00614 m = 0; 00615 degreesOfFreedom = degreesOfFreedom - 1; 00616 am = 1.0; 00617 sum = 1.0 / ao - 1.0; 00618 while(i < degreesOfFreedom) 00619 detail::noncentralChi2OneIteration(arg, lans, dans, pans, i); 00620 degreesOfFreedom = degreesOfFreedom + 1; 00621 density = dans; 00622 probability = pans; 00623 } 00624 // Evaluate successive terms of the expansion 00625 for(++m; m<maxit; ++m) 00626 { 00627 am = b1 * am / m; 00628 detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom); 00629 sum = sum - am; 00630 density = density + am * dans; 00631 hold = am * pans; 00632 probability = probability + hold; 00633 if((pans * sum < eps2) && (hold < eps2)) 00634 break; // converged 00635 } 00636 if(m == maxit) 00637 vigra_fail("noncentralChi2P(): no convergence."); 00638 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability))); 00639 } 00640 00641 } // namespace detail 00642 00643 /*! Chi square distribution. 00644 00645 Computes the density of a chi square distribution with \a degreesOfFreedom 00646 and tolerance \a accuracy at the given argument \a arg 00647 by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>. 00648 00649 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00650 Namespace: vigra 00651 */ 00652 inline double chi2(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7) 00653 { 00654 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).first; 00655 } 00656 00657 /*! Cumulative chi square distribution. 00658 00659 Computes the cumulative density of a chi square distribution with \a degreesOfFreedom 00660 and tolerance \a accuracy at the given argument \a arg, i.e. the probability that 00661 a random number drawn from the distribution is below \a arg 00662 by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>. 00663 00664 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00665 Namespace: vigra 00666 */ 00667 inline double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7) 00668 { 00669 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).second; 00670 } 00671 00672 /*! Non-central chi square distribution. 00673 00674 Computes the density of a non-central chi square distribution with \a degreesOfFreedom, 00675 noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument 00676 \a arg. It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from 00677 http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of 00678 degrees of freedom. 00679 00680 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00681 Namespace: vigra 00682 */ 00683 inline double noncentralChi2(unsigned int degreesOfFreedom, 00684 double noncentrality, double arg, double accuracy = 1e-7) 00685 { 00686 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).first; 00687 } 00688 00689 /*! Cumulative non-central chi square distribution. 00690 00691 Computes the cumulative density of a chi square distribution with \a degreesOfFreedom, 00692 noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument 00693 \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg 00694 It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from 00695 http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of 00696 degrees of freedom (see noncentralChi2CDFApprox() for a constant-time algorithm). 00697 00698 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00699 Namespace: vigra 00700 */ 00701 inline double noncentralChi2CDF(unsigned int degreesOfFreedom, 00702 double noncentrality, double arg, double accuracy = 1e-7) 00703 { 00704 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).second; 00705 } 00706 00707 /*! Cumulative non-central chi square distribution (approximate). 00708 00709 Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom, 00710 and noncentrality parameter \a noncentrality at the given argument 00711 \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg 00712 It uses the approximate transform into a normal distribution due to Wilson and Hilferty 00713 (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32). 00714 The algorithm's running time is independent of the inputs, i.e. is should be used 00715 when noncentralChi2CDF() is too slow, and approximate values are sufficient. The accuracy is only 00716 about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5. 00717 00718 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00719 Namespace: vigra 00720 */ 00721 inline double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg) 00722 { 00723 return detail::noncentralChi2CDFApprox(degreesOfFreedom, noncentrality, arg); 00724 } 00725 00726 00727 00728 namespace detail { 00729 00730 // both f1 and f2 are unsigned here 00731 template<class FPT> 00732 inline 00733 FPT safeFloatDivision( FPT f1, FPT f2 ) 00734 { 00735 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max() 00736 ? NumericTraits<FPT>::max() 00737 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) || 00738 f1 == NumericTraits<FPT>::zero() 00739 ? NumericTraits<FPT>::zero() 00740 : f1/f2; 00741 } 00742 00743 } // namespace detail 00744 00745 /*! Tolerance based floating-point comparison. 00746 00747 Check whether two floating point numbers are equal within the given tolerance. 00748 This is useful because floating point numbers that should be equal in theory are 00749 rarely exactly equal in practice. If the tolerance \a epsilon is not given, 00750 twice the machine epsilon is used. 00751 00752 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00753 Namespace: vigra 00754 */ 00755 template <class T1, class T2> 00756 bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon) 00757 { 00758 typedef typename PromoteTraits<T1, T2>::Promote T; 00759 if(l == 0.0) 00760 return VIGRA_CSTD::fabs(r) <= epsilon; 00761 if(r == 0.0) 00762 return VIGRA_CSTD::fabs(l) <= epsilon; 00763 T diff = VIGRA_CSTD::fabs( l - r ); 00764 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) ); 00765 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) ); 00766 00767 return (d1 <= epsilon && d2 <= epsilon); 00768 } 00769 00770 template <class T1, class T2> 00771 bool closeAtTolerance(T1 l, T2 r) 00772 { 00773 typedef typename PromoteTraits<T1, T2>::Promote T; 00774 return closeAtTolerance(l, r, 2.0 * NumericTraits<T>::epsilon()); 00775 } 00776 00777 //@} 00778 00779 } // namespace vigra 00780 00781 #endif /* VIGRA_MATHUTIL_HXX */
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|