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