libs/math/doc/sf/poisson_optimisation.qbk

   1 [sect:optim Optimisation Examples]
   2
   3 [h4 Poisson Distribution - Optimization and Accuracy is quite complicated.
   4
   5 The general formula for calculating the CDF uses the incomplete gamma thus:
   6
   7   return gamma_q(k+1, mean);
   8
   9 But the case of small integral k is *very* common, so it is worth considering optimisation.
  10
  11 The first obvious step is to use a finite sum of each PDF (Probability *density* function)
  12 for each value of k to build up the CDF (*cumulative* distribution function).
  13
  14 This could be done using the PDF function for the distribution,
  15 for which there are two equivalent formulae:
  16
  17   return exp(-mean + log(mean) * k - lgamma(k+1));
  18
  19   return gamma_p_derivative(k+1, mean);
  20
  21 The pdf would probably be more accurate using `gamma_p_derivative`.
  22
  23 The reason is that the expression:
  24
  25   -mean + log(mean) * k - lgamma(k+1)
  26
  27 Will produce a value much smaller than the largest of the terms, so you get
  28 cancellation error: and then when you pass the result to `exp()` which
  29 converts the absolute error in its argument to a relative error in the
  30 result (explanation available if required), you effectively amplify the
  31 error further still.
  32
  33 `gamma_p_derivative` is just a thin wrapper around some of the internals of
  34 the incomplete gamma, it does its utmost to avoid issues like this, because
  35 this function is responsible for virtually all of the error in the result.
  36 Hopefully further advances in the future might improve things even further
  37 (what is really needed is an 'accurate' `pow(1+x)` function, but that's a whole
  38 other story!).
  39
  40 But calling `pdf` function makes repeated, redundant, checks on the value of `mean` and `k`,
  41
  42   result += pdf(dist, i);
  43
  44 so it may be faster to substitute the formula for the pdf in a summation loop
  45
  46   result += exp(-mean) * pow(mean, i) / unchecked_factorial(i);
  47
  48 (simplified by removing casting from RealType).
  49
  50 Of course, mean is unchanged during this summation,
  51 so exp(mean) should only be calculated once, outside the loop.
  52
  53 Optimising compilers 'might' do this, but one can easily ensure this.
  54
  55 Obviously too, k must be small enough that unchecked_factorial is OK:
  56 34 is an obvious choice as the limit for 32-bit float.
  57 For larger k, the number of iterations is like to be uneconomic.
  58 Only experiment can determine the optimum value of k
  59 for any particular RealType (float, double...)
  60
  61 But also note that
  62
  63 The incomplete gamma already optimises this case
  64 (argument "a" is a small int),
  65 although only when the result q (1-p) would be < 0.5.
  66
  67 And moreover, in the above series, each term can be calculated
  68 from the previous one much more efficiently:
  69
  70   cdf = sum from 0 to k of C[k]
  71
  72 with:
  73
  74   C[0] = exp(-mean)
  75   C[N+1] = C[N] * mean / (N+1)
  76
  77 So hopefully that's:
  78
  79      {
  80        RealType result = exp(-mean);
  81        RealType term = result;
  82        for(int i = 1; i <= k; ++i)
  83        { // cdf is sum of pdfs.
  84           term *= mean / i;
  85           result += term;
  86        }
  87        return result;
  88      }
  89
  90 This is exactly the same finite sum as used by `gamma_p/gamma_q` internally.
  91
  92 As explained previously it's only used when the result
  93
  94   p > 0.5 or 1-p = q < 0.5.
  95
  96 The slight danger when using the sum directly like this, is that if
  97 the mean is small and k is large then you're calculating a value ~1, so
  98 conceivably you might overshoot slightly.  For this and other reasons in the
  99 case when p < 0.5 and q > 0.5 `gamma_p/gamma_q` use a different (infinite but
 100 rapidly converging) sum, so that danger isn't present since you always
 101 calculate the smaller of p and q.
 102
 103 So... it's tempting to suggest that you just call `gamma_p/gamma_q` as
 104 required.  However, there is a slight benefit for the k = 0 case because you
 105 avoid all the internal logic inside `gamma_p/gamma_q` trying to figure out
 106 which method to use etc.
 107
 108 For the incomplete beta function, there are no simple finite sums
 109 available (that I know of yet anyway), so when there's a choice between a
 110 finite sum of the PDF and an incomplete beta call, the finite sum may indeed
 111 win out in that case.
 112
 113 [endsect] [/sect:optim Optimisation Examples]
 114
 115 [/
 116   Copyright 2006 John Maddock and Paul A. Bristow.
 117   Distributed under the Boost Software License, Version 1.0.
 118   (See accompanying file LICENSE_1_0.txt or copy at
 119   http://www.boost.org/LICENSE_1_0.txt).
 120 ]